7 views (last 30 days)

Show older comments

Tommaso Pettinari on 14 May 2023

Commented: Tommaso Pettinari on 26 May 2023

- spectrum_fft_raw.fig
- spectrum_nufft_ot.fig
- spectrum_nufft_raw.fig
- holesAnalysis.m

Hi, I have a couple of questions about FFT and NUFFT.

- 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?
- 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

##### 2 Comments Show NoneHide None

Show NoneHide None

Paul on 14 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2745939

Open in MATLAB Online

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 on 15 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2746834

Edited: Tommaso Pettinari on 15 May 2023

- clean_frames.mat
- raw_trajectory.mat
- raw_trajectory_ot.mat

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!

Sign in to comment.

Sign in to answer this question.

### Answers (1)

Paul on 15 May 2023

Open in MATLAB Online

- clean_frames.mat
- raw_trajectory.mat
- raw_trajectory_ot.mat

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

##### 5 Comments Show 3 older commentsHide 3 older comments

Show 3 older commentsHide 3 older comments

Tommaso Pettinari on 15 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2747179

Edited: Tommaso Pettinari on 15 May 2023

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 on 15 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2747424

Edited: Paul on 17 May 2023

Open in MATLAB Online

- raw_trajectory.mat
- raw_trajectory_ot.mat
- clean_frames.mat

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 on 23 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2756689

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 on 24 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2757529

Edited: Paul on 24 May 2023

Open in MATLAB Online

- clean_frames.mat
- raw_trajectory.mat
- raw_trajectory_ot.mat

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 on 26 May 2023

#### Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/1963154-huge-scaling-factor-on-amplitude-in-nufft-vs-fft#comment_2761189

Thanks for the answers and your time, Paul.

Sign in to comment.

Sign in to answer this question.

### See Also

### Categories

MATLABMathematicsFourier Analysis and Filtering

Find more on **Fourier Analysis and Filtering** in Help Center and File Exchange

### Tags

- fft
- nufft

### Community Treasure Hunt

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

Start Hunting!

**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- Deutsch
- English
- Français

- United Kingdom(English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)

Contact your local office