API tools faq
paste
deltaluca

Untitled

deltaluca
May 1st, 2013
76
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C 7.63 KB | None | 0 0
raw download clone embed print report
  1. #define Xr(i) x[dX*(i)]
  2. #define Xi(i) x[dX*(i)+1]
  3. #define Wr(i) w[dX*(i)]
  4. #define Wi(i) w[dX*(i)+1]
  5. #define Yr(i) y[dY*(i)]
  6. #define Yi(i) y[dY*(i)+1]
  7. #define Wpr(i) Wp[dWp*(i)]
  8. #define Wpi(i) Wp[dWp*(i)+1]
  9. #define Wnr(i) Wn[dWp*(i)]
  10. #define Wni(i) Wn[dWp*(i)+1]
  11. #define Ynr(i) yn[dY*(i)]
  12. #define Yni(i) yn[dY*(i)+1]
  13.  
  14. /* Fast Fourier Series */
  15. void FFT(int N, double *x, double *w, int dX, double *y, int dY, double *Wp, int dWp) {
  16.     int i, Nyd, NWpd;
  17.     double *xe, *ye, *yn, *We, *Wn, re, im, RE, IM, ra, ia;
  18.     double t0, t1, t2, t3;
  19.  
  20.     switch (N) {
  21.     case 1:
  22.         /* CN = [ 1 ] */
  23.         Xr(0) = Yr(0);
  24.         Xi(0) = Yi(0);
  25.         break;
  26. #if silly_fast
  27.     case 2:
  28.         /* CN = [ 1  1 ]
  29.                 [ 1 -1 ] */
  30.         Xr(0) = Yr(0) + Yr(1);
  31.         Xi(0) = Yi(0) + Yi(1);
  32.         Xr(1) = Yr(0) - Yr(1);
  33.         Xi(1) = Yi(0) - Yi(1);
  34.         break;
  35. #endif
  36.  
  37. #if faster
  38.     case 3:
  39.         /* CN = [ 1  1  1 ]
  40.                 [ 1  w  w*]
  41.                 [ 1  w* w ]
  42.         */
  43.         re = Wpr(1);
  44.         im = Wpi(1);
  45.         Xr(0) = Yr(0) + Yr(1) + Yr(2);
  46.         Xi(0) = Yi(0) + Yi(1) + Yi(2);
  47.  
  48.         Xr(1) = Xr(2) = Yr(0) + (Yr(1) + Yr(2))*re;
  49.         Xi(1) = Xi(2) = Yi(0) + (Yi(1) + Yi(2))*re;
  50.         ra = (Yi(2) - Yi(1))*im;
  51.         ia = (Yr(1) - Yr(2))*im;
  52.         Xr(1) += ra;
  53.         Xi(1) += ia;
  54.         Xr(2) -= ra;
  55.         Xi(2) -= ia;
  56.         break;
  57. #endif
  58.  
  59. #if even_faster
  60.     case 4:
  61.         /* CN = [ 1  1  1  1 ]
  62.                 [ 1  i -1 -i ]
  63.                 [ 1 -1  1 -1 ]
  64.                 [ 1 -i -1  i ] */
  65.         Xr(0) = Xr(2) = Yr(0) + Yr(2);
  66.         Xi(0) = Xi(2) = Yi(0) + Yi(2);
  67.         ra = Yr(1) + Yr(3);
  68.         ia = Yi(1) + Yi(3);
  69.         Xr(0) += ra;
  70.         Xi(0) += ia;
  71.         Xr(2) -= ra;
  72.         Xi(2) -= ia;
  73.  
  74.         Xr(1) = Xr(3) = Yr(0) - Yr(2);
  75.         Xi(1) = Xi(3) = Yi(0) - Yi(2);
  76.         ra = Yi(3) - Yi(1);
  77.         ia = Yr(1) - Yr(3);
  78.         if (dWp > 0) { /* FastDFS */
  79.             Xr(1) += ra;
  80.             Xi(1) += ia;
  81.             Xr(3) -= ra;
  82.             Xi(3) -= ia;
  83.         }
  84.         else { /* FastDFT */
  85.             Xr(1) -= ra;
  86.             Xi(1) -= ia;
  87.             Xr(3) += ra;
  88.             Xi(3) += ia;
  89.         }
  90.         break;
  91. #endif
  92.  
  93. #if faster
  94.     case 5:
  95.         /* CN = [ 1  1  1  1  1 ]
  96.                 [ 1  w  q  q* w*]
  97.                 [ 1  q  w* w  q*]
  98.                 [ 1  q* w  w* q ]
  99.                 [ 1  w* q* q  w ]
  100.         q = w^2 */
  101.         re = Wpr(1);
  102.         im = Wpi(1);
  103.         RE = Wpr(2);
  104.         IM = Wpi(2);
  105.         Xr(0) = Yr(0) + Yr(1) + Yr(2) + Yr(3) + Yr(4);
  106.         Xi(0) = Yi(0) + Yi(1) + Yi(2) + Yi(3) + Yi(4);
  107.  
  108.         Xr(1) = Xr(4) = Yr(0) + (t0 = Yr(1) + Yr(4))*re + (t1 = Yr(2) + Yr(3))*RE;
  109.         Xi(1) = Xi(4) = Yi(0) + (t2 = Yi(1) + Yi(4))*re + (t3 = Yi(2) + Yi(3))*RE;
  110.         Xr(2) = Xr(3) = Yr(0) + t1*re + t0*RE;
  111.         Xi(2) = Xi(3) = Yi(0) + t3*re + t2*RE;
  112.  
  113.         ra = (t0 = Yi(4) - Yi(1))*im + (t1 = Yi(3) - Yi(2))*IM;
  114.         ia = (t2 = Yr(1) - Yr(4))*im + (t3 = Yr(2) - Yr(3))*IM;
  115.         Xr(1) += ra;
  116.         Xi(1) += ia;
  117.         Xr(4) -= ra;
  118.         Xi(4) -= ia;
  119.         ra = t0*IM - t1*im;
  120.         ia = t2*IM - t3*im;
  121.         Xr(2) += ra;
  122.         Xi(2) += ia;
  123.         Xr(3) -= ra;
  124.         Xi(3) -= ia;
  125.         break;
  126. #endif
  127.  
  128. #if silly_fast
  129.     case 6:
  130.         /* CN = [ 1  1  1  1  1  1 ]
  131.                 [ 1  w #w -1 -w  w*]
  132.                 [ 1 #w -w  1 #w -w ]
  133.                 [ 1 -1  1 -1  1 -1 ]
  134.                 [ 1 -w #w  1 -w #w ]
  135.                 [ 1 w* -w -1 #w  w ]
  136.         #w = -(w*) */
  137.         re = Wpr(1);
  138.         im = Wpi(1);
  139.  
  140.         Xr(0) = Xr(3) = Yr(0) + Yr(2) + Yr(4);
  141.         Xi(0) = Xi(3) = Yi(0) + Yi(2) + Yi(4);
  142.         ra = Yr(1) + Yr(3) + Yr(5);
  143.         ia = Yi(1) + Yi(3) + Yi(5);
  144.         Xr(0) += ra;
  145.         Xi(0) += ia;
  146.         Xr(3) -= ra;
  147.         Xi(3) -= ia;
  148.  
  149.         Xr(1) = Xr(2) = Xr(4) = Xr(5) = Yr(0) - re*(Yr(2) + Yr(4));
  150.         Xi(1) = Xi(2) = Xi(4) = Xi(5) = Yi(0) - re*(Yi(2) + Yi(4));
  151.         ra = im*(Yi(4) - Yi(2));
  152.         ia = im*(Yr(2) - Yr(4));
  153.         Xr(1) += ra; Xr(4) += ra;
  154.         Xi(1) += ia; Xi(4) += ia;
  155.         Xr(2) -= ra; Xr(5) -= ra;
  156.         Xi(2) -= ia; Xi(5) -= ia;
  157.  
  158.         ra = re*(Yr(1) + Yr(5)) - Yr(3);
  159.         ia = re*(Yi(1) + Yi(5)) - Yi(3);
  160.         Xr(1) += ra; Xr(4) -= ra;
  161.         Xi(1) += ia; Xi(4) -= ia;
  162.         Xr(2) -= ra; Xr(5) += ra;
  163.         Xi(2) -= ia; Xi(5) += ia;
  164.  
  165.         ra = im*(Yi(5) - Yi(1));
  166.         ia = im*(Yr(1) - Yr(5));
  167.         Xr(1) += ra; Xr(4) -= ra;
  168.         Xi(1) += ia; Xi(4) -= ia;
  169.         Xr(2) += ra; Xr(5) -= ra;
  170.         Xi(2) += ia; Xi(5) -= ia;
  171.         break;
  172. #endif
  173.     default:
  174.         if ((N%2) != 0) {
  175.             Nyd = N*dY;
  176.             NWpd = N*dWp;
  177.             ye = y + Nyd, We = Wp + NWpd;
  178.             Xr(0) = Yr(0);
  179.             Xi(0) = Yi(0);
  180.             for (yn = y+dY; yn < ye; yn += dY) {
  181.                 Xr(0) += Ynr(0);
  182.                 Xi(0) += Yni(0);
  183.             }
  184.             x += dX;
  185.             for (Wn = Wp, i = dWp; i*dWp < NWpd*dWp; i += dWp, x += dX, Wn = Wp) {
  186.                 Xr(0) = Yr(0);
  187.                 Xi(0) = Yi(0);
  188.                 Wn += i;
  189.                 for (yn = y+dY; yn < ye; yn += dY) {
  190.                     Xr(0) += Wnr(0)*Ynr(0) - Wni(0)*Yni(0);
  191.                     Xi(0) += Wnr(0)*Yni(0) + Wni(0)*Ynr(0);
  192.                     Wn += i;
  193.                     if (((size_t)Wn)*dWp >= ((size_t)We)*dWp)
  194.                         Wn -= NWpd;
  195.                 }
  196.             }
  197.         }
  198.         else {
  199.             FFT(N/2, w,       x,       dX, y,   dY*2, Wp,dWp*2);
  200.             FFT(N/2, w+N*dX/2,x+N*dX/2,dX, y+dY,dY*2, Wp,dWp*2);
  201. #if even_faster
  202.             xe = x + N*dX/2;
  203.             Xr(0)   = Wr(0) + Wr(N/2);
  204.             Xi(0)   = Wi(0) + Wi(N/2);
  205.             Xr(N/2) = Wr(0) - Wr(N/2);
  206.             Xi(N/2) = Wi(0) - Wi(N/2);
  207.             x += dX;
  208.             w += dX;
  209.             Wp += dWp;
  210.             for (; x < xe; x += dX, w += dX, Wp += dWp) {
  211.                 re = Wpr(0)*Wr(N/2) - Wpi(0)*Wi(N/2);
  212.                 im = Wpr(0)*Wi(N/2) + Wpi(0)*Wr(N/2);
  213.                 Xr(0)   = Wr(0) + re;
  214.                 Xi(0)   = Wi(0) + im;
  215.                 Xr(N/2) = Wr(0) - re;
  216.                 Xi(N/2) = Wi(0) - im;
  217.             }
  218. #else
  219.             for (Wn = Wp, xe = x + N*dX/2; x < xe; x += dX, w += dX, Wn += dWp) {
  220.                 /* FastDFT calls FFT with Wp pointing to 'just past' the end of the
  221.                    true Wp array, and accessing Wp[N] for positive N is not valid.
  222.                    At this point, we've not optimised FFT algorithm to 'know' that
  223.                    Wp[0] is actually 1, and it would mess up the operation counts
  224.                    to include this optimisation at this point in time.
  225.  
  226.                    So... I do it like this until later in the exercise
  227.                    when -D even_faster is defined. */
  228.                 re = (Wn==Wp?1:Wnr(0))*Wr(N/2) - (Wn==Wp?0:Wni(0))*Wi(N/2);
  229.                 im = (Wn==Wp?1:Wnr(0))*Wi(N/2) + (Wn==Wp?0:Wni(0))*Wni(0)*Wr(N/2);
  230.                 Xr(0)   = Wr(0) + re;
  231.                 Xi(0)   = Wi(0) + im;
  232.                 Xr(N/2) = Wr(0) - re;
  233.                 Xi(N/2) = Wi(0) - im;
  234.             }
  235. #endif
  236.         }
  237.     }
  238. }
  239.  
  240. void FastDFS(double *x, double *y, double *w, double *Wp, int N, int skip) {
  241.     FFT(N, x,w,skip, y,skip, Wp,2);
  242. }
  243.  
  244. void FastDFT(double *x, double *y, double *w, double *Wp, int N, int skip) {
  245.     int i;
  246.     double iN = 1./N;
  247.     double *ye;
  248.     FFT(N, y,w,skip, x,skip, Wp+N*2,-2);
  249.     ye = y + N*skip;
  250.     for (ye = y + N*skip; y < ye; y += skip) {
  251.         y[0] *= iN;
  252.         y[1] *= iN;
  253.     }
  254. }
Advertisement
Add Comment
Please, Sign In to add comment
Public Pastes
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!