# nufft wrong result ?

조회 수: 5(최근 30일)
Kraka Doros 2022년 5월 8일
편집: Kraka Doros 2022년 5월 14일
% 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)

댓글을 달려면 로그인하십시오.

### 채택된 답변

Paul 2022년 5월 8일
편집: Paul 2022년 5월 8일
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]

댓글을 달려면 로그인하십시오.

### 추가 답변(3개)

Kraka Doros 2022년 5월 8일
편집: Kraka Doros 2022년 5월 8일
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표시숨기기 없음
Paul 2022년 5월 8일
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.

댓글을 달려면 로그인하십시오.

Kraka Doros 2022년 5월 9일
편집: Kraka Doros 2022년 5월 9일
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표시숨기기 이전 댓글 수: 1
Kraka Doros 2022년 5월 10일
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 ?

댓글을 달려면 로그인하십시오.

Kraka Doros 2022년 5월 14일
편집: Kraka Doros 2022년 5월 14일
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....

댓글을 달려면 로그인하십시오.

### Community Treasure Hunt

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

Start Hunting!

Translated by