Huge scaling factor on amplitude in NUFFT vs FFT

10 vues (au cours des 30 derniers jours)
Tommaso Pettinari le 14 Mai 2023
Hi, I have a couple of questions about FFT and NUFFT.
1. I tried to compute a spectral analysis of my data (vectors of more than 30000 samples, due to a sampling frequency of 240 Hz), using both FFT and NUFFT algorithm, on the same vector of sampled data. The results had the same shape and were identical, except for a high scaling factor (10^4). Moreover, using only NUFFT on a non-uniform vector with only data exceeding a certain threshold (0.95, in my case) from the same previous array, the scaling factor is obviously different, but still present. Is there a quick way to find and possibly correct this scaling factor?
2. On the NUFFT's spectral analysis, there's a higher number of harmonics at low frequencies. Could this be due to the fact that, with the threshold stripping, I have cleaned the signal, so I can see better the low-frequency components of spectum?
2 commentairesAfficher AucuneMasquer Aucune
Paul le 14 Mai 2023
Hi Tommaso,
I looked at both m-files and didn't see any calls to fft or nufft. Maybe I missed it.
For the first part of item 1, can you provide the data vector in a .mat file and show the fft and nufft and plotting commands you're using on that data? Preferably directly in the Answers editor rather than attaching an m-file. Copy/paste if easier.
fft and nufft should return the same result with uniformly spaced data
rng(100);
x = rand(1,10);
X1 = fft(x);
X2 = nufft(x);
norm(X2-X1)
ans = 0
% assuming 240 Hz and explicitly defining the t and f inputs to nufft
Ts = 1/240;
t = (0:9)*Ts;
Fs = 240;
f = (0:9)/10*Fs;
X2 = nufft(x,t,f);
norm(X2-X1)
ans = 0
Tommaso Pettinari le 15 Mai 2023
Modifié(e) : Tommaso Pettinari le 15 Mai 2023
Hi Paul.
I am sending you a simplified version of the code in the space below and three variables to open in the workspace before running the code.
[m, n]=size(raw_trajectory);
fs = 240;
%% Spectral analisys with nufft
frequencyPSD = zeros(n/2, 3);
for i=1:1
% Signal's amplitude spectrum
dft_trajectory = abs(nufft(raw_trajectory_ot(~isnan(raw_trajectory_ot(:,i)),i),clean_frames(i,~isnan(raw_trajectory_ot(:,i)))));
spectrum = dft_trajectory(1:round(m/2)+1);
spectrum(2:end-1) = 2*spectrum(2:end-1);
f = 0:fs/m:fs/2;
figure
subplot(2,1,1)
plot(f(2:end), spectrum(2:end));
title('Amplitude spectrum with nufft'); xlabel('Frequency (Hz)'); ylabel('Modulo del Segnale'); axis tight; grid on;
% Signal's power spectrum
psd = dft_trajectory(1:round(m/2)+1); % amplitude spectrum
psd = 1/(fs*m) * psd.^2; % abs of amplitude spectrum
psd(2:end-1) = 2 * psd(2:end-1); % multiplication by 2 of positive power spectrum part
subplot(2,1,2)
plot(f(2:end), psd(2:end)); % plot of power spectrum
title('Power spectrum with nufft'); xlabel('Frequency (Hz)'); ylabel('Power (W)'); axis tight; grid on;
hold off
end
%% Analisi spettrale con fft
frequencyPSD = zeros(n/2, 3);
figure
for i=1:1
% Amplitude spectrum
dft_trajectory = abs(fft(raw_trajectory(:,i))/m);
spectrum(:,i) = dft_trajectory(1:round(m/2)+1);
spectrum(2:end-1,i) = 2*spectrum(2:end-1,i);
f = 0:fs/m:fs/2;
subplot(211)
plot(f(2:end), spectrum(2:end,i));
title('Ampitude spectrum with fft'); xlabel('Frequency (Hz)'); ylabel('Modulo del Segnale'); axis tight; grid on;
end
for i=1:1
% Power spectrum
dft_trajectory = abs(fft(raw_trajectory(:,i))/m);
psd(:,i) = dft_trajectory(1:round(m/2)+1);
psd(:,i) = 1/(fs*m) * psd(:,i).^2;
psd(2:end-1,i) = 2 * psd(2:end-1,i);
subplot(212)
plot(f(2:end), psd(2:end,i));
title('Power spectrum with fft'); xlabel('Frequency (Hz)'); ylabel('Power (W)'); axis tight; grid on;
end
%% Spectral analysis with nufft on raw trajectories
frequencyPSD = zeros(n/2, 3);
for i=1:1
% Amplitude spectrum
dft_trajectory = abs(nufft(raw_trajectory(:,i)));
spectrum = dft_trajectory(1:round(m/2)+1);
spectrum(2:end-1) = 2*spectrum(2:end-1);
f = 0:fs/m:fs/2;
figure
subplot(2,1,1)
plot(f(2:end), spectrum(2:end));
title('Amplitude spectrum with nufft on raw trajectory'); xlabel('Frequency (Hz)'); ylabel('Modulo del Segnale'); axis tight; grid on;
% Power spectrum
psd = dft_trajectory(1:round(m/2)+1);
psd = 1/(fs*m) * psd.^2;
psd(2:end-1) = 2 * psd(2:end-1);
subplot(2,1,2)
plot(f(2:end), psd(2:end));
title('Power spectrum with nufft on raw trajectory'); xlabel('Frequency (Hz)'); ylabel('Power (W)'); axis tight; grid on;
hold off
end
FFT and NUFFT are called at lines 7, 31, 42 and 56.

