CLAM-Development
1.1
|
00001 /* 00002 * Copyright (c) 2001-2004 MUSIC TECHNOLOGY GROUP (MTG) 00003 * UNIVERSITAT POMPEU FABRA 00004 * 00005 * 00006 * This program is free software; you can redistribute it and/or modify 00007 * it under the terms of the GNU General Public License as published by 00008 * the Free Software Foundation; either version 2 of the License, or 00009 * (at your option) any later version. 00010 * 00011 * This program is distributed in the hope that it will be useful, 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00014 * GNU General Public License for more details. 00015 * 00016 * You should have received a copy of the GNU General Public License 00017 * along with this program; if not, write to the Free Software 00018 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00019 * 00020 */ 00021 00022 00023 #include "FFT_ooura.hxx" 00024 00025 #include "Assert.hxx" 00026 #include "Audio.hxx" 00027 #include "Spectrum.hxx" 00028 #include "SpectrumConfig.hxx" 00029 #include "CLAM_Math.hxx" 00030 #include "ProcessingFactory.hxx" 00031 00032 namespace CLAM 00033 { 00034 00035 namespace Hidden 00036 { 00037 static const char * metadata[] = { 00038 "key", "FFT_ooura", 00039 "category", "Analysis", 00040 "description", "FFT_ooura", 00041 0 00042 }; 00043 static FactoryRegistrator<ProcessingFactory, FFT_ooura> reg = metadata; 00044 } 00045 00046 bool FFT_ooura::ConcreteConfigure(const ProcessingConfig& c) { 00047 int oldSize=mSize; 00048 FFT_base::ConcreteConfigure(c); 00049 if ( !isPowerOfTwo( mSize ) ) 00050 { 00051 AddConfigErrorMessage("Configure failed: FFT Ooura algorithm only works for input buffers that are a power of two!"); 00052 return false; 00053 } 00054 if (mSize>0) { 00055 if (mSize != oldSize) { 00056 ReleaseMemory(); 00057 SetupMemory(); 00058 } 00059 return true; 00060 } 00061 ReleaseMemory(); 00062 return false; 00063 } 00064 00065 void FFT_ooura::ReleaseMemory() { 00066 if (ip) { delete[] ip; ip = 0; } 00067 if (w) { delete[] w; w = 0; } 00068 } 00069 00070 void FFT_ooura::SetupMemory() { 00071 int ipSize = (int)(2+(1<<(int)(log(mSize/2+0.5)/log(2.0))/2)); 00072 ip = new int[ipSize]; 00073 for (int i=0; i<ipSize; i++) ip[i] = 0; 00074 00075 int wSize = (int)(mSize*5/8-1); 00076 w = new TData[wSize]; 00077 for (int i=0; i<wSize; i++) w[i] = 0; 00078 } 00079 00080 FFT_ooura::FFT_ooura() 00081 : ip(0), w(0) 00082 { 00083 Configure(FFTConfig()); 00084 } 00085 00086 FFT_ooura::FFT_ooura(const FFTConfig &c) throw(ErrDynamicType) 00087 : ip(0), w(0) 00088 { 00089 Configure(c); 00090 }; 00091 00092 FFT_ooura::~FFT_ooura() 00093 { 00094 ReleaseMemory(); 00095 } 00096 00097 bool FFT_ooura::Do() 00098 { 00099 mOutput.GetData().SetSize( mInput.GetSize()/2+1); 00100 bool toReturn = Do(mInput.GetAudio(), mOutput.GetData()); 00101 mInput.Consume(); 00102 mOutput.Produce(); 00103 return toReturn; 00104 } 00105 00106 bool FFT_ooura::Do(const Audio& in, Spectrum &out){ 00107 TData *inbuffer; 00108 00109 CLAM_DEBUG_ASSERT(IsRunning(), 00110 "FFT_ooura: Do(): Not in execution mode"); 00111 CLAM_DEBUG_ASSERT(isPowerOfTwo(mSize), 00112 "FFT_ooura: Do(): Not a power of two"); 00113 00114 out.SetSpectralRange(in.GetSampleRate()/2); 00115 00116 switch(mState) { 00117 case sComplex: 00118 inbuffer = in.GetBuffer().GetPtr(); 00119 // Buffer dump. This is a kludge; the right way to do this 00120 // is using a non-inplace version of rdft (which would 00121 // not reduce performance). 00122 for (int i=0; i<mSize; i++) 00123 fftbuffer[i] = inbuffer[i]; 00124 rdft(mSize, 1, fftbuffer, ip, w); 00125 ToComplex(out); 00126 break; 00127 case sComplexSync: 00128 inbuffer = in.GetBuffer().GetPtr(); 00129 // Buffer dump. This is a kludge; the right way to do this 00130 // is using a non-inplace version of rdft (which would 00131 // not reduce performance). 00132 for (int i=0; i<mSize; i++) 00133 fftbuffer[i]=inbuffer[i]; 00134 rdft(mSize, 1, fftbuffer, ip, w); 00135 ToComplex(out); 00136 out.SynchronizeTo(mComplexflags); 00137 break; 00138 case sOther: 00139 CheckTypes(in,out); 00140 inbuffer = in.GetBuffer().GetPtr(); 00141 // Buffer dump. This is a kludge; the right way to do this 00142 // is using a non-inplace version of rdft (which would 00143 // not reduce performance). 00144 for (int i=0; i<mSize; i++) 00145 fftbuffer[i]=inbuffer[i]; 00146 rdft(mSize, 1, fftbuffer, ip, w); 00147 00148 ToOther(out); 00149 break; 00150 default: 00151 CLAM_ASSERT(false, "FFT_ooura: Do(): Inconsistent state"); 00152 } 00153 00154 return true; 00155 } 00156 00157 00158 void FFT_ooura::ToComplex(Spectrum &out) { 00159 Array<Complex>& outbuffer = out.GetComplexArray(); 00160 00161 outbuffer[0].SetReal(fftbuffer[0]); // Real Values 00162 outbuffer[0].SetImag(0); // Real Values 00163 outbuffer[mSize/2].SetReal(fftbuffer[1]); 00164 outbuffer[mSize/2].SetImag(0); 00165 00166 for (int i=1; i<mSize/2; i++) { 00167 outbuffer[i].SetReal(fftbuffer[2*i]); 00168 outbuffer[i].SetImag(-fftbuffer[2*i+1]); 00169 } 00170 00171 outbuffer.SetSize(mSize/2+1); 00172 } 00173 00174 void FFT_ooura::ToOther(Spectrum &out) 00175 { 00176 if(out.HasComplexArray()) { 00177 ToComplex(out); 00178 out.SynchronizeTo(mComplexflags); 00179 } 00180 else { 00181 SpecTypeFlags flags; 00182 out.GetType(flags); 00183 00184 00185 ToComplex(mComplexSpectrum); 00186 out.SynchronizeTo(mComplexSpectrum); 00187 } 00188 } 00189 00190 00192 // Original Ooura FFT functions, modified to accept TData instead of double 00193 00194 void FFT_ooura::rdft(int n, int isgn, TData *a, int *ip, TData *w) 00195 { 00196 int nw, nc; 00197 TData xi; 00198 00199 nw = ip[0]; 00200 if (n > (nw << 2)) { 00201 nw = n >> 2; 00202 makewt(nw, ip, w); 00203 } 00204 nc = ip[1]; 00205 if (n > (nc << 2)) { 00206 nc = n >> 2; 00207 makect(nc, ip, w + nw); 00208 } 00209 if (isgn >= 0) { 00210 if (n > 4) { 00211 bitrv2(n, ip + 2, a); 00212 cftfsub(n, a, w); 00213 rftfsub(n, a, nc, w + nw); 00214 } else if (n == 4) { 00215 cftfsub(n, a, w); 00216 } 00217 xi = a[0] - a[1]; 00218 a[0] += a[1]; 00219 a[1] = xi; 00220 } else { 00221 a[1] = 0.5 * (a[0] - a[1]); 00222 a[0] -= a[1]; 00223 if (n > 4) { 00224 rftbsub(n, a, nc, w + nw); 00225 bitrv2(n, ip + 2, a); 00226 cftbsub(n, a, w); 00227 } else if (n == 4) { 00228 cftfsub(n, a, w); 00229 } 00230 } 00231 } 00232 00233 // workspace setup functions... 00234 00235 void FFT_ooura::makewt(int nw, int *ip, TData *w) 00236 { 00237 int j, nwh; 00238 TData delta, x, y; 00239 00240 ip[0] = nw; 00241 ip[1] = 1; 00242 if (nw > 2) { 00243 nwh = nw >> 1; 00244 delta = CLAM_atan(1.0) / nwh; 00245 w[0] = 1; 00246 w[1] = 0; 00247 w[nwh] = CLAM_cos(delta * nwh); 00248 w[nwh + 1] = w[nwh]; 00249 if (nwh > 2) { 00250 for (j = 2; j < nwh; j += 2) { 00251 x = CLAM_cos(delta * j); 00252 y = CLAM_sin(delta * j); 00253 w[j] = x; 00254 w[j + 1] = y; 00255 w[nw - j] = y; 00256 w[nw - j + 1] = x; 00257 } 00258 bitrv2(nw, ip + 2, w); 00259 } 00260 } 00261 } 00262 00263 00264 void FFT_ooura::makect(int nc, int *ip, TData *c) 00265 { 00266 int j, nch; 00267 TData delta; 00268 00269 ip[1] = nc; 00270 if (nc > 1) { 00271 nch = nc >> 1; 00272 delta = CLAM_atan(1.0) / nch; 00273 c[0] = CLAM_cos(delta * nch); 00274 c[nch] = 0.5 * c[0]; 00275 for (j = 1; j < nch; j++) { 00276 c[j] = 0.5 * CLAM_cos(delta * j); 00277 c[nc - j] = 0.5 * CLAM_sin(delta * j); 00278 } 00279 } 00280 } 00281 00282 void FFT_ooura::bitrv2(int n, int *ip, TData *a) 00283 { 00284 int j, j1, k, k1, l, m, m2; 00285 TData xr, xi, yr, yi; 00286 00287 ip[0] = 0; 00288 l = n; 00289 m = 1; 00290 while ((m << 3) < l) { 00291 l >>= 1; 00292 for (j = 0; j < m; j++) { 00293 ip[m + j] = ip[j] + l; 00294 } 00295 m <<= 1; 00296 } 00297 m2 = 2 * m; 00298 if ((m << 3) == l) { 00299 for (k = 0; k < m; k++) { 00300 for (j = 0; j < k; j++) { 00301 j1 = 2 * j + ip[k]; 00302 k1 = 2 * k + ip[j]; 00303 xr = a[j1]; 00304 xi = a[j1 + 1]; 00305 yr = a[k1]; 00306 yi = a[k1 + 1]; 00307 a[j1] = yr; 00308 a[j1 + 1] = yi; 00309 a[k1] = xr; 00310 a[k1 + 1] = xi; 00311 j1 += m2; 00312 k1 += 2 * m2; 00313 xr = a[j1]; 00314 xi = a[j1 + 1]; 00315 yr = a[k1]; 00316 yi = a[k1 + 1]; 00317 a[j1] = yr; 00318 a[j1 + 1] = yi; 00319 a[k1] = xr; 00320 a[k1 + 1] = xi; 00321 j1 += m2; 00322 k1 -= m2; 00323 xr = a[j1]; 00324 xi = a[j1 + 1]; 00325 yr = a[k1]; 00326 yi = a[k1 + 1]; 00327 a[j1] = yr; 00328 a[j1 + 1] = yi; 00329 a[k1] = xr; 00330 a[k1 + 1] = xi; 00331 j1 += m2; 00332 k1 += 2 * m2; 00333 xr = a[j1]; 00334 xi = a[j1 + 1]; 00335 yr = a[k1]; 00336 yi = a[k1 + 1]; 00337 a[j1] = yr; 00338 a[j1 + 1] = yi; 00339 a[k1] = xr; 00340 a[k1 + 1] = xi; 00341 } 00342 j1 = 2 * k + m2 + ip[k]; 00343 k1 = j1 + m2; 00344 xr = a[j1]; 00345 xi = a[j1 + 1]; 00346 yr = a[k1]; 00347 yi = a[k1 + 1]; 00348 a[j1] = yr; 00349 a[j1 + 1] = yi; 00350 a[k1] = xr; 00351 a[k1 + 1] = xi; 00352 } 00353 } else { 00354 for (k = 1; k < m; k++) { 00355 for (j = 0; j < k; j++) { 00356 j1 = 2 * j + ip[k]; 00357 k1 = 2 * k + ip[j]; 00358 xr = a[j1]; 00359 xi = a[j1 + 1]; 00360 yr = a[k1]; 00361 yi = a[k1 + 1]; 00362 a[j1] = yr; 00363 a[j1 + 1] = yi; 00364 a[k1] = xr; 00365 a[k1 + 1] = xi; 00366 j1 += m2; 00367 k1 += m2; 00368 xr = a[j1]; 00369 xi = a[j1 + 1]; 00370 yr = a[k1]; 00371 yi = a[k1 + 1]; 00372 a[j1] = yr; 00373 a[j1 + 1] = yi; 00374 a[k1] = xr; 00375 a[k1 + 1] = xi; 00376 } 00377 } 00378 } 00379 } 00380 00381 void FFT_ooura::cftfsub(int n, TData *a, TData *w) 00382 { 00383 int j, j1, j2, j3, l; 00384 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 00385 00386 l = 2; 00387 if (n > 8) { 00388 cft1st(n, a, w); 00389 l = 8; 00390 while ((l << 2) < n) { 00391 cftmdl(n, l, a, w); 00392 l <<= 2; 00393 } 00394 } 00395 if ((l << 2) == n) { 00396 for (j = 0; j < l; j += 2) { 00397 j1 = j + l; 00398 j2 = j1 + l; 00399 j3 = j2 + l; 00400 x0r = a[j] + a[j1]; 00401 x0i = a[j + 1] + a[j1 + 1]; 00402 x1r = a[j] - a[j1]; 00403 x1i = a[j + 1] - a[j1 + 1]; 00404 x2r = a[j2] + a[j3]; 00405 x2i = a[j2 + 1] + a[j3 + 1]; 00406 x3r = a[j2] - a[j3]; 00407 x3i = a[j2 + 1] - a[j3 + 1]; 00408 a[j] = x0r + x2r; 00409 a[j + 1] = x0i + x2i; 00410 a[j2] = x0r - x2r; 00411 a[j2 + 1] = x0i - x2i; 00412 a[j1] = x1r - x3i; 00413 a[j1 + 1] = x1i + x3r; 00414 a[j3] = x1r + x3i; 00415 a[j3 + 1] = x1i - x3r; 00416 } 00417 } else { 00418 for (j = 0; j < l; j += 2) { 00419 j1 = j + l; 00420 x0r = a[j] - a[j1]; 00421 x0i = a[j + 1] - a[j1 + 1]; 00422 a[j] += a[j1]; 00423 a[j + 1] += a[j1 + 1]; 00424 a[j1] = x0r; 00425 a[j1 + 1] = x0i; 00426 } 00427 } 00428 } 00429 00430 void FFT_ooura::cftbsub(int n, TData *a, TData *w) 00431 { 00432 int j, j1, j2, j3, l; 00433 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 00434 00435 l = 2; 00436 if (n > 8) { 00437 cft1st(n, a, w); 00438 l = 8; 00439 while ((l << 2) < n) { 00440 cftmdl(n, l, a, w); 00441 l <<= 2; 00442 } 00443 } 00444 if ((l << 2) == n) { 00445 for (j = 0; j < l; j += 2) { 00446 j1 = j + l; 00447 j2 = j1 + l; 00448 j3 = j2 + l; 00449 x0r = a[j] + a[j1]; 00450 x0i = -a[j + 1] - a[j1 + 1]; 00451 x1r = a[j] - a[j1]; 00452 x1i = -a[j + 1] + a[j1 + 1]; 00453 x2r = a[j2] + a[j3]; 00454 x2i = a[j2 + 1] + a[j3 + 1]; 00455 x3r = a[j2] - a[j3]; 00456 x3i = a[j2 + 1] - a[j3 + 1]; 00457 a[j] = x0r + x2r; 00458 a[j + 1] = x0i - x2i; 00459 a[j2] = x0r - x2r; 00460 a[j2 + 1] = x0i + x2i; 00461 a[j1] = x1r - x3i; 00462 a[j1 + 1] = x1i - x3r; 00463 a[j3] = x1r + x3i; 00464 a[j3 + 1] = x1i + x3r; 00465 } 00466 } else { 00467 for (j = 0; j < l; j += 2) { 00468 j1 = j + l; 00469 x0r = a[j] - a[j1]; 00470 x0i = -a[j + 1] + a[j1 + 1]; 00471 a[j] += a[j1]; 00472 a[j + 1] = -a[j + 1] - a[j1 + 1]; 00473 a[j1] = x0r; 00474 a[j1 + 1] = x0i; 00475 } 00476 } 00477 } 00478 00479 void FFT_ooura::cft1st(int n, TData *a, TData *w) 00480 { 00481 int j, k1, k2; 00482 TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 00483 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 00484 00485 x0r = a[0] + a[2]; 00486 x0i = a[1] + a[3]; 00487 x1r = a[0] - a[2]; 00488 x1i = a[1] - a[3]; 00489 x2r = a[4] + a[6]; 00490 x2i = a[5] + a[7]; 00491 x3r = a[4] - a[6]; 00492 x3i = a[5] - a[7]; 00493 a[0] = x0r + x2r; 00494 a[1] = x0i + x2i; 00495 a[4] = x0r - x2r; 00496 a[5] = x0i - x2i; 00497 a[2] = x1r - x3i; 00498 a[3] = x1i + x3r; 00499 a[6] = x1r + x3i; 00500 a[7] = x1i - x3r; 00501 wk1r = w[2]; 00502 x0r = a[8] + a[10]; 00503 x0i = a[9] + a[11]; 00504 x1r = a[8] - a[10]; 00505 x1i = a[9] - a[11]; 00506 x2r = a[12] + a[14]; 00507 x2i = a[13] + a[15]; 00508 x3r = a[12] - a[14]; 00509 x3i = a[13] - a[15]; 00510 a[8] = x0r + x2r; 00511 a[9] = x0i + x2i; 00512 a[12] = x2i - x0i; 00513 a[13] = x0r - x2r; 00514 x0r = x1r - x3i; 00515 x0i = x1i + x3r; 00516 a[10] = wk1r * (x0r - x0i); 00517 a[11] = wk1r * (x0r + x0i); 00518 x0r = x3i + x1r; 00519 x0i = x3r - x1i; 00520 a[14] = wk1r * (x0i - x0r); 00521 a[15] = wk1r * (x0i + x0r); 00522 k1 = 0; 00523 for (j = 16; j < n; j += 16) { 00524 k1 += 2; 00525 k2 = 2 * k1; 00526 wk2r = w[k1]; 00527 wk2i = w[k1 + 1]; 00528 wk1r = w[k2]; 00529 wk1i = w[k2 + 1]; 00530 wk3r = wk1r - 2 * wk2i * wk1i; 00531 wk3i = 2 * wk2i * wk1r - wk1i; 00532 x0r = a[j] + a[j + 2]; 00533 x0i = a[j + 1] + a[j + 3]; 00534 x1r = a[j] - a[j + 2]; 00535 x1i = a[j + 1] - a[j + 3]; 00536 x2r = a[j + 4] + a[j + 6]; 00537 x2i = a[j + 5] + a[j + 7]; 00538 x3r = a[j + 4] - a[j + 6]; 00539 x3i = a[j + 5] - a[j + 7]; 00540 a[j] = x0r + x2r; 00541 a[j + 1] = x0i + x2i; 00542 x0r -= x2r; 00543 x0i -= x2i; 00544 a[j + 4] = wk2r * x0r - wk2i * x0i; 00545 a[j + 5] = wk2r * x0i + wk2i * x0r; 00546 x0r = x1r - x3i; 00547 x0i = x1i + x3r; 00548 a[j + 2] = wk1r * x0r - wk1i * x0i; 00549 a[j + 3] = wk1r * x0i + wk1i * x0r; 00550 x0r = x1r + x3i; 00551 x0i = x1i - x3r; 00552 a[j + 6] = wk3r * x0r - wk3i * x0i; 00553 a[j + 7] = wk3r * x0i + wk3i * x0r; 00554 wk1r = w[k2 + 2]; 00555 wk1i = w[k2 + 3]; 00556 wk3r = wk1r - 2 * wk2r * wk1i; 00557 wk3i = 2 * wk2r * wk1r - wk1i; 00558 x0r = a[j + 8] + a[j + 10]; 00559 x0i = a[j + 9] + a[j + 11]; 00560 x1r = a[j + 8] - a[j + 10]; 00561 x1i = a[j + 9] - a[j + 11]; 00562 x2r = a[j + 12] + a[j + 14]; 00563 x2i = a[j + 13] + a[j + 15]; 00564 x3r = a[j + 12] - a[j + 14]; 00565 x3i = a[j + 13] - a[j + 15]; 00566 a[j + 8] = x0r + x2r; 00567 a[j + 9] = x0i + x2i; 00568 x0r -= x2r; 00569 x0i -= x2i; 00570 a[j + 12] = -wk2i * x0r - wk2r * x0i; 00571 a[j + 13] = -wk2i * x0i + wk2r * x0r; 00572 x0r = x1r - x3i; 00573 x0i = x1i + x3r; 00574 a[j + 10] = wk1r * x0r - wk1i * x0i; 00575 a[j + 11] = wk1r * x0i + wk1i * x0r; 00576 x0r = x1r + x3i; 00577 x0i = x1i - x3r; 00578 a[j + 14] = wk3r * x0r - wk3i * x0i; 00579 a[j + 15] = wk3r * x0i + wk3i * x0r; 00580 } 00581 } 00582 00583 00584 void FFT_ooura::cftmdl(int n, int l, TData *a, TData *w) 00585 { 00586 int j, j1, j2, j3, k, k1, k2, m, m2; 00587 TData wk1r, wk1i, wk2r, wk2i, wk3r, wk3i; 00588 TData x0r, x0i, x1r, x1i, x2r, x2i, x3r, x3i; 00589 00590 m = l << 2; 00591 for (j = 0; j < l; j += 2) { 00592 j1 = j + l; 00593 j2 = j1 + l; 00594 j3 = j2 + l; 00595 x0r = a[j] + a[j1]; 00596 x0i = a[j + 1] + a[j1 + 1]; 00597 x1r = a[j] - a[j1]; 00598 x1i = a[j + 1] - a[j1 + 1]; 00599 x2r = a[j2] + a[j3]; 00600 x2i = a[j2 + 1] + a[j3 + 1]; 00601 x3r = a[j2] - a[j3]; 00602 x3i = a[j2 + 1] - a[j3 + 1]; 00603 a[j] = x0r + x2r; 00604 a[j + 1] = x0i + x2i; 00605 a[j2] = x0r - x2r; 00606 a[j2 + 1] = x0i - x2i; 00607 a[j1] = x1r - x3i; 00608 a[j1 + 1] = x1i + x3r; 00609 a[j3] = x1r + x3i; 00610 a[j3 + 1] = x1i - x3r; 00611 } 00612 wk1r = w[2]; 00613 for (j = m; j < l + m; j += 2) { 00614 j1 = j + l; 00615 j2 = j1 + l; 00616 j3 = j2 + l; 00617 x0r = a[j] + a[j1]; 00618 x0i = a[j + 1] + a[j1 + 1]; 00619 x1r = a[j] - a[j1]; 00620 x1i = a[j + 1] - a[j1 + 1]; 00621 x2r = a[j2] + a[j3]; 00622 x2i = a[j2 + 1] + a[j3 + 1]; 00623 x3r = a[j2] - a[j3]; 00624 x3i = a[j2 + 1] - a[j3 + 1]; 00625 a[j] = x0r + x2r; 00626 a[j + 1] = x0i + x2i; 00627 a[j2] = x2i - x0i; 00628 a[j2 + 1] = x0r - x2r; 00629 x0r = x1r - x3i; 00630 x0i = x1i + x3r; 00631 a[j1] = wk1r * (x0r - x0i); 00632 a[j1 + 1] = wk1r * (x0r + x0i); 00633 x0r = x3i + x1r; 00634 x0i = x3r - x1i; 00635 a[j3] = wk1r * (x0i - x0r); 00636 a[j3 + 1] = wk1r * (x0i + x0r); 00637 } 00638 k1 = 0; 00639 m2 = 2 * m; 00640 for (k = m2; k < n; k += m2) { 00641 k1 += 2; 00642 k2 = 2 * k1; 00643 wk2r = w[k1]; 00644 wk2i = w[k1 + 1]; 00645 wk1r = w[k2]; 00646 wk1i = w[k2 + 1]; 00647 wk3r = wk1r - 2 * wk2i * wk1i; 00648 wk3i = 2 * wk2i * wk1r - wk1i; 00649 for (j = k; j < l + k; j += 2) { 00650 j1 = j + l; 00651 j2 = j1 + l; 00652 j3 = j2 + l; 00653 x0r = a[j] + a[j1]; 00654 x0i = a[j + 1] + a[j1 + 1]; 00655 x1r = a[j] - a[j1]; 00656 x1i = a[j + 1] - a[j1 + 1]; 00657 x2r = a[j2] + a[j3]; 00658 x2i = a[j2 + 1] + a[j3 + 1]; 00659 x3r = a[j2] - a[j3]; 00660 x3i = a[j2 + 1] - a[j3 + 1]; 00661 a[j] = x0r + x2r; 00662 a[j + 1] = x0i + x2i; 00663 x0r -= x2r; 00664 x0i -= x2i; 00665 a[j2] = wk2r * x0r - wk2i * x0i; 00666 a[j2 + 1] = wk2r * x0i + wk2i * x0r; 00667 x0r = x1r - x3i; 00668 x0i = x1i + x3r; 00669 a[j1] = wk1r * x0r - wk1i * x0i; 00670 a[j1 + 1] = wk1r * x0i + wk1i * x0r; 00671 x0r = x1r + x3i; 00672 x0i = x1i - x3r; 00673 a[j3] = wk3r * x0r - wk3i * x0i; 00674 a[j3 + 1] = wk3r * x0i + wk3i * x0r; 00675 } 00676 wk1r = w[k2 + 2]; 00677 wk1i = w[k2 + 3]; 00678 wk3r = wk1r - 2 * wk2r * wk1i; 00679 wk3i = 2 * wk2r * wk1r - wk1i; 00680 for (j = k + m; j < l + (k + m); j += 2) { 00681 j1 = j + l; 00682 j2 = j1 + l; 00683 j3 = j2 + l; 00684 x0r = a[j] + a[j1]; 00685 x0i = a[j + 1] + a[j1 + 1]; 00686 x1r = a[j] - a[j1]; 00687 x1i = a[j + 1] - a[j1 + 1]; 00688 x2r = a[j2] + a[j3]; 00689 x2i = a[j2 + 1] + a[j3 + 1]; 00690 x3r = a[j2] - a[j3]; 00691 x3i = a[j2 + 1] - a[j3 + 1]; 00692 a[j] = x0r + x2r; 00693 a[j + 1] = x0i + x2i; 00694 x0r -= x2r; 00695 x0i -= x2i; 00696 a[j2] = -wk2i * x0r - wk2r * x0i; 00697 a[j2 + 1] = -wk2i * x0i + wk2r * x0r; 00698 x0r = x1r - x3i; 00699 x0i = x1i + x3r; 00700 a[j1] = wk1r * x0r - wk1i * x0i; 00701 a[j1 + 1] = wk1r * x0i + wk1i * x0r; 00702 x0r = x1r + x3i; 00703 x0i = x1i - x3r; 00704 a[j3] = wk3r * x0r - wk3i * x0i; 00705 a[j3 + 1] = wk3r * x0i + wk3i * x0r; 00706 } 00707 } 00708 } 00709 00710 00711 void FFT_ooura::rftfsub(int n, TData *a, int nc, TData *c) 00712 { 00713 int j, k, kk, ks, m; 00714 TData wkr, wki, xr, xi, yr, yi; 00715 00716 m = n >> 1; 00717 ks = 2 * nc / m; 00718 kk = 0; 00719 for (j = 2; j < m; j += 2) { 00720 k = n - j; 00721 kk += ks; 00722 wkr = 0.5 - c[nc - kk]; 00723 wki = c[kk]; 00724 xr = a[j] - a[k]; 00725 xi = a[j + 1] + a[k + 1]; 00726 yr = wkr * xr - wki * xi; 00727 yi = wkr * xi + wki * xr; 00728 a[j] -= yr; 00729 a[j + 1] -= yi; 00730 a[k] += yr; 00731 a[k + 1] -= yi; 00732 } 00733 } 00734 00735 00736 void FFT_ooura::rftbsub(int n, TData *a, int nc, TData *c) 00737 { 00738 int j, k, kk, ks, m; 00739 TData wkr, wki, xr, xi, yr, yi; 00740 00741 a[1] = -a[1]; 00742 m = n >> 1; 00743 ks = 2 * nc / m; 00744 kk = 0; 00745 for (j = 2; j < m; j += 2) { 00746 k = n - j; 00747 kk += ks; 00748 wkr = 0.5 - c[nc - kk]; 00749 wki = c[kk]; 00750 xr = a[j] - a[k]; 00751 xi = a[j + 1] + a[k + 1]; 00752 yr = wkr * xr + wki * xi; 00753 yi = wkr * xi - wki * xr; 00754 a[j] -= yr; 00755 a[j + 1] = yi - a[j + 1]; 00756 a[k] += yr; 00757 a[k + 1] = yi - a[k + 1]; 00758 } 00759 a[m + 1] = -a[m + 1]; 00760 } 00761 00762 00763 }; // CLAM 00764