nufft wrong result ?
16 views (last 30 days)
Show older comments
% my code :
% Create 256 not equally spaced time points between 1 and 1000 :
n=256;
t=sort(randperm(1000,n)');
% I received a column of 256 points with first and last points :
% t(1), t(n)
% I want the first point to be 0 :
t=t-t(1);
% Take values of x=sin(2*pi*0.025*t) on the 256 time points :
x=sin(2*pi*0.025*t);
% Pretend that these 256 values are my samples that i took in non equally
% spaced time points, without knowing the frequency and want to use
% nufft to calculate f of the sin (f=0.025) :
F=nufft(x,t,(0:n-1)/n);
% plot abs of F, over frequency bins with fresolution=1/T
% were T is the total time of sampling, t(256) :
T=t(256);
fresolution=1/T;
plot((0:n-1)/T,abs(F)/2*n)
% receiving wrong answer (f=0.006)
0 Comments
Accepted Answer
Paul
on 8 May 2022
Edited: Paul
on 8 May 2022
Hello Kraka,
Have you tried setting up your query points to exactly include f = 0.025? Also, it wasn't clear why the plot command had x-values (0:(n-1))/T?
% Create 256 not equally spaced time points between 1 and 1000 :
n=256;
rng(100);
t=sort(randperm(1000,n)');
% I received a column of 256 points with first and last points :
% t(1), t(n)
% I want the first point to be 0 :
t=t-t(1);
% Take values of x=sin(2*pi*0.025*t) on the 256 time points :
x=sin(2*pi*0.025*t);
% Pretend that these 256 values are my samples that i took in non equally
% spaced time points, without knowing the frequency and want to use
% nufft to calculate f of the sin (f=0.025) :
f = (0:n-1)/n;
F=nufft(x,t,f);
% plot abs of F, over frequency bins with fresolution=1/T
% were T is the total time of sampling, t(256) :
% original plot, changed to stem
T=t(256);
fresolution=1/T;
figure
stem((0:n-1)/T,abs(F)/2*n)
xlim([0 0.1])
% plot again, use the actual query frequencies, spike is close to .025
figure
stem((0:n-1)/n,abs(F)/2*n)
xlim([0 .1])
% exact query frequency at 0.025
f = 0:.005:.5;
F=nufft(x,t,f);
figure
stem(f,abs(F)/2*n)
xlim([0 .1]
0 Comments
More Answers (3)
Kraka Doros
on 8 May 2022
Edited: Kraka Doros
on 8 May 2022
1 Comment
Paul
on 8 May 2022
I'm afraid I can't help with this. Have you looked at the doc page for nufft()? It shows the equation used to compute the nufft that should be very straightforward to implement in Matlab. I'm sure that there are lots of issues to consider wrt to precision and efficiency, but maybe the simple approach would be suitable for your application.
See Also
Categories
Find more on Multirate Signal Processing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!