MATLAB IFFT doesn't match the analytical one

  • MATLAB
  • Thread starter tworitdash
  • Start date
  • #1
68
14
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)));


Numerical_IFFT.png

Analytical_IFFT.png
 

Answers and Replies

  • #2
jasonRF
Science Advisor
Gold Member
1,462
519
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
 
  • Like
Likes hutchphd and Twigg
  • #3
13,030
6,916
Also it would be nice if the diagrams were labeled so we know if they came from MATLAB or Mathematica.
 
  • #4
68
14
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?
 
  • #5
Delta2
Homework Helper
Insights Author
Gold Member
4,532
1,835
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.
 
  • #6
marcusl
Science Advisor
Gold Member
2,774
432
Shouldn’t you ifftshift after transforming, not before?
 
  • #7
68
14
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,

[tex] S(v) = \frac{1}{\sqrt(2\pi \sigma^2)} e^{-\frac{(v - \mu)^2}{2\sigma^2}} [/tex]

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,

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

The N is the number of pulses, the [itex] \delta t [/itex] is the pulse repetition time, [itex] \lambda [/itex] 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.

Screen Shot 2021-06-28 at 9.20.37 AM.png





The expression I used is this,

[tex] 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} )}[/tex]


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.


Screen Shot 2021-06-28 at 9.43.26 AM.png



Screen Shot 2021-06-28 at 9.43.56 AM.png



Screen Shot 2021-06-28 at 9.44.19 AM.png
 

Attachments

  • Screen Shot 2021-06-28 at 9.13.06 AM.png
    Screen Shot 2021-06-28 at 9.13.06 AM.png
    18.5 KB · Views: 26
  • #8
marcusl
Science Advisor
Gold Member
2,774
432
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.
 
  • #9
68
14
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'})
 
  • #10
marcusl
Science Advisor
Gold Member
2,774
432
Dusan Zrnic is a giant in the weather radar field (he literally wrote "the book" on Doppler weather radar), so you are looking at a good source.
 
  • #11
68
14
Dusan Zrnic is a giant in the weather radar field (he literally wrote "the book" on Doppler weather radar), so you are looking at a good source.
Yes, it is a very good source indeed. My actual model is based on his idea from that paper. The problem I am facing is in the next step. This model of Zrnic is valid when the radar beam stays still and points to one direction in space. I am trying to create a model which considers the rotation of the beam with a beam width of roughly two degrees in azimuth. I can deduce the time domain signal for a monochromatic (Dirac) weather model in case of a rotating radar, but I am unable to do it for a Gaussian spectrum. I don't know how to modulate the time domain signal from this paper of Zrnic with a Gaussian Power spectrum to accommodate the rotation of the radar. Do you have any idea where can I find good sources for this kind of problem? This is the first step in my Ph.D. and I am unable to find a good reference for that. All references I found are more or less related to a rotating target with constant velocity, but no one talks about a distribution.
 

Related Threads on MATLAB IFFT doesn't match the analytical one

  • Last Post
Replies
3
Views
2K
Replies
4
Views
10K
Replies
0
Views
7K
Replies
2
Views
6K
Replies
3
Views
4K
Replies
2
Views
9K
Replies
1
Views
2K
Replies
7
Views
2K
Top