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 |
|
|