# Huge scaling factor on amplitude in NUFFT vs FFT

5 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 1 commentaire plus ancienMasquer 1 commentaire plus ancien
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 4 commentaires plus anciensMasquer 4 commentaires plus anciens
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!