CSL 5.2

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

FFT_Wrapper.cpp
Go to the documentation of this file.
1 //
2 // FFT_Wrapper.cpp -- wrapper class for FFTs that hides implementation details
3 // This file includes the 3 standard concrete subclasses:
4 // FFTW (assumes fftw3f),
5 // RealFFT (RealFFT code included in CSL), and
6 // KISS FFT (included, untested).
7 //
8 // See the copyright notice and acknowledgment of authors in the file COPYRIGHT
9 //
10
11 #include "FFT_Wrapper.h"
12 #include "math.h"
13 #include <string.h> // for memcpy
14
15 using namespace csl;
16
17 //-------------------------------------------------------------------------------------------------//
18
19 #ifdef USE_FFTW // the FFTW-based version
20
21 // FFTW_Wrapper = FFTW-based concrete implementation
22
23 FFTW_Wrapper::FFTW_Wrapper(unsigned size, CSL_FFTType type, CSL_FFTDir forward)
24 : Abst_FFT_W(size, type, forward) {
25 // allocate buffers
26 mSampBuf = (SampleBuffer) fftwf_malloc(size * sizeof(sample));
27 mSpectBuf = (fftwf_complex *) fftwf_malloc(mCSize * sizeof(fftwf_complex));
28
29 if (mDirection == CSL_FFT_FORWARD) { // create FFT plans
30 // (int n, float *in, fftw_complex *out, unsigned flags);
31 mPlan = fftwf_plan_dft_r2c_1d(size, mSampBuf, mSpectBuf, FFTWF_FLAGS);
32 } else { // inverse FFT
33 // int n, fftw_complex *in, float *out, unsigned flags);
34 mPlan = fftwf_plan_dft_c2r_1d(size, mSpectBuf, mSampBuf, FFTWF_FLAGS);
35 }
36 // logMsg("FFTW %d", size);
37 }
38
39 FFTW_Wrapper::~FFTW_Wrapper() {
40 fftwf_destroy_plan(mPlan);
41 fftwf_free(mSampBuf);
42 fftwf_free(mSpectBuf);
43 }
44
45 /// run the transform
46
47 void FFTW_Wrapper::nextBuffer(Buffer & in, Buffer & out) throw (CException) {
48
49 if (mDirection == CSL_FFT_FORWARD) { // mDirection == CSL_FFT_FORWARD
50 // copy input into sample buffer
51 memcpy(mSampBuf, in.buffer(0), mSize * sizeof(sample));
52 // printf("\t%x - %x\n", in.buffer(0), out.buffer(0));
53
54 fftwf_execute(mPlan); //// GO ////
55
56 SampleBuffer ioPtr = out.buffer(0); // out buffer
57 fftwf_complex * spPtr = mSpectBuf; // spectrum ptr
58
59 if (mType == CSL_FFT_REAL) { // real: write magnitude to mono output
60 for (unsigned j = 0; j < mCSize; j++, spPtr++)
61 *ioPtr++ = hypotf((*spPtr)[0], (*spPtr)[0]);
62
63 } else if (mType == CSL_FFT_COMPLEX) { // copy complex ptr to output
64 memcpy(out.buffer(0), mSpectBuf, (mCSize * sizeof(fftwf_complex)));
65
66 } else if (mType == CSL_FFT_MAGPHASE) { // write mag/phase to buffer[0]/[1]
67 SampleBuffer inOutPh = out.buffer(1);
68 for (unsigned j = 0; j < in.mNumFrames; j++) {
69 // fprintf(stderr, "re:%f cx:%f ", (*spPtr)[0], (*spPtr)[1]);
70 *ioPtr++ = hypotf((*spPtr)[0], (*spPtr)[1]); // write magnitude
71 spPtr++;
72 if ((*spPtr)[0] == 0.0f) {
73 if ((*spPtr)[1] >= 0.0f)
74 *inOutPh++ = CSL_PIHALF;
75 else
76 *inOutPh++ = CSL_PI + CSL_PIHALF;
77 } else {
78 *inOutPh++ = atan((*spPtr)[1] / (*spPtr)[0]); // write phase
79 }
80 }
81 }
82 } else { // mDirection == CSL_FFT_INVERSE
83 // copy data into spectrum
84 memcpy(mSpectBuf, in.buffer(0), (mCSize * sizeof(fftwf_complex)));
85
86 fftwf_execute(mPlan); // GO
87 // copy real output
88 memcpy(out.buffer(0), mSampBuf, (mSize * sizeof(sample)));
89 }
90 }
91
92 #endif // FFTW
93
94 //-------------------------------------------------------------------------------------------------//
95
96 #ifdef USE_FFTREAL // the FFTReal-based version
97
98 // FFTR_Wrapper = FFTReal-based concrete implementation
99
100 FFTR_Wrapper::FFTR_Wrapper(unsigned size, CSL_FFTType type, CSL_FFTDir forward)
101 : Abst_FFT_W(size, type, forward), mFFT(size) {
102 SAFE_MALLOC(mTempBuf, sample, size + 1);
103 // logMsg("FFTReal %d", size);
104 }
105
106 FFTR_Wrapper::~FFTR_Wrapper() {
107 SAFE_FREE(mTempBuf);
108 }
109
110 // execute = run the transform
111
112 void FFTR_Wrapper::nextBuffer(Buffer & in, Buffer & out) throw (CException) {
113 if (mDirection == CSL_FFT_FORWARD) { // mDirection == CSL_FFT_FORWARD
114 SampleBuffer ioPtr = in.buffer(0); // set input data ptr
115 // printf("\t%x - %x\n", in.buffer(0), out.buffer(0));
116
117 mFFT.do_fft(mTempBuf, ioPtr); // perform FFT
118
119 // do_fft (flt_t f [], const flt_t x []) -- this is the fcn signature
120 // - x: pointer on the source array (time).
121 // - f: pointer on the destination array (frequencies).
122 // f [0...length(x)/2] = real values,
123 // f [length(x)/2+1...length(x)-1] = imaginary values of coeff 1...length(x)/2-1.
124
125 if (mType == CSL_FFT_COMPLEX) { // raw: copy complex points to output
126 SampleComplexPtr cxPtr = (SampleComplexPtr) out.buffer(0);
127 SampleBuffer rPtr = mTempBuf;
128 SampleBuffer iPtr = mTempBuf + (mSize / 2) ; // + 1;
129 float normFactor = 1.0 / sqrt((double) mSize);
130 for (unsigned j = 0; j < mSize / 2; j++) {
131 ComplexPtr cplx = cxPtr[j];
132 cx_r(cplx) = *rPtr++ * normFactor;
133 cx_i(cplx) = *iPtr++ * normFactor;
134 }
135 }
136 else if (mType == CSL_FFT_REAL) { // real: write magnitude to mono output
137 ioPtr = out.buffer(0); // output pointer
138 SampleBuffer rPtr = mTempBuf;
139 SampleBuffer iPtr = mTempBuf + (mSize / 2) ; // + 1;
140 for (unsigned j = 0; j < mSize / 2; j++)
141 *ioPtr++ = hypotf(*rPtr++, *iPtr++);
142 }
143 else if (mType == CSL_FFT_MAGPHASE) { // complex: write mag/phase to buffer[0]/[1]
144 ioPtr = out.buffer(0); // output pointer
145 SampleBuffer rPtr = mTempBuf;
146 SampleBuffer iPtr = mTempBuf + (mSize / 2);
147 SampleBuffer inOutPh = out.buffer(1);
148 for (unsigned j = 0; j < mSize / 2; j++) {
149 *ioPtr++ = hypotf(*rPtr, *iPtr);
150 if (*rPtr == 0.0f) {
151 if (*iPtr >= 0.0f)
152 *inOutPh++ = CSL_PIHALF;
153 else
154 *inOutPh++ = CSL_PI + CSL_PIHALF;
155 } else {
156 *inOutPh++ = atan(*iPtr / *rPtr);
157 }
158 rPtr++;
159 iPtr++;
160 } // end of loop
161 }
162 } else { // mDirection == CSL_FFT_INVERSE
163 // assume CSL_FFT_COMPLEX format
164 // copy complex spectrum to un-packed IFFT format
165 SampleComplexPtr ioPtr = (SampleComplexPtr) in.buffer(0);
166 SampleBuffer rPtr = mTempBuf; // copy data to FFTReal format
167 SampleBuffer iPtr = mTempBuf + (mSize / 2);
168
169 for (unsigned j = 0; j < mSize / 2; j++) { // loop to unpack complex array
170 ComplexPtr cplx = ioPtr[j];
171 *rPtr++ = cx_r(cplx);
172 *iPtr++ = cx_i(cplx);
173 }
174 SampleBuffer oPtr = out.buffer(0); // output pointer
175
176 mFFT.do_ifft(mTempBuf, oPtr);
177
178 // do_ifft (const flt_t f [], flt_t x [])
179 // - f: pointer on the source array (frequencies).
180 // f [0...length(x)/2] = real values,
181 // f [length(x)/2+1...length(x)] = imaginary values of coeff 1...length(x)-1.
182 // - x: pointer on the destination array (time).
183
184 }
185 }
186
187 #endif
188
189 //-------------------------------------------------------------------------------------------------//
190
191
192 #ifdef USE_KISSFFT
193
194 // KISSFFT_Wrapper = KISS FFT-based concrete implementation
195
196 KISSFFT_Wrapper::KISSFFT_Wrapper(unsigned size, CSL_FFTType type, CSL_FFTDir forward)
197 : Abst_FFT_W(size, type, forward) {
198 int dir = (forward == CSL_FFT_FORWARD) ? 0 : 1;
199 mFFT = kiss_fft_alloc(size, dir, NULL, NULL);
200 SAFE_MALLOC(inBuf, SampleComplex, size);
201 SAFE_MALLOC(outBuf, SampleComplex, size);
202 }
203
204 KISSFFT_Wrapper::~KISSFFT_Wrapper() {
205 free(mFFT);
206 kiss_fft_cleanup();
207 }
208
209 // run the transform between in and out
210
211 void KISSFFT_Wrapper::nextBuffer(Buffer & in, Buffer & out) throw (CException) {
212
213 if (mDirection == CSL_FFT_FORWARD) { // mDirection == CSL_FFT_FORWARD
214 SampleBuffer ioPtr = in.buffer(0); // input data ptr
215 SampleComplexPtr cxPtr = inBuf;
216 for (int j = 0; j < mSize; j++) { // loop to pack complex array
217 *cxPtr[0] = *ioPtr++;
218 cxPtr++;
219 }
220
221 // kiss_fft(kiss_fft_cfg cfg, const kiss_fft_cpx *fin, kiss_fft_cpx *fout)
222 // Perform an FFT on a complex input buffer.
223 // for a forward FFT,
224 // fin should be f[0] , f[1] , ... ,f[nfft-1]
225 // fout will be F[0] , F[1] , ... ,F[nfft-1]
226 // Note that each element is complex and can be accessed like f[k].r and f[k].i
227
228 kiss_fft(mFFT, (const kiss_fft_cpx *) inBuf, (kiss_fft_cpx *) outBuf);
229
230 ioPtr = out.buffer(0); // output pointer
231 if (mType == CSL_FFT_REAL) { // real: write magnitude to mono output
232 SampleBuffer rPtr = mTempBuf;
233 SampleBuffer iPtr = mTempBuf + (mSize / 2);
234 *ioPtr++ = *rPtr++;
235 *ioPtr++;
236 for (int j = 1; j < mSize; j++)
237 *ioPtr++ = hypotf(*rPtr++, *iPtr++);
238 } else if (mType == CSL_FFT_COMPLEX) { // raw: copy complex points to output
239 memcpy(ioPtr, outBuf, mSize * sizeof(SampleComplex));
240 }
241 } else { // mDirection == CSL_FFT_INVERSE
242
243 kiss_fft(mFFT, (const kiss_fft_cpx *) in.buffer(0), (kiss_fft_cpx *) outBuf);
244
245 SampleComplexPtr cxPtr = outBuf;
246 SampleBuffer ioPtr = out.buffer(0);
247 for (int j = 0; j < mSize; j++) // loop to unpack complex array
248 *ioPtr++ = cxPtr[j][0];
249 }
250 }
251
252 #endif