Azzera filtri
Azzera filtri

A weird feature in the obtained function from calculating the inverse Fourier transform using ifft

26 visualizzazioni (ultimi 30 giorni)
Hey,
I am calculating the inverse Fourier transform of a function and obtaining the result in the time domain using the ifft function in matlab. When I perform the Fourier transform and plot the absolute value of the obtained function in time I see a little figure at the other side of the time axis like a tail for the function. I am wondering does it mean there is an issue in the inverse Fourier transform taken?
I have attached the pictures for the whole function I obtain in the time domain using the ifft function and the two chunks of the function at the two ends of the time axis. My question is about the small chunk on the right of the whole function's picture shown.
I appreciate your comments.
  4 Commenti

Accedi per commentare.

Risposte (2)

Matt J
Matt J il 15 Apr 2023
Modificato: Matt J il 15 Apr 2023
IFFTs are N-periodic, so if you shift a lobe far enough to the left or right of a length N sampling, one side of the lobe will start to wrap around circulantly to the other side.
  9 Commenti
Shaily_T
Shaily_T il 16 Apr 2023
Thanks @Image Analyst! I would say if I have only the output on the left of the whole obtained function in time and not the small chunk on the right it is reasonable and as expected. However, I do have that small chunk on the right of the obtained function in time when I calculate the inverse Fourier transform using the ifft() function. No, I don't use ifftshift. And the plot is for the absolute value of the the inverse fourier transform result using ifft() versus time.
Shaily_T
Shaily_T il 16 Apr 2023
@Matt J Thanks! So, I'll play around with the frequency sampling interval to obtain the desirable time horizon of width T to see what is the behaviour of the small chunk on the right of the whole function image.

Accedi per commentare.


Paul
Paul il 16 Apr 2023
Hi Shaily_T,
It looks like the purpose of the code is to linearly convolve an input signal with an impulse response via mulitiplication in the frequency domain. However, multiplying the DFT's of two signals is equivalent to a circular convolution. Linear convolution requires zero padding the input and the impulse response in the time domain before taking their DFT's.
Are h and FourierInputField developed from taking the FFT with zero padding?
  9 Commenti
Shaily_T
Shaily_T il 19 Apr 2023
@Paul @Walter Roberson It is my fault. I deleted the .png files and the code as I have put the exact same images and file in my new question here https://www.mathworks.com/matlabcentral/answers/1949413-the-difference-between-the-obtained-result-using-fftshift-when-employing-the-ifft-function?s_tid=mlc_ans_email_ques. But it seems it was not an issue to have them in both questions. So, I reattached them. Sorry for the confusion.
Paul
Paul il 22 Apr 2023
Without knowing much about FourierInputField, h, or nu .....
Load and plot the data
fnames = unzip('ifft.zip');
load(fnames{2},'nu','h');
load(fnames{4},'FourierInputField');
figure
plot(subplot(411),nu,real(FourierInputField))
plot(subplot(412),nu,imag(FourierInputField))
plot(subplot(413),nu,real(h))
plot(subplot(414),nu,imag(h))
We see that the data is not symmetric around 0, so I'll chop the data off on the right to make the data symmetric around zero.
ii = find(nu > -nu(1));
FourierInputField(ii) = [];
nu(ii) = [];
h(ii) = [];
Plot the truncated data.
figure
plot(subplot(411),nu,real(FourierInputField))
plot(subplot(412),nu,imag(FourierInputField))
plot(subplot(413),nu,real(h))
plot(subplot(414),nu,imag(h))
It was stated in this comment that FourierInputField are samples of the CTFT of a Gaussian pulse. If this is the case, FourierInputField should have linear phase
figure
plot(nu,unwrap(angle(FourierInputField))),grid
The phase looks linear, but it should go through zero at nu = 0. I think it doesn't in this plot because of unwrapping confusion, but it really does
figure
plot(nu,angle(FourierInputField)),grid
xlim([-.1 .1]*1e8)
Because FourierInputField is on both sides of zero, we first iffshift, then take the ifft, and then fftshift to get the time domain result centered on n = 0
TimeInputField = fftshift(ifft(ifftshift(FourierInputField)));
Now make a discrete-time vector. No element in the nu vector is equal to zero, which seems odd
min(abs(nu))
ans = 1.3526e+04
so we don't have the typical case where we'd have nu = n*Fs; Instead I'll approximate the time vector as
N = numel(nu)
N = 8683
Fs = mean(diff(nu))*N;
Ts = 1/Fs;
t = (-8682/2:8682/2)*Ts;
Plot TimeInputField
figure
plot(t,abs(TimeInputField))
xline(34*Ts)
xlim([0 50]*Ts)
It's basically a pulse centered on t = 34*Ts. If it's a Gaussian pulse, the phase should be zero
figure
plot(t,angle(TimeInputField))
xlim([0 50]*Ts)
The phase is zero for values of t where where the pulse amplitude is not insignificant.
Now compute the output, which is the circular convolution of the TimeInputField and Timeh, whatever that is. Even though h those sharp discontinuites and flat spots, they occur out on the tails of FourierInputField, so maybe those features of h don't matter.
y = fftshift(ifft(ifftshift(FourierInputField.*h)));
figure
plot(t,abs(y)),grid
xlim([-50 100]*Ts)
y does have some peaking in negative time. Those peaks whould show up on the far right of the result if we hadn't applied fftshift
ynoshift = (ifft(ifftshift(FourierInputField.*h)));
figure
plot((0:numel(ynoshift)-1)*Ts,abs(ynoshift))
xlim([numel(ynoshift)-50, numel(ynoshift)]*Ts)
As done in the code provided with the question, we can zero pad in the frequency domain to interpolate the time domain result. When plotting we have to scale Ts to account for the additional samples added in the frequency domain.
y = fftshift(ifft(ifftshift([zeros(1,4e4),FourierInputField, zeros(1,4e4)].*[zeros(1,4e4), h, zeros(1,4e4)])));
Tsinterp = Ts*N/numel(y);
figure
t = (-88682/2:88682/2)*Tsinterp;
plot(t,abs(y)),grid
xline(34*Ts)
xlim([-500 1000]*Tsinterp)

Accedi per commentare.

Prodotti


Release

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by