~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

Linux Cross Reference
Tina6/tina-libs/tina/math/mathNum_fourier.c

Version: ~
Architecture: ~ [ i386 ] ~ [ alpha ] ~ [ m68k ] ~ [ mips ] ~ [ ppc ] ~ [ sparc ] ~ [ sparc64 ] ~

  1 /**********
  2  *
  3  * This file is part of the TINA Open Source Image Analysis Environment
  4  * henceforth known as TINA
  5  *
  6  * TINA is free software; you can redistribute it and/or modify
  7  * it under the terms of the GNU General Public License as
  8  * published by the Free Software Foundation.
  9  *
 10  * TINA is distributed in the hope that it will be useful,
 11  * but WITHOUT ANY WARRANTY; without even the implied warranty of
 12  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 13  * GNU General Public License for more details.
 14  *
 15  * You should have received a copy of the GNU General Public License
 16  * along with TINA; if not, write to the Free Software Foundation, Inc.,
 17  * 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 18  *
 19  * ANY users of TINA who require exemption from the existing licence must
 20  * negotiate a new licence with Dr. Neil.A.Thacker, the sole agent for
 21  * the University of Manchester.
 22  *
 23  **********
 24  *
 25  * Program :    TINA
 26  * File    :  $Source: /home/tina/cvs/tina-libs/tina/math/mathNum_fourier.c,v $
 27  * Date    :  $Date: 2005/01/23 19:10:21 $
 28  * Version :  $Revision: 1.6 $
 29  * CVS Id  :  $Id: mathNum_fourier.c,v 1.6 2005/01/23 19:10:21 paul Exp $
 30  *
 31  * Author  :  Legacy TINA
 32  *
 33  * Notes : Fast Fourier Transform
 34  *
 35  *********
 36 */
 37 /** 
 38  *  @file
 39  *  @brief  Fast Fourier Transform functions.   
 40  *
 41  *  See Numerical Recipes (3rd ed) Ch 12.2 pp 504-508 for more details regarding the algorithm.
 42  * 
 43 */
 44 
 45 #include "mathNum_fourier.h"
 46 
 47 #if HAVE_CONFIG_H
 48 #include <config.h>
 49 #endif
 50 
 51 #include <math.h>
 52 #include <tina/sys/sysDef.h>
 53 #include <tina/sys/sysPro.h>
 54 #include <tina/math/math_GenDef.h>
 55 
 56 void            fft_cmplx_inplace(Complex * z, int n)   /** assumes n is a power of 2 **/
 57 
 58 
 59 {
 60         int             i;
 61         double         *y = dvector_alloc(0, 2 * n);
 62 
 63         for (i = 0; i < n; i++)
 64         {
 65                 y[2 * i] = cmplx_re(z[i]);
 66                 y[2 * i + 1] = cmplx_im(z[i]);
 67         }
 68 
 69         fourier(y - 1, n, 1);
 70 
 71 
 72         for (i = 0; i < n; i++)
 73         {
 74                 cmplx_re(z[i]) = y[2 * i];
 75                 cmplx_im(z[i]) = y[2 * i + 1];
 76         }
 77 
 78         dvector_free(y, 0);
 79 }
 80 
 81 void            fft_inverse_cmplx_inplace(Complex * z, int n)   /** assumes n is a power of 2 **/
 82 
 83 
 84 {
 85         int             i;
 86         double         *y = dvector_alloc(0, 2 * n);
 87 
 88         for (i = 0; i < n; i++)
 89         {
 90                 y[2 * i] = cmplx_re(z[i]);
 91                 y[2 * i + 1] = cmplx_im(z[i]);
 92         }
 93 
 94         fourier(y - 1, n, -1);
 95 
 96         for (i = 0; i < n; i++)
 97         {
 98                 cmplx_re(z[i]) = y[2 * i] / n;
 99                 cmplx_im(z[i]) = y[2 * i + 1] / n;
100         }
101 
102         dvector_free(y, 0);
103 }
104 
105 /** assumes complex data is interlaced in data[1] ... data[2n] **/
106 
107 void            fourier(double *data, int nn, int isign)
108 {
109         int             n, mmax, m, j, istep, i;
110         double          wtemp, wr, wpr, wpi, wi, theta;
111         double          tempr, tempi;
112 
113         n = nn << 1;
114         j = 1;
115         for (i = 1; i < n; i += 2)
116         {
117                 if (j > i)
118                 {
119                         SWAP(double, data[j], data[i]);
120                         SWAP(double, data[j + 1], data[i + 1]);
121                 }
122                 m = n >> 1;
123                 while (m >= 2 && j > m)
124                 {
125                         j -= m;
126                         m >>= 1;
127                 }
128                 j += m;
129         }
130         mmax = 2;
131         while (n > mmax)
132         {
133                 istep = 2 * mmax;
134                 theta = TWOPI / (isign * mmax);
135                 wtemp = sin(0.5 * theta);
136                 wpr = -2.0 * wtemp * wtemp;
137                 wpi = sin(theta);
138                 wr = 1.0;
139                 wi = 0.0;
140                 for (m = 1; m < mmax; m += 2)
141                 {
142                         for (i = m; i <= n; i += istep)
143                         {
144                                 j = i + mmax;
145                                 tempr = wr * data[j] - wi * data[j + 1];
146                                 tempi = wr * data[j + 1] + wi * data[j];
147                                 data[j] = data[i] - tempr;
148                                 data[j + 1] = data[i + 1] - tempi;
149                                 data[i] += tempr;
150                                 data[i + 1] += tempi;
151                         }
152                         wr = (wtemp = wr) * wpr - wi * wpi + wr;
153                         wi = wi * wpr + wtemp * wpi + wi;
154                 }
155                 mmax = istep;
156         }
157 }
158 

~ [ source navigation ] ~ [ diff markup ] ~ [ identifier search ] ~ [ freetext search ] ~ [ file search ] ~

This page was automatically generated by the LXR engine.
Visit the LXR main site for more information.