How do I normalize the output of nufft?

9 次查看（过去 30 天）
chris crowley
chris crowley 2021-10-26
编辑： Chris Turnes 2021-11-4
I am trying to figure out what normilization of the output of nufft(x,t) satisfied Parseval's theorem and am unable to figure it out. I tested uniformly sampled data and can see that nufft clearly is not normalized the same way as fft, but am unable to figure out correct normalization.
To illustrate my confusion let me first define a signal that I can calculate the energy of analytically.
%% define a signal
dt = 1/1000; % [Hz]
T = 3; % [s]
t = 0:dt:T; % [s]
f1 = 123; % [Hz]
A1 = 1/sqrt(2); % [-]
f2 = 40; % [Hz]
A2 = 1/sqrt(2); % [-]
sig = sqrt(2)*(A1.*cos(t*2*pi*f1) +A2.*cos(t*2*pi*f2));
L = length(sig);
% True energy in the time domain
E_true = A1.^2 + A2.^2;
Now, lets compute the non-uniform FFT of this uniform data useing nufft
%% Do the FFT
% Take the non-uniform FFT
NUFT = nufft(sig,t);
% Take the standard FFT
UFT = fft(sig);
Once I have computer the approximate Fourier Transform, let's calculate the power spectrum and normalize it the way FFT would require.
% calculate power spectrum
P_NU = abs(NUFT(round(1:L/2+1))).^2;
P_U = abs(UFT(round(1:L/2+1))).^2;
% create the freqs in freq domain
f = (0:L/2)/(L*dt);
% normalize both as if it were FFT
P_NU = P_NU*2/(L^2);
P_U = P_U*2/(L^2);
Now, lets look at Parsevals theroem
%% check Parseval's theorem
% Parseval's theorem is just energy conservation
% Energy in the time domain = energy in the freq domain
% The analytical Energy
disp(['The true energy = ',num2str(E_true)])
% Energy in the time signal
E_t = trapz(t,sig.^2)/(dt*(L-1));
disp(['Energy in time domain = ',num2str(E_t)])
% Energy in the freq domain
E_f_U = trapz(f,P_U)*(dt*(L-1));
disp(['Energy in freq domain (fft)= ',num2str(E_f_U)])
% Energy in the freq domain
E_f_NU = trapz(f,P_NU)*(dt*(L-1));
disp(['Energy in freq domain (nufft)= ',num2str(E_f_NU)])
disp(' ')
Running this code yeilds:
The true energy = 1
Energy in time domain = 1
Energy in freq domain (fft)= 1.0007
Energy in freq domain (nufft)= 0.0006591
How am I supposed to normalize the output of nufft?
  2 个评论
chris crowley
chris crowley 2021-10-29
It seems to work great if dt = 1, but if dt<1 ploting the abs(nufft(sig,t)) looks very different from abs(fft(sig)). By maping the time variable to a new time variable with integer dt seems to fix it.
NUFT = nufft(sig,t/min(diff(t)));
Chris Turnes
Chris Turnes 2021-11-4
编辑：Chris Turnes 2021-11-4
I suspect that the answer I posted on the question linked below is essentially the answer to your question:
https://www.mathworks.com/matlabcentral/answers/1453969-can-we-have-a-very-simple-and-clear-comparison-example-on-the-functions-fft-and-nufft-works?s_tid=srchtitle
The short version is that the default arguments for the sample/frequency points of nufft are 0:(N-1) and (0:(N-1))/N, respectively. If your time vector has a different scale than "sample units", then you have to manually scale your time units (as you've discovered).

请先登录，再进行评论。

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Fourier Analysis and Filtering 的更多信息

标签

产品


版本

R2021a

Community Treasure Hunt

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

Start Hunting!

Translated by