/[marvin]/marvin/src/libsip/algo/fft.cc
ViewVC logotype

Diff of /marvin/src/libsip/algo/fft.cc

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.2 by hannibal, Mon Aug 11 22:38:01 2003 UTC revision 1.3 by exos, Thu Sep 11 05:00:21 2003 UTC
# Line 40  Line 40 
40  *******************************************************************/  *******************************************************************/
41  #include "fft.hh"  #include "fft.hh"
42    
43  fft::fft ( ContSignal signal ){  fft::fft ( DiscSignal &signal, bool inv){
44    
45      double angle = -2.0 * PI;
46      unsigned nbi;
47      unsigned ns = signal.length();
48      double *result_r_tmp = (double *)malloc(sizeof(double) * ns);
49      double *result_i_tmp = (double *)malloc(sizeof(double) * ns);
50    
51    
52      if ( !puiss_de_2(ns) )
53        {                              
54          std::cerr << "Not a two power length: "<<ns << std::endl;
55          exit(0);
56        }
57    
58      if (inv)
59        angle *= -1;
60      
61      nbi = NumberOfBitsNeeded (ns);
62    
63      int j, i, k, n;
64      double tr, ti;
65      
66      for (i = 0; i < ns; i++ )
67        {
68          j = ReverseBits (i, nbi);
69          result_r_tmp[j] = signal[i];
70          result_i_tmp[j] = 0.0;
71          }
72      
73      result_r = new DiscSignal(result_r_tmp, ns);
74      result_i = new DiscSignal(result_i_tmp, ns);
75      
76      
77      unsigned ben, bsz;
78      for ( ben = 1, bsz = 2; bsz <= ns; bsz <<= 1 )
79        {
80            double delta_angle = angle / (double)bsz;
81            double sm2 = sin ( -2 * delta_angle );
82            double sm1 = sin ( -delta_angle );
83            double cm2 = cos ( -2 * delta_angle );
84            double cm1 = cos ( -delta_angle );
85            double w = 2 * cm1;
86            double ar[3], ai[3];
87    
88            for ( i=0; i < ns; i += bsz )
89            {
90                ar[2] = cm2;
91                ar[1] = cm1;
92    
93                ai[2] = sm2;
94                ai[1] = sm1;
95    
96                for ( j=i, n=0; n < ben; j++, n++ )
97                {
98                    ar[0] = w*ar[1] - ar[2];
99                    ar[2] = ar[1];
100                    ar[1] = ar[0];
101    
102                    ai[0] = w*ai[1] - ai[2];
103                    ai[2] = ai[1];
104                    ai[1] = ai[0];
105    
106                    k = j + ben;
107                    tr = ar[0]*(*result_r)[k] - ai[0]*(*result_i)[k];
108                    ti = ar[0]*(*result_i)[k] + ai[0]*(*result_r)[k];
109    
110                    (*result_r)[k] = (*result_r)[j] - tr;
111                    (*result_i)[k] = (*result_i)[j] - ti;
112    
113                    (*result_r)[j] += tr;
114                    (*result_i)[j] += ti;
115                }
116            }
117    
118            ben = bsz;
119        }
120    
121        if (inv)
122        {
123            double denom = (double)ns;
124    
125            for ( i=0; i < ns; i++ )
126            {
127                (*result_r)[i] /= denom;
128                (*result_i)[i] /= denom;
129            }
130        }
131        
132  }  }
133    
134  fft::fft ( DiscSignal signal ){  DiscSignal &fft::get_magnitude()
135    {
136      DiscSignal *result;
137    
138      result = new DiscSignal((double *)malloc(sizeof(double) * result_i->length()), result_i->length());
139      for (int i = 0; i < result_i->length(); ++i)
140        (*result)[i] = sqrt ((*result_i)[i] * (*result_i)[i] + (*result_r)[i] * (*result_r)[i]);
141    
142      return *result;
143  }  }
144    
145  fft::fft ( ContLogSignal signal ){  
146    unsigned fft::NumberOfBitsNeeded ( unsigned PowerOfTwo )
147    {
148        unsigned i;
149    
150        if ( PowerOfTwo < 2 )
151        {
152          std::cerr<<"Not a good signal."<<std::endl;
153          exit(1);
154        }
155    
156        for ( i=0; ; i++ )
157        {
158            if ( PowerOfTwo & (1 << i) )
159                return i;
160        }
161    }
162    
163    
164    unsigned fft::ReverseBits ( unsigned index, unsigned NumBits )
165    {
166        unsigned i, rev;
167    
168        for ( i=rev=0; i < NumBits; i++ )
169        {
170            rev = (rev << 1) | (index & 1);
171            index >>= 1;
172        }
173    
174        return rev;
175    }
176    
177    bool fft::puiss_de_2 ( unsigned x )
178    {
179        if ( x < 2 )
180            return false;
181    
182        if ( x & (x-1) )      
183            return false;
184    
185        return true;
186    }
187    
188    
189    
190    fft::fft ( ContSignal signal ){
191  }  }
192    
193  fft::fft ( DiscLogSignal signal ){  fft::fft ( DiscLogSignal signal ){
194  }  }
195    
196    fft::fft ( ContLogSignal signal ){
197    }
198    

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.3

savannah-hackers-public@gnu.org
ViewVC Help
Powered by ViewVC 1.1.26