nufft wrong result ?

9 views (last 30 days)
Kraka Doros
Kraka Doros on 8 May 2022
Edited: Kraka Doros on 14 May 2022
% 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)

Accepted Answer

Paul
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]

More Answers (3)

Kraka Doros
Kraka Doros on 8 May 2022
Edited: Kraka Doros on 8 May 2022
Hi Paul.
Thanks for your detailed and well presented answer. You are exactly correct. After i read comment, retried and received the correct answer (as already you show in your example).
I am wondering if you can help me with the same above example (trying find the 0.025 frequency), but in Matlab2015. I used Matlab online for my question example. This will end in a few days. In my PC (has a 32-bit Windows OS), i am using Matlab2015, which do not include nufft. Instead i compile a nufft algorythm by Matthew Ferrara. This algorythm has these inputs :
f: frequency-domain data (a complex Mx1 vector) unwrapped from a matrix knots: k-space locations at which the data were measured (an Mx2 vector). This data must be in double format.
knots: the frequency locations of the data points. These locations should be normalized to correspond to the grid boundaries [-N/2, N/2 -1/N], N even. If the knots are not scaled properly, they will be shifted and scaled into this normalized form.
N = [Nx]: the 1x1 vector denoting the size of the spatial grid in the image domain. This parameter will determine the spatial extent of the image.
Optional Inputs (highly recommended):
accuracy: a positive integer indicating the desired number of digits of accuracy
GridListx:(column vector) The x-locations of the frequency grid (not yet scaled by c or 2*pi) onto which the data should be interpolated.
I do not know how to match-correspond my data with these inputs. So far i tried : t for f input, x for knots input and 256 for N input. Also 4 for accuracy input. But received wrong results.
Any help will be appreciate.
  1 Comment
Paul
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.

Sign in to comment.


Kraka Doros
Kraka Doros on 9 May 2022
Edited: Kraka Doros on 9 May 2022
I can see that NUFFT is given from the equation :
Unfortunately, i cannot even use the simple approach. I am trying to guess, what are represented by Xn , pn and fk. I completely, have no idea. And I do not know how to use this equation, for to take 256 complex numbers.
In description says that p are sample points from 0 to 1 and f are frequencies from 0 to N. I tried t (time points, from my previous example) for p, and x (values of sin(2*pi*0.025*t)) for x. For f, i tried the bins of query frequencies. Confused...
  2 Comments
Kraka Doros
Kraka Doros on 10 May 2022
Edited: Kraka Doros on 10 May 2022
This is very impressive for me. It was a lot of work to understood the codes (i am new in Matlab), but it worth. I like it !
So, after the computation of F1, i made :
plot((0:n-1)/n,abs(F1)/2*n)
were previously the f, has defined as f=(0:n-1)/n , and received this :
with peak very close to 0.025 :
But if i change f to f=(0:n-1)/T , were T=t(n) , and plot to :
plot((0:n-1)/T,abs(F1)/2*n)
i received this :
which, i thing, is more clear and its peak is a little bit closer to 0.025 :
Anyway, thank you very, very much for your spending time, to help me understand.
If you please i want a clarification on something :
1) numel(f) and n, are not the same ? Then why you used :
kk=1:numel(f) , jj=1:n and not
kk=1:numel(f) , jj=1:numel(f) or
kk=1:n , jj=1:n ?
2) Why did you add F1(kk) in the left of the right part of
F1(kk) = x(jj)*exp(-1j*2*pi*t(jj)*f(kk)); to make it like :
(F1(kk) = F1(kk) + x(jj)*exp(-1j*2*pi*t(jj)*f(kk));) ?
3) Why, for x and t, use (jj), and for f use (kk) ? Cannot use (jj) or (kk) for all ?

Sign in to comment.


Kraka Doros
Kraka Doros on 14 May 2022
Edited: Kraka Doros on 14 May 2022
OK. No need for clarification. It was very complicated, but i got it . I understood why you add F1(kk) in the left of the right part of the equation and why you used both kk and jj.
The only thing i dont understand is the 1) : Is there a case that numel(f) not equal to n ?
I suspect that you used it like so, because you want to use a general form, when somebody do not want to take reasults for n "f bins", but for something else, for example, n/2....

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by