nufft wrong result ?

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

Answer 1

[Paul]

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]

Answer 2

[Kraka Doros]

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.

[Paul]

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.

Answer 3

[Kraka Doros]

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

[Paul]

Example data from original question

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

Take the nufft of the sequence

f = (0:n-1)/n;
F = nufft(x,t,f);

Implenment the nufft using the defining sum for the nufft() doc page

F1 = zeros(n,1);

for kk = 1:numel(f)
    for jj = 1:n
        F1(kk) = F1(kk) +  x(jj)*exp(-1j*2*pi*t(jj)*f(kk));
    end
end

Compare the results

figure
subplot(211);
plot(real(F1-F),'-x');
subplot(212);
plot(imag(F1-F),'-o');

Not surprised that they are not the same because Matlab's actual algoriithm for sure won't be an implementation of that double loop. I suspect that at some point this simple approach will break down. But this at least gets you somewhere. In principal, it can be be implemented in the Symbolic Math Toolbox using vpa math as well. Hope it helps.

[Kraka Doros]

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 ?

Answer 4

[Kraka Doros]

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