Connectez-vous pour commenter.

Réponses (1)

Paul le 15 Mai 2023
In the fft section, the variable dft_trajectory is defined as
% dft_trajectory = abs(fft(raw_trajectory(:,i))/m);
but the nufft section does not include the division by m in its calculation of dft_trajectory. That alone accounts for a big difference.
Also, the raw_trajectory data has some outliers and data drops that should be cleaned up
[m, n]=size(raw_trajectory);
fs = 240;
figure
tvec = (0:m-1)/fs;
% I was playing aroung here to deal with the outliers, but the data drop at
% around t = 1 second last for quite a long duration.
% raw_trajectory(:,1) = filloutliers(raw_trajectory(:,1),'linear','movmedian',100);
plot(subplot(211),clean_frames(1,:),raw_trajectory_ot(:,1))
plot(subplot(212),tvec,raw_trajectory(:,1))
5 commentairesAfficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens
Paul le 24 Mai 2023
Modifié(e) : Paul le 24 Mai 2023
fft and nufft do not use different scaling factors. However, nufft can be very sensitive to a non-zero mean when there are lots of gaps in the data, as is the case in this problem.
[m, n] = size(raw_trajectory);
fs = 240;
DFT of original data
i = 1;
dft_fft = fft(raw_trajectory(:,i));
f_fft = (0:numel(dft_fft)-1)/numel(dft_fft)*fs;
NUFFT of the "cleaned data"
ytemp = raw_trajectory_ot(~isnan(raw_trajectory_ot(:,i)),i);
ttemp = -1/240+clean_frames(i,~isnan(raw_trajectory_ot(:,i)));
f_nufft = f_fft;
dft_nufft = nufft(ytemp,ttemp,f_nufft);
Compare their magnitude responses
figure
plot(f_nufft,abs(dft_nufft),f_fft,abs(dft_fft)),grid
xlim([5 235]),ylim([0 1e5])
legend('nufft','fft')
The question is, "why does the nufft of the cleaned data have larger magnitude than the fft of the raw data?"
First, let's show that the nufft is essentially the same as the fft when the bad records that are nan-filled are instead zero-filled.
Fill the bad records with zeros and compute the DFT
ytemp1 = raw_trajectory_ot(:,1);
ytemp1(isnan(ytemp1)) = 0;
ttemp1 = (0:m-1)/fs;
dft_fft1 = fft(ytemp1);
Compare the relative error between the nufft of the nan-filled signal and the fft of the zero-filled signal. I'm pretty sure that that with full precision the difference would be identically zero, but with rounding errors the relative error is always less than 4e-9, and usually much less than that.
figure
plot(abs((dft_nufft-dft_fft1)./dft_fft1))
So, the equivalent question would be, "why does the fft of the zero-filled signal" have such large amplitude?"
Plot the zero-filled signal
figure
plot(ttemp1,ytemp1)
The zero-filling causes lots of transients, because the mean of the clean signal is around 1200 and there are lots of bad records that are being filled.
To get rid of the transients, let's subtract the mean of the clean data, then zero-fill the bad records, then add back the mean, then take the DFT
ytemp2 = raw_trajectory_ot(:,1);
mytemp2 = mean(ytemp2,'omitnan');
ytemp2 = ytemp2 - mytemp2;
ytemp2(isnan(ytemp2)) = 0;
ytemp2 = ytemp2 + mytemp2;
ttemp2 = (0:m-1)/fs;
dft_fft2 = fft(ytemp2);
figure
plot(f_nufft,abs(dft_nufft),f_fft,abs(dft_fft),f_fft,abs(dft_fft2)),grid
xlim([5 235]),ylim([0 1e5])
legend('nufft','fft','fft2')
We see that subtracting the mean, followed by zero-filling, followed by adding back the mean results in a much lower amplitude spectrum. This suggests that we should subtract the mean before taking the nufft if we don't want to do zero-filling
dft_nufft1 = nufft(ytemp-mean(ytemp),ttemp,f_nufft);
figure
plot(f_nufft,abs(dft_nufft),f_fft,abs(dft_fft),f_nufft,abs(dft_nufft1))
ylim([0 1e5])
legend('nufft','fft','nufft1')
Just compare the fft of the original signal and the nufft of the mean-subtracted signal
figure
plot(f_nufft,nan+abs(dft_nufft),f_fft,abs(dft_fft),f_nufft,abs(dft_nufft1))
xlim([5 235])
legend('nufft','fft','nufft1')
Tommaso Pettinari le 26 Mai 2023

Connectez-vous pour commenter.

Catégories

En savoir plus sur Fourier Analysis and Filtering dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by