How do I normalize the output of nufft?

4 views (last 30 days)
chris crowley
chris crowley on 26 Oct 2021
Edited: Chris Turnes on 4 Nov 2021
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 Comments
chris crowley
chris crowley on 29 Oct 2021
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 on 4 Nov 2021
Edited: Chris Turnes on 4 Nov 2021
I suspect that the answer I posted on the question linked below is essentially the answer to your question:
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).

Sign in to comment.

Answers (0)

Products


Release

R2021a

Community Treasure Hunt

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

Start Hunting!