# MATLAB IFFT doesn't match the analytical one

• MATLAB
tworitdash
I have a Gaussian shape frequency domain spectrum of which I am calculating the Inverse Fourier transform. I use both IFFT of MATLAB and also an analytical expression of Mathematica. They are not the same.

I don't know where it is wrong. I have pasted both the figures. The numerical one has been zoomed in actually because the time domain pulse came out to be too thin. Any help is appreciated.

The code is below:
MATLAB code for IFFT with Analytical equation:
clear;
% close all;

lambda = 0.03;

PRT = 1e-3;
mu = 4; % Mean Doppler
sigma = 0.2;
v_amb = 7.5;

N = 60000; % Total number of points in time axis

t1 = 0:PRT:(N - 1)*PRT; % Time axis
vel_axis = linspace(-v_amb, v_amb, N); % velocity axis for the entire rotation

% vel_axis_hs = linspace(-v_amb, v_amb, hs); % velocity axis for one beamwidth
% t1 = -N/2*PRT:PRT:(N/2 - 1)*PRT; % Time axis
% s_analyt_o = 1./sqrt(2*pi*(sigma).^2) .* exp(-(vel_axis - mu).^2./(2*(sigma).^2));

% SNR = 10^(30/10);
% X = rand(1, N);
% Theta = 2 .* pi * rand(1, N);

S_ = gaussmf(vel_axis, [sigma, mu]);

% Noise = sum(S_) ./ (N .* SNR);
% s_analyt_o = -(S_ + Noise) .* log(X);

s_analyt_o = S_;
s_num = ifft(ifftshift(sqrt(s_analyt_o))); % Numerical IFFT

s_analyt = ifftshift((2/pi)^(1/4) * sqrt(sigma) .* exp(-t1 .* ((sigma)^2 .* t1 + 1j .* mu))); % analytical IFFT

figure; plot(t1, abs((s_analyt)));

figure; plot(t1, abs(ifftshift(s_num)));

Delta2

Gold Member
This kind of problem depends on the details of what you are doing. If you provide a more clear description of exactly what you are trying to do, including equations and how you used Mathematica to compute the IFFT, it would be much easier to help. It is easy to have missing factors of ##2 \pi##, ##N##, etc, in amplitudes or function arguments. And I suspect that you used Mathematica to comput the inverse Fourier transform, not the inverse fast Fourier transform (IFFT) - yes?
Jason

hutchphd and Twigg
Mentor
Also it would be nice if the diagrams were labeled so we know if they came from MATLAB or Mathematica.

tworitdash
tworitdash
This kind of problem depends on the details of what you are doing. If you provide a more clear description of exactly what you are trying to do, including equations and how you used Mathematica to compute the IFFT, it would be much easier to help. It is easy to have missing factors of ##2 \pi##, ##N##, etc, in amplitudes or function arguments. And I suspect that you used Mathematica to comput the inverse Fourier transform, not the inverse fast Fourier transform (IFFT) - yes?
Jason
I want to update some things in the question. I will make it clear, but I can't see the edit button anymore. Should I update the things in a new answer?

Homework Helper
Gold Member
I want to update some things in the question. I will make it clear, but I can't see the edit button anymore. Should I update the things in a new answer?
I think the edit option is lost in 24 hours after you initially post the message.
So
Either message a moderator (e.g. @Greg Bernhardt ) to restore the edit option (if that is possible, i am not sure)
Or just write a new post here.

tworitdash
Gold Member
Shouldn’t you ifftshift after transforming, not before?

tworitdash
tworitdash
Now, I will write a more detailed explanation of what I was trying to do. I have considered a Gaussian-shaped spectrum for Doppler velocities that follows the following equation,

$$S(v) = \frac{1}{\sqrt(2\pi \sigma^2)} e^{-\frac{(v - \mu)^2}{2\sigma^2}}$$

Here,

It is assumed to be a group of scatters of different shapes and sizes moving with different velocities.

I want to see how does it look in the time domain when the radar is receiving this signal in slow time. As this is a power spectrum, for the time-domain we need the IDFT of the square root of this power spectrum (voltage equivalent).

The pulse repetition time is 1 [ms], which corresponds to a pulse repetition frequency of 1 [kHz].

The velocity axis is,

$$v = [-\frac{\lambda}{4 \delta t}, \frac{\lambda}{4 \delta t})$$ with steps of $\frac{\lambda}{2 N \delta t}$.

The N is the number of pulses, the $\delta t$ is the pulse repetition time, $\lambda$ is the operating wavelength.

I approached it in two different ways. First, in MATLAB numerically using ifft function. Second, using MATLAB but analytical expression obtained from WolframAlpha (The mathematica engine). The screenshot is attached from Wolfram Alpha. As, in WolframAlpha it just converted from frequency to time, I used some modifications which may make it a velocity-time conversion.

The expression I used is this,

$$s(t) = (\frac{2}{\pi})^{1/4} \sqrt{\frac{2\sigma}{\lambda}} e^{-t ( (\sigma \frac{4\pi}{\lambda})^2 t + j \mu \frac{4\pi}{\lambda} )}$$

The code is shown below.

IDFT:
clear;
% close all;

lambda = 0.03;

PRT = 1e-3;
mu = 4; % Mean Doppler
sigma = 0.2;
v_amb = 7.5;

N = 60000; % Total number of points in time axis

t1 = 0:PRT:(N - 1)*PRT; % Time axis
vel_axis = linspace(-v_amb, v_amb, N); % velocity axis for the entire rotation

