home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 22 gnu / 22-gnu.zip / mesch12a.zip / fft.c < prev    next >
C/C++ Source or Header  |  1994-01-13  |  4KB  |  145 lines
  1.  
  2. /**************************************************************************
  3. **
  4. ** Copyright (C) 1993 David E. Steward & Zbigniew Leyk, all rights reserved.
  5. **
  6. **                 Meschach Library
  7. ** 
  8. ** This Meschach Library is provided "as is" without any express 
  9. ** or implied warranty of any kind with respect to this software. 
  10. ** In particular the authors shall not be liable for any direct, 
  11. ** indirect, special, incidental or consequential damages arising 
  12. ** in any way from use of the software.
  13. ** 
  14. ** Everyone is granted permission to copy, modify and redistribute this
  15. ** Meschach Library, provided:
  16. **  1.  All copies contain this copyright notice.
  17. **  2.  All modified copies shall carry a notice stating who
  18. **      made the last modification and the date of such modification.
  19. **  3.  No charge is made for this software or works derived from it.  
  20. **      This clause shall not be construed as constraining other software
  21. **      distributed on the same medium as this software, nor is a
  22. **      distribution fee considered a charge.
  23. **
  24. ***************************************************************************/
  25.  
  26.  
  27. /*
  28.     Fast Fourier Transform routine
  29.     Loosely based on the Fortran routine in Rabiner & Gold's
  30.     "Digital Signal Processing"
  31. */
  32.  
  33. static char rcsid[] = "$Id: fft.c,v 1.3 1994/01/13 05:45:33 des Exp $";
  34.  
  35. #include        <stdio.h>
  36. #include        <math.h>
  37. #include        "matrix.h"
  38. #include        "matrix2.h"
  39.  
  40.  
  41. /* fft -- d.i.t. fast Fourier transform 
  42.         -- radix-2 FFT only
  43.         -- vector extended to a power of 2 */
  44. void    fft(x_re,x_im)
  45. VEC     *x_re, *x_im;
  46. {
  47.     int         i, ip, j, k, li, n, length;
  48.     Real      *xr, *xi;
  49.     Real    theta, pi = 3.1415926535897932384;
  50.     Real      w_re, w_im, u_re, u_im, t_re, t_im;
  51.     Real      tmp, tmpr, tmpi;
  52.  
  53.     if ( ! x_re || ! x_im )
  54.         error(E_NULL,"fft");
  55.     if ( x_re->dim != x_im->dim )
  56.         error(E_SIZES,"fft");
  57.  
  58.     n = 1;
  59.     while ( x_re->dim > n )
  60.         n *= 2;
  61.     x_re = v_resize(x_re,n);
  62.     x_im = v_resize(x_im,n);
  63.     printf("# fft: x_re =\n");  v_output(x_re);
  64.     printf("# fft: x_im =\n");  v_output(x_im);
  65.     xr   = x_re->ve;
  66.     xi   = x_im->ve;
  67.  
  68.     /* Decimation in time (DIT) algorithm */
  69.     j = 0;
  70.     for ( i = 0; i < n-1; i++ )
  71.     {
  72.         if ( i < j )
  73.         {
  74.             tmp = xr[i];
  75.             xr[i] = xr[j];
  76.             xr[j] = tmp;
  77.             tmp = xi[i];
  78.             xi[i] = xi[j];
  79.             xi[j] = tmp;
  80.         }
  81.         k = n / 2;
  82.         while ( k <= j )
  83.         {
  84.             j -= k;
  85.             k /= 2;
  86.         }
  87.         j += k;
  88.     }
  89.  
  90.     /* Actual FFT */
  91.     for ( li = 1; li < n; li *= 2 )
  92.     {
  93.         length = 2*li;
  94.         theta  = pi/li;
  95.         u_re   = 1.0;
  96.         u_im   = 0.0;
  97.         if ( li == 1 )
  98.         {
  99.             w_re = -1.0;
  100.             w_im =  0.0;
  101.         }
  102.         else if ( li == 2 )
  103.         {
  104.             w_re =  0.0;
  105.             w_im =  1.0;
  106.         }
  107.         else
  108.         {
  109.             w_re = cos(theta);
  110.             w_im = sin(theta);
  111.         }
  112.         for ( j = 0; j < li; j++ )
  113.         {
  114.             for ( i =  j; i < n; i += length )
  115.             {
  116.                 ip = i + li;
  117.                 /* step 1 */
  118.                 t_re = xr[ip]*u_re - xi[ip]*u_im;
  119.                 t_im = xr[ip]*u_im + xi[ip]*u_re;
  120.                 /* step 2 */
  121.                 xr[ip] = xr[i]  - t_re;
  122.                 xi[ip] = xi[i]  - t_im;
  123.                 /* step 3 */
  124.                 xr[i] += t_re;
  125.                 xi[i] += t_im;
  126.             }
  127.             tmpr = u_re*w_re - u_im*w_im;
  128.             tmpi = u_im*w_re + u_re*w_im;
  129.             u_re = tmpr;
  130.             u_im = tmpi;
  131.         }
  132.     }
  133. }
  134.  
  135. /* ifft -- inverse FFT using the same interface as fft() */
  136. void    ifft(x_re,x_im)
  137. VEC    *x_re, *x_im;
  138. {
  139.     /* we just use complex conjugates */
  140.  
  141.     sv_mlt(-1.0,x_im,x_im);
  142.     fft(x_re,x_im);
  143.     sv_mlt(-1.0/((double)(x_re->dim)),x_im,x_im);
  144. }
  145.