FFT convolution and iFFT process

48 views (last 30 days)
Yassine
Yassine on 2 Aug 2022
Answered: Paul on 2 Aug 2022
I'm having a problem with conv, fft and ifft MATLAB functions.
So basically I have a transfer function FT=X/C in frequency domain of size N_FT and a step df (can be variable), and an input signal C(t) of size N and with a time step 1/Fs (given) and want to get the output X(t) in time domain.
The maneuvre is rather simple : C(t) to frequency domain via FFT to obtain C(f) ==> Convolute C(f) and FT(f) to obtain X(f) ==> Transfer X(f) back to time domain via iFFT. SO X(t)=ifft(conv(FT,fft(C(t)))
Constructing the frequency vector respecting the rule f=step*(0:L-1)/L aka start at 0, have a length of L=N=length(C(t)) and a step=1/dt=Fs/N
Below is a simplified version of my code :
f=Fs*(0:N-1)/N; %making the frequency vector for C(f)
Cf=real(fft(Ct,N)); % getting C(f)
Xf=conv(Cf,FT)*df; %convoluting and multiplying by df because I saw this somewhere else
Xt=real(ifft(Xf)); %getting Xt
plot(t,Xt) %plotting Xt with the time vector of C(t) (because logically they have to be the same)
My problem is, the output X(t) curve depends on the size of C(t)=N and doesnt look much like the expected result. (see screenshot 1)
HOWEVER, when I try to set another N_fft to assure I have a step of df equivalent of that to FT and a I get better result (see screenshot 2)
N_fft=round(4*Fs/df); %making Cf vector have the same step as FT vector // multiplying by 4 makes it the closest in this case
f=linspace(0,Fs,N_fft+1); %making the frequency vector for C(f)
Cf=real(fft(Ct,N_fft)); % getting C(f)
Xf=conv(Cf,FT)*df; %multiplying by df because I saw this somewhere else
Xt=real(ifft(Xf)); %getting Xt
plot(t,Xt)
I'm wondering how to choose N_fft correctly to get the right signal.
I appreciate any advice.
  3 Comments
Yassine
Yassine on 2 Aug 2022
So I don't have to worry about FT(f) and C(f) not having the same frequency step ?
However, choosing N_fft like you said gives the same result as the first one :(
Star Strider
Star Strider on 2 Aug 2022
The sampling frequency of the filter design must be the same as the sampling frequency of the original signal (and the sampling intervals must be constant).
Using the ‘NFFT’ definition I suggested may not change the results, however it will increase the efficiency of the fft calculation, and its inverse.

Sign in to comment.

Answers (1)

Paul
Paul on 2 Aug 2022
Hi Yassine,
What does this mean: "df (can be variable),"? df should be a constant for all of FT. Is that not the case?
Let's assume that all signals are sampled in time with the same sampling period (Ts) or sampling frequency Fs = 1/Ts.
In order to perform linear convolution in time via multiplication of DFTs, we first have zero pad in the time domain. So starting in the time domain for now, the sequence is: zeropad -> FFT -> muliply -> IFFT.
For example, let ft be the time domain samples corresponding to FT, i.e., the finite duration impulse response of transfer function
rng(100);
ft = rand(1,5);
Let ct be the time domain input sequence
ct = rand(1,7);
Then, we know the output should be
xt = conv(ft,ct)
xt = 1×11
0.0661 0.3983 0.6871 0.6916 1.2684 1.4033 0.7253 0.9231 0.8445 0.1809 0.0010
We can also obtain this output as
N = numel(ct) + numel(ft) - 1;
ifft(fft(zeropad(ft,N)).*fft(zeropad(ct,N)))
ans = 1×11
0.0661 0.3983 0.6871 0.6916 1.2684 1.4033 0.7253 0.9231 0.8445 0.1809 0.0010
Now, in this Question, if I understand it correctly, we start with ct and FT, where FT is
FT = fft(ft); % no zero padding
So if we want to get xt, we need first get ft, then zero pad, and then proceed as above
N = numel(FT) + numel(ct) - 1;
x = ifft(fft(zeropad(ifft(FT),N)).*fft(zeropad(ct,N)))
x = 1×11
0.0661 0.3983 0.6871 0.6916 1.2684 1.4033 0.7253 0.9231 0.8445 0.1809 0.0010
Finally, if we want the samples of x to approximate samples of the continuous-time convolution, we neet to multilply by Ts
x = x*Ts;
function out = zeropad(in,N)
out = [in zeros(1,N-numel(in))]; % assumes in is a row vector
end

Categories

Find more on Fourier Analysis and Filtering in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!