% vel_axis_hs = linspace(-v_amb, v_amb, hs); % velocity axis for one beamwidth
% t1 = -N/2*PRT:PRT:(N/2 - 1)*PRT; % Time axis
% s_analyt_o = 1./sqrt(2*pi*(sigma).^2) .* exp(-(vel_axis - mu).^2./(2*(sigma).^2));

% SNR = 10^(30/10);
% X = rand(1, N);
% Theta = 2 .* pi * rand(1, N);

S_ = gaussmf(vel_axis, [sigma, mu]);

% Noise = sum(S_) ./ (N .* SNR);
% s_analyt_o = -(S_ + Noise) .* log(X);

s_analyt_o = S_;
s_num = ifft(ifftshift(sqrt(s_analyt_o))); % Numerical IFFT

s_analyt = ((2/pi)^(1/4) * sqrt(sigma * 4 * pi / lambda) .* exp(-t1 .* ((sigma * 4 * pi / lambda)^2 .* t1 + 1j .* mu .* 4 * pi / lambda))); % analytical IFFT

figure; plot(t1, abs((s_analyt))/max(abs(s_analyt)), 'LineWidth', 2);

hold on; plot(t1, abs((s_num))/max(abs(s_num)), '*');
grid on;

xlim([-10 80]);

legend({'Wolfram Alpha Expression', 'MATLAB numerical'})

The results are shown. Here, I don't use any ifftshift. I have put one zoomed in image near time 0 and one at time 60 [sec]. These are the extreme time values. In the MATLAB one we see at both edges we have a pulse. However, we don't have it in the analytical one. It has the pulse only at 0.

#### Attachments

• Screen Shot 2021-06-28 at 9.13.06 AM.png
18.5 KB · Views: 75
Gold Member
The Doppler shift from a moving target shows up in slow time as a phase shift from pulse to pulse, which you lose due to a conceptual issue--by starting with a power spectrum, you have no phase information at all. As a result, your estimate cannot reproduce the most important feature of the slow time data--the Doppler-induced phase shift.

tworitdash
tworitdash
The Doppler shift from a moving target shows up in slow time as a phase shift from pulse to pulse, which you lose due to a conceptual issue--by starting with a power spectrum, you have no phase information at all. As a result, your estimate cannot reproduce the most important feature of the slow time data--the Doppler-induced phase shift.
Good point. That I did only for the numerical one. Now, I multiplied a uniformly distributed random phase from ##0## to ##2\pi## with the square root of the power spectrum. I still get the same result because I am plotting the absolute values. To check if I have good phase information, I did the DFT again from these time-domain signals. I get the correct amplitude spectrum for both the analytical one and the numerical one.

My main problem is to find the time-domain equivalent of a Gaussian velocity spectrum. I proceed with this way but I am still not sure. That is why I checked for analytical solutions so that I can make a phase modulation later. I follow this paper to do the same. They start from a power spectrum. Then, they use a random phase that is uniformly distributed and multiply it with the square root of the power spectrum and then they do an IDFT to find the I and Q component of the time domain data.

https://journals.ametsoc.org/view/journals/apme/14/4/1520-0450_1975_014_0619_sowdsa_2_0_co_2.xml

The code with the random phase in the amplitude spectrum is given below. It is a modified version.

Modified:
clear;
% close all;

lambda = 0.03;

PRT = 1e-3;
mu = 5; % Mean Doppler
sigma = 0.2;
v_amb = 7.5;

N = 60000; % Total number of points in time axis

t1 = 0:PRT:(N - 1)*PRT; % Time axis
vel_axis = linspace(-v_amb, v_amb, N); % velocity axis for the entire rotation

% vel_axis_hs = linspace(-v_amb, v_amb, hs); % velocity axis for one beamwidth
% t1 = -N/2*PRT:PRT:(N/2 - 1)*PRT; % Time axis
% s_analyt_o = 1./sqrt(2*pi*(sigma).^2) .* exp(-(vel_axis - mu).^2./(2*(sigma).^2));

% SNR = 10^(30/10);
% X = rand(1, N);
% Theta = 2 .* pi * rand(1, N);

S_ = gaussmf(vel_axis, [sigma, mu]);

% Noise = sum(S_) ./ (N .* SNR);
% s_analyt_o = -(S_ + Noise) .* log(X);

s_analyt_o = S_;
s_num = ifft(ifftshift(sqrt(s_analyt_o) .* exp(1j .* 2 * pi * rand(1, N)))); % Numerical IFFT

s_analyt = ((2/pi)^(1/4) * sqrt(sigma * 4 * pi / lambda) .* exp(-t1 .* ((sigma * 4 * pi / lambda)^2 .* t1 - 1j .* mu .* 4 * pi / lambda))); % analytical IFFT

figure;
plot(t1, abs((s_analyt))/max(abs(s_analyt)), 'LineWidth', 2);

hold on; plot(t1, abs((s_num))/max(abs(s_num)), '*');
grid on;

xlim([-10 80]);

legend({'Wolfram Alpha Expression', 'MATLAB numerical'})

figure;

plot(vel_axis, abs(fftshift(fft(s_analyt)))/max(abs(fft(s_analyt)))); hold on;
plot(vel_axis, abs(fftshift(fft(s_num)))/max(abs(fft(s_num))));

legend({'Wolfram Alpha Expression', 'MATLAB numerical'})