# nufft wrong result ?

12 vues (au cours des 30 derniers jours)
Kraka Doros le 8 Mai 2022
Modifié(e) : Kraka Doros le 14 Mai 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)
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

Paul le 8 Mai 2022
Modifié(e) : Paul le 8 Mai 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 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Plus de réponses (3)

Kraka Doros le 8 Mai 2022
Modifié(e) : Kraka Doros le 8 Mai 2022
Hi Paul.
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 commentaireAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens
Paul le 8 Mai 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.

Connectez-vous pour commenter.

Kraka Doros le 9 Mai 2022
Modifié(e) : Kraka Doros le 9 Mai 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 commentairesAfficher AucuneMasquer Aucune
Paul le 9 Mai 2022
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 le 10 Mai 2022
Modifié(e) : Kraka Doros le 10 Mai 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)
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 ?

Connectez-vous pour commenter.

Kraka Doros le 14 Mai 2022
Modifié(e) : Kraka Doros le 14 Mai 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....
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Develop Apps Using App Designer dans Help Center et File Exchange

### Community Treasure Hunt

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

Start Hunting!

Translated by