CODEWITHSUNDEEP
matlabfft
Medulla Oblongata
Medulla Oblongata

Reputation: 3961

inverse fourier transform of data does not give correct amplitude

I'm trying to calculate the inverse Fourier transform of some data using Matlab. I start with raw data in the frequency domain and want to visualise the data in the time domain. Here is my MWE:

a = 1.056;

% frequency of data (I cannot change this)
w = linspace(-100,100,1e6);

L = length(w); % no. sample points
ts = L/1000; % time sampling
Ts = ts/L; % sampling rate
Fs = 1/Ts; % sampling freq

t = (-L/2:L/2-1)*ts/L; % time

Y = sqrt(pi/a)*exp(-w.^2/(4*a)); % my data

yn = Fs*ifftshift(ifft(fftshift(Y(end:-1:1)))) % numerical soln
ya = exp(-a*t.^2); % analytic solution

figure; hold on
plot(t,yn,'.')
plot(t,ya,'-')
xlabel('time, t')
legend('numerical','analytic')
xlim([-5,5])

I have adapted the code from this question however the amplitude is too large:

enter image description here

Can you please tell me what I'm doing wrong?

Upvotes: 1

Views: 582

Answers (1)

Cris Luengo
Cris Luengo

Reputation: 60504

There are three issues with your code:

  1. You define ts = L/1000 and then compute Fs, which gives you 1000. But Fs is given by the w array you've set up: the full range of w is 2*pi*Fs:

    Fs = -w(1)/pi;  % sampling freq
Ts = 1/Fs;      % sampling rate

    Or, equivalently, set Fs = mean(diff(w))*L / (2*pi)

  2. w is defined, but does not include 0. Just like you define t carefully to include the 0 in just the right place, so should you define w to include 0 in just the right place. One simple way to do this is to define it with one more value, then delete the last value:

    w = linspace(-100,100,1e6+1);
w(end) = [];

    If your input data does not include the 0 frequency, you should resample it so that it does. The DFT (FFT) expects a 0 frequency bin.

  3. You're using ifftshift and fftshift reversed: fftshift shifts the origin from the leftmost array element to the middle, and ifftshift shifts it from the middle to the left. You define your signal with the origin in the middle, so you need to use ifftshift on it to move the origin where the fft and ifft functions expect it. Use fftshift on the output of these two functions to center the origin for display. Because your data is even sized, these two functions do exactly the same thing, and you will not notice the difference. But if the data were odd sized, you'd see the difference.

The following code gives a perfect match:

a = 1.056;

% frequency of data (I cannot change this)
w = linspace(-100,100,1e6+1); w(end) = [];

L = length(w); % no. sample points
Fs = -w(1)/pi; % sampling freq
Ts = 1/Fs; % sampling rate

t = (-L/2:L/2-1)*Ts; % time

Y = sqrt(pi/a)*exp(-w.^2/(4*a)); % my data

yn = Fs*fftshift(ifft(ifftshift(Y(end:-1:1)))); % numerical soln
ya = exp(-a*t.^2); % analytic solution

figure; hold on
plot(t,yn,'.')
plot(t,ya,'-')
xlabel('time, t')
legend('numerical','analytic')
xlim([-5,5])

Upvotes: 1

Related Questions