A weird feature in the obtained function from calculating the inverse Fourier transform using ifft
38 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(2 个)
Paul
2023-4-16
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 个评论
Paul
2023-4-22
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))
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)
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)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!