Huge scaling factor on amplitude in NUFFT vs FFT

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?

Thanks in advance

[Paul]

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]

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.

Thanks for your help!

Answer 1

[Paul]

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

load raw_trajectory
load clean_frames
load raw_trajectory_ot

[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))

[Tommaso Pettinari]

Yes, dividing by m the problem about the scaling factor between NUFFT and FFT results on "raw_trajectory" is solved, but the real problem is the different scaling using NUFFT on "raw_trajectory_ot" as data and "clean_frames" as samples (the cleaned version of "raw_trajectory", NaN values are excluded in NUFFT calculation); there is a scaling factor of about 10^2. Is this due to the holes (the frames I excluded)?

NUFFT on cleaned data:

dft_trajectory = abs(nufft(raw_trajectory_ot(~isnan(raw_trajectory_ot(:,i)),i)/m,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));

FFT on raw data:

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));

Is there any explanation on the second question?

Thanks for your help!

[Paul]

Ok. Now I think I understand. raw_trajectory_ot has NaNs filled in for the bad frames, but is otherwise the the same as raw_trajectory. Just looking by eye, it seemed that some frames in raw_trajectory_ot were nan-filled when the raw_trajectory data looked pretty good.

Having said that, the best I can offer, for now, is the following.

load raw_trajectory.mat
load raw_trajectory_ot.mat
load clean_frames.mat

[m, n]=size(raw_trajectory);
fs = 240;

For some reason, clean_frames starts at 1/240, so I subtracted 1/240 below to start it at zero. I don't know if that was correct. If not, then the fft calculation would have to be adjusted to account for the fact that the first sample doesn't start at 0. That would only impact the phase.

clean_frames(1,1:5)

ans = 1×5

0.0042 0.0083 0.0125 0.0167 0.0208

DFT of raw_trajectory and its associated frequency vector

i = 1;
dft_fft   = abs(fft(raw_trajectory(:,i)));
f_fft = (0:numel(dft_fft)-1)/numel(dft_fft)*fs;

NUDFT of raw_trajectory_ot. Here, it seems to me that if the sample points are in units of seconds, the query points have to be in units of Hz. So I'll use the f_fft vector as the query points for nufft.

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 = abs(nufft(ytemp,ttemp,f_nufft));


figure
plot(subplot(211),f_nufft,abs(dft_nufft)),grid,ylim([0 1e7])
plot(subplot(212),f_fft,abs(dft_fft)),grid,xlim([1 240])

figure
plot(f_nufft,abs(dft_nufft),f_fft,abs(dft_fft)),grid,xlim([5 235]),ylim([0 1e5])
legend('nufft','fft')

The nufft plot does not reveal any frequency structure in the data. Is that expected? I don't know why the nufft data is so much larger than the fft data.

Another way to look at this is to use spectrogram to get a frequency/time plot.

First, plot the data

tvec = (0:m-1)/fs;
figure
plot(tvec/60,raw_trajectory(:,1))

Now, the spectrogram

figure
nspec = 1024;
spectrogram(raw_trajectory(:,1),nspec,nspec/2,nspec,fs)

Other than the horizontal, all-frequency lines at the times where there are massive transients in the data, it doesn't look theres any particular frequencies that stand out. Maybe playing around with the spectrogram inputs can tease something out of the data.

Is there any particular frequency content that is expected?

[Tommaso Pettinari]

Hi, thanks for the explanation. For what concern the frequency structure of the data, it isn't expected any particular structure, or better it's to search if there is .

As regards a particular expected frequence content, these data are from human baby movements, so the signal sholud be below the 5-10 Hz maximum; over, it should be noise.

In any case, the strange thing is the big scaling factor between fft and nufft, that you highlighted in the subplot as well, of which we haven't found any explanation up to now.

[Paul]

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.

load raw_trajectory.mat
load raw_trajectory_ot.mat
load clean_frames.mat
[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]

Thanks for the answers and your time, Paul.