CVApp Forum@View topic - frequency resolution of non-uniform FFT: NUFFT

frequency resolution of non-uniform FFT: NUFFT

Eigen,
Scilab,
Maxima(Mathematica-like)
R,Octave(matlib-like)
and all open source tools for mathematic and scientific fields.

frequency resolution of non-uniform FFT: NUFFT

Postby thomas » 2015年 Jun 24日, 18:23

Selection of convolution kernel in non-uniform fast Fourier transform for Fourier domain optical coherence tomography
https://www.osapublishing.org/vjbo/fulltext.cfm?uri=oe-19-27-26891&id=225906

https://www.facebook.com/groups/taiwan.open.platform.club/496540490495443/?notif_t=group_comment
http://www.math.nyu.edu/faculty/greengar/glee_nufft_sirev.pdf
http://math.hawaii.edu/home/theses/MA_2013_Brown.pdf
1. non-uniform fft,如果沒有zero padding, 表示NUFFT的整個觀察週期不變(沒有因為interpolation而改變觀察週期),
如果需要zero padding,觀察週期就變長, frequency resolution就變高,
但zero padding並不是真正提高frequency resolution.
詳情: http://www.bitweenie.com/listings/fft-zero-padding/

2. nufft本來有N個非等距取樣點,透過gaussain gridding找出2M個等距取樣點,
因為觀察週期不變, frequncy domain的 frequency resolution也不變.
而因為取樣點數變了(從N->2M), 所以最高工作頻率會改變,
這個時候要注意的反而是工作頻率會不會有alias.

3.透過測試程式跑模擬,一個簡單的consine function的FT, 確實是不管M變多變少(只改變sampling rate),
spectra peak都是出現在固定的整數上.
這表示改變M只改變最高的工作頻率,而frequency resolution沒有改變.

445.pdf
(628.75 KiB) Downloaded 372 times

446-447.pdf
(716.04 KiB) Downloaded 348 times
thomas
 
Posts: 534
Joined: 2013年 May 4日, 09:52

Output fk0. fk1 of non-uniform FFT(NUFFT)

Postby thomas » 2015年 Jun 25日, 10:59

type1 NUFFT
1. The image part of output fk0, fk1 will be "zero" if the xj[] is euqi-spaced.
2. The image part of output fk0, fk1 will NOT be "zero" if the xj[] is non euqi-spaced.

3. The simulation uses "cosine and sine" with the time position is disturbed by a random shift.

Code: Select all
float r_s=M_PI*2.0/(nj) * ((rand()%10 - 5.0 )/5.0f);
xj[j] = M_PI * k1/(nj/2) + r_s;//generate a random non-equispaced


The sine/cosine simulation:
Code: Select all
//http://forum.cvapp.org/viewtopic.php?f=21&t=431&sid=4e8c86381ce6b3a43a53d7d8f8b22058
void nufft_sine(void)
{
   int ier,iflag,j,k1,ms,nj;
   float xj[MX];
   float eps=1e-5;
   double err;
   COMPLEX8 cj[MX];
   COMPLEX8 fk0[MX],fk1[MX];
   double observingT=0.0f;
//
//     --------------------------------------------------
//     create some test data
//     --------------------------------------------------
   ms = 64;   //oversampling 2*ms in time domain
   nj = 77;
   //fortran array index from "1", but c is from "0"
   printf("input xj, nj and cj\n");
   printf("nj bin=%f\n", M_PI*2.0/(nj));
   for(k1 = -nj/2; k1 <= (nj-1)/2;k1++){
      j = k1+nj/2;
      //xj[j] = M_PI * cos(-M_PI*(j+1)/nj);
      //float r_s=M_PI*2.0/(nj*2) * ((rand()%10 - 5 )/5.0f);
      float r_s=M_PI*2.0/(nj) * ((rand()%10 - 5.0 )/5.0f);
      xj[j] = M_PI * k1/(nj/2) + r_s;//generate a random non-equispaced
      //cj[j] = cmplxf(cos(5.0*xj[j])+cos(1.0*xj[j]),
      //            sin(2.0*xj[j]));
      cj[j].r = cos(5.0*xj[j])+cos(1.0*xj[j]);
      cj[j].i = sin(2.0*xj[j]);
      
      printf("%d,%d:%f, %f,(%f,%f)\n", k1, j, r_s, xj[j], cj[j].r, cj[j].i);
   }
   printf("dirft1d1f_\n");
   iflag=-1;   //forward
   dirft1d1f_(&nj, xj, cj, &iflag, &ms, fk0);
   printf("nufft1d1ff90_ffte_\n");
   iflag=-1;   //forward
   nufft1d1ff90_ffte_(&nj, xj, cj, &iflag, &eps, &ms, fk1, &ier);
   errcompf(fk0, fk1, ms, &err);
   printf("ier=%d, err=%e\n\n", ier, err);
   printf("nj=%d, ms=%d\n",nj,ms);
   for(int i = 0; i < ms; i ++){
      printf("%d:fk0=(%f,%f),fk1=(%f,%f)\n", i, fk0[i].r,fk0[i].i,
         fk1[i].r,fk1[i].i);
      printf("(%f,%f)\n",cabsf_sq(fk0[i]), cabsf_sq(fk1[i]));
   }
}
thomas
 
Posts: 534
Joined: 2013年 May 4日, 09:52


Return to Mathematic, Engineering and Science tools

Who is online

Users browsing this forum: No registered users and 1 guest

cron