annotate fft.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix
author Chris Cannam
date Mon, 05 Sep 2011 15:22:35 +0100
parents 977f541d6683
children
rev   line source
xue@11 1 /*
xue@11 2 Harmonic sinusoidal modelling and tools
xue@11 3
xue@11 4 C++ code package for harmonic sinusoidal modelling and relevant signal processing.
xue@11 5 Centre for Digital Music, Queen Mary, University of London.
xue@11 6 This file copyright 2011 Wen Xue.
xue@11 7
xue@11 8 This program is free software; you can redistribute it and/or
xue@11 9 modify it under the terms of the GNU General Public License as
xue@11 10 published by the Free Software Foundation; either version 2 of the
xue@11 11 License, or (at your option) any later version.
xue@11 12 */
xue@1 13 //---------------------------------------------------------------------------
xue@1 14
Chris@2 15 #include <string.h>
xue@1 16 #include <stdlib.h>
xue@1 17 #include "fft.h"
xue@1 18
Chris@5 19 /** \file fft.h */
Chris@5 20
xue@1 21 //---------------------------------------------------------------------------
Chris@5 22 /**
xue@1 23 function Atan2: (0, 0)-safe atan2
xue@1 24
xue@1 25 Returns 0 is x=y=0, atan2(x,y) otherwise.
xue@1 26 */
xue@1 27 double Atan2(double y, double x)
xue@1 28 {
xue@1 29 if (x==0 && y==0) return 0;
xue@1 30 else return atan2(y, x);
xue@1 31 }//Atan2
xue@1 32
Chris@5 33 /**
xue@11 34 function Log2: Log2
xue@11 35
xue@11 36 Returns x for Log2(2^x) if x is integer.
xue@11 37 */
xue@11 38 int Log2(int x)
xue@11 39 {
xue@11 40 return floor(log(x)/log(2)+0.5);
xue@11 41 }//Log2
xue@11 42
xue@11 43 /**
xue@1 44 function BitInv: inverse bit order of Value within an $Order-bit expression.
xue@1 45
xue@1 46 In: integer Value smaller than 2^Order
xue@1 47
xue@1 48 Returns an integer whose lowest Order bits are the lowest Order bits of Value in reverse order.
xue@1 49 */
xue@1 50 int BitInv(int Value, int Order)
xue@1 51 {
xue@1 52 int Result;
xue@1 53 Result=0;
xue@1 54 for (int i=0;i<Order;i++)
xue@1 55 {
xue@1 56 Result=(Result<<1)+(Value&0x00000001);
xue@1 57 Value=Value>>1;
xue@1 58 }
xue@1 59 return Result;
xue@1 60 }//BitInv
xue@1 61
Chris@5 62 /**
xue@1 63 function SetTwiddleFactors: fill w[N/2] with twiddle factors used in N-point complex FFT.
xue@1 64
xue@1 65 In: N
xue@1 66 Out: array w[N/2] containing twiddle factors
xue@1 67
xue@1 68 No return value.
xue@1 69 */
xue@1 70 void SetTwiddleFactors(int N, cdouble* w)
xue@1 71 {
xue@1 72 double ep=-M_PI*2/N;
xue@1 73 for (int i=0; i<N/2; i++)
xue@1 74 {
xue@1 75 double tmp=ep*i;
xue@1 76 w[i].x=cos(tmp), w[i].y=sin(tmp);
xue@1 77 }
xue@1 78 }//SetTwiddleFactors
xue@1 79
xue@1 80 //---------------------------------------------------------------------------
Chris@5 81 /**
xue@1 82 function CFFTCbii: basic complex DIF-FFT module, applied after bit-inversed ordering of inputs
xue@1 83
xue@1 84 In: Order: integer, equals log2(Wid)
xue@1 85 W[Wid/2]: twiddle factors
xue@1 86 X[Wid]: complex waveform
xue@1 87 Out: X[Wid]: complex spectrum
xue@1 88
xue@1 89 No return value.
xue@1 90 */
xue@1 91 void CFFTCbii(int Order, cdouble* W, cdouble* X)
xue@1 92 {
xue@1 93 int i, j, k, ElemsPerGroup, Groups, X0, X1, X2;
xue@1 94 cdouble Temp;
xue@1 95 for (i=0; i<Order; i++)
xue@1 96 {
xue@1 97 ElemsPerGroup=1<<i;
xue@1 98 Groups=1<<(Order-i-1);
xue@1 99 X0=0;
xue@1 100 for (j=0; j<Groups; j++)
xue@1 101 {
xue@1 102 for (k=0; k<ElemsPerGroup; k++)
xue@1 103 {
xue@1 104 int kGroups=k*Groups;
xue@1 105 X1=X0+k;
xue@1 106 X2=X1+ElemsPerGroup;
xue@1 107 //X(X2)<-X(X2)*W
xue@1 108 Temp.x=X[X2].x*W[kGroups].x-X[X2].y*W[kGroups].y,
xue@1 109 X[X2].y=X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1 110 X[X2].x=Temp.x;
xue@1 111 Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1 112 X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1 113 X[X1]=Temp;
xue@1 114 }
xue@1 115 X0+=ElemsPerGroup*2;
xue@1 116 }
xue@1 117 }
xue@1 118 }//CFFTCbii
xue@1 119
Chris@5 120 /**
xue@1 121 function CFFTC: in-place complex FFT
xue@1 122
xue@1 123 In: Order: integer, equals log2(Wid)
xue@1 124 W[Wid/2]: twiddle factors
xue@1 125 X[Wid]: complex waveform
xue@1 126 bitinv[Wid]: bit-inversion table
xue@1 127 Out: X[Wid]: complex spectrum
xue@1 128
xue@1 129 No return value.
xue@1 130 */
xue@1 131 void CFFTC(int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 132 {
xue@1 133 int N=1<<Order, i, jj;
xue@1 134 cdouble Temp;
xue@1 135 int* bitinv1=bitinv;
xue@1 136 if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1 137 for (i=1; i<N-1; i++)
xue@1 138 {
xue@1 139 jj=bitinv[i];
xue@1 140 if (i<jj)
xue@1 141 {
xue@1 142 Temp=X[i];
xue@1 143 X[i]=X[jj];
xue@1 144 X[jj]=Temp;
xue@1 145 }
xue@1 146 }
xue@1 147 if (!bitinv1) free(bitinv);
xue@1 148 CFFTCbii(Order, W, X);
xue@1 149 }//CFFTC
xue@1 150
Chris@5 151 /**
xue@1 152 function CFFTC: complex FFT
xue@1 153
xue@1 154 In: Input[Wid]: complex waveform
xue@1 155 Order: integer, equals log2(Wid)
xue@1 156 W[Wid/2]: twiddle factors
xue@1 157 bitinv[Wid]: bit-inversion table
xue@1 158 Out:X[Wid]: complex spectrum
xue@1 159 Amp[Wid]: amplitude spectrum
xue@1 160 Arg[Wid]: phase spectrum
xue@1 161
xue@1 162 No return value.
xue@1 163 */
xue@1 164 void CFFTC(cdouble* Input, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 165 {
xue@1 166 int i, N=1<<Order;
xue@1 167
xue@1 168 if (X!=Input) memcpy(X, Input, sizeof(cdouble)*N);
xue@1 169 CFFTC(Order, W, X, bitinv);
xue@1 170
xue@1 171 if (Amp) for (i=0; i<N; i++) Amp[i]=sqrt(X[i].x*X[i].x+X[i].y*X[i].y);
xue@1 172 if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(X[i].y, X[i].x);
xue@1 173 }//CFFTC
xue@1 174
xue@1 175 //---------------------------------------------------------------------------
Chris@5 176 /**
xue@1 177 function CIFFTCbii: basic complex IFFT module, applied after bit-inversed ordering of inputs
xue@1 178
xue@1 179 In: Order: integer, equals log2(Wid)
xue@1 180 W[Wid/2]: twiddle factors
xue@1 181 X[Wid]: complex spectrum
xue@1 182 Out: X[Wid]: complex waveform
xue@1 183
xue@1 184 No return value.
xue@1 185 */
xue@1 186 void CIFFTCbii(int Order, cdouble* W, cdouble* X)
xue@1 187 {
xue@1 188 int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2;
xue@1 189 cdouble Temp;
xue@1 190
xue@1 191 for (i=0; i<Order; i++)
xue@1 192 {
xue@1 193 ElemsPerGroup=1<<i;
xue@1 194 Groups=1<<(Order-i-1);
xue@1 195 X0=0;
xue@1 196 for (j=0; j<Groups; j++)
xue@1 197 {
xue@1 198 for (k=0; k<ElemsPerGroup; k++)
xue@1 199 {
xue@1 200 int kGroups=k*Groups;
xue@1 201 X1=X0+k;
xue@1 202 X2=X1+ElemsPerGroup;
xue@1 203 Temp.x=X[X2].x*W[kGroups].x+X[X2].y*W[kGroups].y,
xue@1 204 X[X2].y=-X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1 205 X[X2].x=Temp.x;
xue@1 206 Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1 207 X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1 208 X[X1]=Temp;
xue@1 209 }
xue@1 210 X0=X0+ElemsPerGroup*2;
xue@1 211 }
xue@1 212 }
xue@1 213 for (i=0; i<N; i++)
xue@1 214 {
xue@1 215 X[i].x/=N;
xue@1 216 X[i].y/=N;
xue@1 217 }
xue@1 218 }//CIFFTCbii
xue@1 219
Chris@5 220 /**
xue@1 221 function CIFFTC: in-place complex IFFT
xue@1 222
xue@1 223 In: Order: integer, equals log2(Wid)
xue@1 224 W[Wid/2]: twiddle factors
xue@1 225 X[Wid]: complex spectrum
xue@1 226 bitinv[Wid]: bit-inversion table
xue@1 227 Out: X[Wid]: complex waveform
xue@1 228
xue@1 229 No return value.
xue@1 230 */
xue@1 231 void CIFFTC(int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 232 {
xue@1 233 int N=1<<Order, i, jj, *bitinv1=bitinv;
xue@1 234 cdouble Temp;
xue@1 235 if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1 236 for (i=1; i<N-1; i++)
xue@1 237 {
xue@1 238 jj=bitinv[i];
xue@1 239 if (i<jj)
xue@1 240 {
xue@1 241 Temp=X[i];
xue@1 242 X[i]=X[jj];
xue@1 243 X[jj]=Temp;
xue@1 244 }
xue@1 245 }
xue@1 246 if (!bitinv1) free(bitinv);
xue@1 247 CIFFTCbii(Order, W, X);
xue@1 248 }//CIFFTC
xue@1 249
Chris@5 250 /**
xue@1 251 function CIFFTC: complex IFFT
xue@1 252
xue@1 253 In: Input[Wid]: complex spectrum
xue@1 254 Order: integer, equals log2(Wid)
xue@1 255 W[Wid/2]: twiddle factors
xue@1 256 bitinv[Wid]: bit-inversion table
xue@1 257 Out:X[Wid]: complex waveform
xue@1 258
xue@1 259 No return value.
xue@1 260 */
xue@1 261 void CIFFTC(cdouble* Input, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 262 {
xue@1 263 int N=1<<Order;
xue@1 264 if (X!=Input) memcpy(X, Input, sizeof(cdouble)*N);
xue@1 265 if (bitinv) CIFFTC(Order, W, X, bitinv);
xue@1 266 else CIFFTC(Order, W, X);
xue@1 267 }//CIFFTC
xue@1 268
xue@1 269 //---------------------------------------------------------------------------
Chris@5 270 /**
xue@1 271 function CIFFTR: complex-to-real IFFT
xue@1 272
xue@1 273 In: Input[Wid/2+1]: complex spectrum, imaginary parts of Input[0] and Input[Wid/2] are ignored
xue@1 274 Order: integer, equals log2(Wid)
xue@1 275 W[Wid/2]: twiddle factors
xue@1 276 hbitinv[Wid/2]: half bit-inversion table
xue@1 277 Out:X[Wid]: real waveform
xue@1 278
xue@1 279 No return value.
xue@1 280 */
xue@1 281 void CIFFTR(cdouble* Input, int Order, cdouble* W, double* X, int* hbitinv)
xue@1 282 {
xue@1 283 int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2, *hbitinv1=hbitinv;
xue@1 284 cdouble Temp;
xue@1 285
xue@1 286 Order--; N/=2;
xue@1 287 if (!hbitinv) hbitinv=CreateBitInvTable(Order);
xue@1 288
xue@1 289 cdouble* Xc=(cdouble*)X;
xue@1 290
xue@1 291 Xc[0].y=0.5*(Input[0].x-Input[N].x);
xue@1 292 Xc[0].x=0.5*(Input[0].x+Input[N].x);
xue@1 293 for (int i=1; i<N/2; i++)
xue@1 294 {
xue@1 295 double frp=Input[i].x+Input[N-i].x, frn=Input[i].x-Input[N-i].x,
xue@1 296 fip=Input[i].y+Input[N-i].y, fin=Input[i].y-Input[N-i].y;
xue@1 297 Xc[i].x=0.5*(frp+frn*W[i].y-fip*W[i].x);
xue@1 298 Xc[i].y=0.5*(fin+frn*W[i].x+fip*W[i].y);
xue@1 299 Xc[N-i].x=0.5*(frp-frn*W[i].y+fip*W[i].x);
xue@1 300 Xc[N-i].y=0.5*(-fin+frn*W[i].x+fip*W[i].y);
xue@1 301 }
xue@1 302 Xc[N/2].x=Input[N/2].x;
xue@1 303 Xc[N/2].y=-Input[N/2].y;
xue@1 304
xue@1 305 ElemsPerGroup=1<<Order;
xue@1 306 Groups=1;
xue@1 307
xue@1 308 for (i=0; i<Order; i++)
xue@1 309 {
xue@1 310 ElemsPerGroup/=2;
xue@1 311 X0=0;
xue@1 312 for (j=0; j<Groups; j++)
xue@1 313 {
xue@1 314 int kGroups=hbitinv[j];
xue@1 315 for (k=0; k<ElemsPerGroup; k++)
xue@1 316 {
xue@1 317 X1=X0+k;
xue@1 318 X2=X1+ElemsPerGroup;
xue@1 319 Temp.x=Xc[X2].x*W[kGroups].x+Xc[X2].y*W[kGroups].y,
xue@1 320 Xc[X2].y=-Xc[X2].x*W[kGroups].y+Xc[X2].y*W[kGroups].x;
xue@1 321 Xc[X2].x=Temp.x;
xue@1 322 Temp.x=Xc[X1].x+Xc[X2].x, Temp.y=Xc[X1].y+Xc[X2].y;
xue@1 323 Xc[X2].x=Xc[X1].x-Xc[X2].x, Xc[X2].y=Xc[X1].y-Xc[X2].y;
xue@1 324 Xc[X1].x=Temp.x, Xc[X1].y=Temp.y;
xue@1 325 }
xue@1 326 X0=X0+(ElemsPerGroup<<1);
xue@1 327 }
xue@1 328 Groups*=2;
xue@1 329 }
xue@1 330
xue@1 331 for (i=0; i<N; i++)
xue@1 332 {
xue@1 333 int jj=hbitinv[i];
xue@1 334 if (i<jj)
xue@1 335 {
xue@1 336 Temp=Xc[i];
xue@1 337 Xc[i]=Xc[jj];
xue@1 338 Xc[jj]=Temp;
xue@1 339 }
xue@1 340 }
xue@1 341 double norm=1.0/N;;
xue@1 342 N*=2;
xue@1 343 for (int i=0; i<N; i++) X[i]*=norm;
xue@1 344 if (!hbitinv1) free(hbitinv);
xue@1 345 }//CIFFTR
xue@1 346
xue@1 347 //---------------------------------------------------------------------------
Chris@5 348 /**
xue@1 349 function CreateBitInvTable: creates a table of bit-inversed integers
xue@1 350
xue@1 351 In: Order: interger, equals log2(size of table)
xue@1 352
xue@1 353 Returns a pointer to a newly allocated array containing the table. The returned pointer must be freed
xue@1 354 using free(), NOT delete[].
xue@1 355 */
xue@1 356 int* CreateBitInvTable(int Order)
xue@1 357 {
xue@1 358 int N=1<<Order;
xue@1 359 int* result=(int*)malloc(sizeof(int)*N);
xue@1 360 for (int i=0; i<N; i++) result[i]=BitInv(i, Order);
xue@1 361 return result;
xue@1 362 }//CreateBitInvTable
xue@1 363
xue@1 364
xue@1 365 //---------------------------------------------------------------------------
Chris@5 366 /**
xue@1 367 function RFFTC_ual: unaligned real-to-complex FFT
xue@1 368
xue@1 369 In: Input[Wid]: real waveform
xue@1 370 Order; integer, equals log2(Wid)
xue@1 371 W[Wid/2]: twiddle factors
xue@1 372 hbitinv[Wid/2]: half bit-inversion table
xue@1 373 Out:X[Wid}: complex spectrum
xue@1 374 Amp[Wid]: amplitude spectrum
xue@1 375 Arg[Wid]: phase spetrum
xue@1 376
xue@1 377 No return value.
xue@1 378 */
xue@1 379 void RFFTC_ual(double* Input, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1 380 {
xue@1 381 int N=1<<Order, i, j, k, *hbitinv1=hbitinv, Groups, ElemsPerGroup, X0, X1, X2;
xue@1 382 cdouble Temp, zp, zn;
xue@1 383
xue@1 384 N/=2; Order--;
xue@1 385
xue@1 386 //Input being NULL implies external initialization of X. This is used by RFFTCW but is not
xue@1 387 //recommended for external use.
xue@1 388 if (Input)
xue@1 389 {
xue@1 390 if (!hbitinv) hbitinv=CreateBitInvTable(Order);
xue@1 391
xue@1 392 if (Input==(double*)X)
xue@1 393 {
xue@1 394 //Input being identical to X is not recommended for external use.
xue@1 395 for (int i=0; i<N; i++)
xue@1 396 {
xue@1 397 int bi=hbitinv[i];
xue@1 398 if (i<bi) {cdouble tmp=X[i]; X[i]=X[bi]; X[bi]=tmp;}
xue@1 399 }
xue@1 400 }
xue@1 401 else
xue@1 402 {
xue@1 403 for (i=0; i<N; i++) X[i]=((cdouble*)Input)[hbitinv[i]];
xue@1 404 }
xue@1 405 if (!hbitinv1) free(hbitinv);
xue@1 406 }
xue@1 407
xue@1 408 for (i=0; i<Order; i++)
xue@1 409 {
xue@1 410 ElemsPerGroup=1<<i;
xue@1 411 Groups=1<<(Order-i-1);
xue@1 412 X0=0;
xue@1 413 for (j=0; j<Groups; j++)
xue@1 414 {
xue@1 415 for (k=0; k<ElemsPerGroup; k++)
xue@1 416 {
xue@1 417 X1=X0+k;
xue@1 418 X2=X1+ElemsPerGroup;
xue@1 419 int kGroups=k*2*Groups;
xue@1 420 Temp.x=X[X2].x*W[kGroups].x-X[X2].y*W[kGroups].y,
xue@1 421 X[X2].y=X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1 422 X[X2].x=Temp.x;
xue@1 423 Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1 424 X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1 425 X[X1]=Temp;
xue@1 426 }
xue@1 427 X0=X0+(ElemsPerGroup<<1);
xue@1 428 }
xue@1 429 }
xue@1 430 zp.x=X[0].x+X[0].y, zn.x=X[0].x-X[0].y;
xue@1 431 X[0].x=zp.x;
xue@1 432 X[0].y=0;
xue@1 433 X[N].x=zn.x;
xue@1 434 X[N].y=0;
xue@1 435 for (int k=1; k<N/2; k++)
xue@1 436 {
xue@1 437 zp.x=X[k].x+X[N-k].x, zn.x=X[k].x-X[N-k].x,
xue@1 438 zp.y=X[k].y+X[N-k].y, zn.y=X[k].y-X[N-k].y;
xue@1 439 X[k].x=0.5*(zp.x+W[k].y*zn.x+W[k].x*zp.y);
xue@1 440 X[k].y=0.5*(zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1 441 X[N-k].x=0.5*(zp.x-W[k].y*zn.x-W[k].x*zp.y);
xue@1 442 X[N-k].y=0.5*(-zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1 443 }
xue@1 444 //X[N/2].x=X[N/2].x;
xue@1 445 X[N/2].y=-X[N/2].y;
xue@1 446 N*=2;
xue@1 447
xue@1 448 for (int k=N/2+1; k<N; k++) X[k].x=X[N-k].x, X[k].y=-X[N-k].y;
xue@1 449 if (Amp) for (i=0; i<N; i++) Amp[i]=sqrt(X[i].x*X[i].x+X[i].y*X[i].y);
xue@1 450 if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(X[i].x, X[i].y);
xue@1 451 }//RFFTC_ual
xue@1 452
xue@1 453 //---------------------------------------------------------------------------
Chris@5 454 /**
xue@1 455 function RFFTCW: real-to-complex FFT with window
xue@1 456
xue@1 457 In: Input[Wid]: real waveform
xue@1 458 Window[Wid]: window function
xue@1 459 Order; integer, equals log2(Wid)
xue@1 460 W[Wid/2]: twiddle factors
xue@1 461 hbitinv[Wid/2]: half bit-inversion table
xue@1 462 Out:X[Wid}: complex spectrum
xue@1 463 Amp[Wid]: amplitude spectrum
xue@1 464 Arg[Wid]: phase spetrum
xue@1 465
xue@1 466 No return value.
xue@1 467 */
xue@1 468 void RFFTCW(double* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1 469 {
xue@1 470 int N=1<<Order, *hbitinv1=hbitinv;
xue@1 471 if (!hbitinv) hbitinv=CreateBitInvTable(Order-1);
xue@1 472 N/=2;
xue@1 473
xue@1 474 if (Input==(double*)X)
xue@1 475 { //so that X[n].x IS Input[2n], X[n].y IS Input[2n+1]
xue@1 476 for (int n=0; n<N; n++)
xue@1 477 {
xue@1 478 int bi=hbitinv[n], n2=n+n, bi2=bi+bi;
xue@1 479 if (n<bi)
xue@1 480 {
xue@1 481 double tmp=X[n].x*Window[n2]; X[n].x=X[bi].x*Window[bi2]; X[bi].x=tmp;
xue@1 482 tmp=X[n].y*Window[n2+1]; X[n].y=X[bi].y*Window[bi2+1]; X[bi].y=tmp;
xue@1 483 }
xue@1 484 else if (n==bi)
xue@1 485 { //so that X[n].x IS Input[bi]
xue@1 486 X[n].x*=Window[bi2], X[n].y*=Window[bi2+1];
xue@1 487 }
xue@1 488 }
xue@1 489 }
xue@1 490 else
xue@1 491 {
xue@1 492 for (int n=0; n<N; n++)
xue@1 493 {
xue@1 494 int bi=hbitinv[n], bi2=bi+bi; X[n].x=Input[bi2]*Window[bi2], X[n].y=Input[bi2+1]*Window[bi2+1];
xue@1 495 }
xue@1 496 }
xue@1 497
xue@1 498 RFFTC_ual(0, Amp, Arg, Order, W, X, hbitinv);
xue@1 499 if (!hbitinv1) free(hbitinv);
xue@1 500 }//RFFTCW
xue@1 501
Chris@5 502 /**
xue@1 503 function RFFTCW: real-to-complex FFT with window
xue@1 504
xue@1 505 In: Input[Wid]: real waveform as _int16 array
xue@1 506 Window[Wid]: window function
xue@1 507 Order; integer, equals log2(Wid)
xue@1 508 W[Wid/2]: twiddle factors
xue@1 509 hbitinv[Wid/2]: half bit-inversion table
xue@1 510 Out:X[Wid}: complex spectrum
xue@1 511 Amp[Wid]: amplitude spectrum
xue@1 512 Arg[Wid]: phase spetrum
xue@1 513
xue@1 514 No return value.
xue@1 515 */
xue@1 516 void RFFTCW(__int16* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1 517 {
xue@1 518 int N=1<<Order, *hbitinv1=hbitinv;
xue@1 519
xue@1 520 N/=2;
xue@1 521 if (!hbitinv) hbitinv=CreateBitInvTable(Order-1);
xue@1 522 for (int n=0; n<N; n++)
xue@1 523 {
xue@1 524 int bi=hbitinv[n], bi2=bi+bi; X[n].x=Input[bi2]*Window[bi2], X[n].y=Input[bi2+1]*Window[bi2+1];
xue@1 525 }
xue@1 526
xue@1 527 RFFTC_ual(0, Amp, Arg, Order, W, X, hbitinv);
xue@1 528 if (!hbitinv1) free(hbitinv);
xue@1 529 }//RFFTCW
xue@1 530
xue@1 531 //---------------------------------------------------------------------------
Chris@5 532 /**
xue@1 533 function CFFTCW: complex FFT with window
xue@1 534
xue@1 535 In: Window[Wid]: window function
xue@1 536 Order: integer, equals log2(Wid)
xue@1 537 W[Wid/2]: twiddle factors
xue@1 538 X[Wid]: complex waveform
xue@1 539 bitinv[Wid]: bit-inversion table
xue@1 540 Out:X[Wid], complex spectrum
xue@1 541
xue@1 542 No return value.
xue@1 543 */
xue@1 544 void CFFTCW(double* Window, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 545 {
xue@1 546 int N=1<<Order;
xue@1 547 for (int i=0; i<N; i++) X[i].x*=Window[i], X[i].y*=Window[i];
xue@1 548 CFFTC(Order, W, X, bitinv);
xue@1 549 }//CFFTCW
xue@1 550
Chris@5 551 /**
xue@1 552 function CFFTCW: complex FFT with window
xue@1 553
xue@1 554 In: Input[Wid]: complex waveform
xue@1 555 Window[Wid]: window function
xue@1 556 Order: integer, equals log2(Wid)
xue@1 557 W[Wid/2]: twiddle factors
xue@1 558 X[Wid]: complex waveform
xue@1 559 bitinv[Wid]: bit-inversion table
xue@1 560 Out:X[Wid], complex spectrum
xue@1 561 Amp[Wid], amplitude spectrum
xue@1 562 Arg[Wid], phase spectrum
xue@1 563
xue@1 564 No return value.
xue@1 565 */
xue@1 566 void CFFTCW(cdouble* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1 567 {
xue@1 568 int N=1<<Order;
xue@1 569 for (int i=0; i<N; i++) X[i].x=Input[i].x*Window[i], X[i].y=Input[i].y*Window[i];
xue@1 570 CFFTC(X, Amp, Arg, Order, W, X, bitinv);
xue@1 571 }//CFFTCW
xue@1 572
xue@1 573 //---------------------------------------------------------------------------
Chris@5 574 /**
xue@1 575 function RDCT1: fast DCT1 implemented using FFT. DCT-I has the time scale 0.5-sample shifted from the DFT.
xue@1 576
xue@1 577 In: Input[Wid]: real waveform
xue@1 578 Order: integer, equals log2(Wid)
xue@1 579 qW[Wid/8]: quarter table of twiddle factors
xue@1 580 qX[Wid/4]: quarter FFT data buffer
xue@1 581 qbitinv[Wid/4]: quarter bit-inversion table
xue@1 582 Out:Output[Wid]: DCT-I of Input.
xue@1 583
xue@1 584 No return value. Content of qW is destroyed on return. On calling the fft buffers should be allocated
xue@1 585 to size 0.25*Wid.
xue@1 586 */
xue@1 587 void RDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv)
xue@1 588 {
xue@1 589 const double lmd0=sqrt(0.5);
xue@1 590 if (Order==0)
xue@1 591 {
xue@1 592 Output[0]=Input[0]*lmd0;
xue@1 593 return;
xue@1 594 }
xue@1 595 if (Order==1)
xue@1 596 {
xue@1 597 double tmp1=(Input[0]+Input[1])*lmd0;
xue@1 598 Output[1]=(Input[0]-Input[1])*lmd0;
xue@1 599 Output[0]=tmp1;
xue@1 600 return;
xue@1 601 }
xue@1 602 int order=Order-1, N=1<<(Order-1), C=1;
xue@1 603 bool createbitinv=!qbitinv;
xue@1 604 if (createbitinv) qbitinv=CreateBitInvTable(Order-2);
xue@1 605 double *even=(double*)malloc8(sizeof(double)*N*2);
xue@1 606 double *odd=&even[N];
xue@1 607 //data pass from Input to Output, combined with the first sequence split
xue@1 608 for (int i=0, N2=N*2; i<N; i++)
xue@1 609 {
xue@1 610 even[i]=Input[i]+Input[N2-1-i];
xue@1 611 odd[i]=Input[i]-Input[N2-1-i];
xue@1 612 }
xue@1 613 for (int i=0; i<N; i++) Output[i*2]=even[i], Output[i*2+1]=odd[i];
xue@1 614 while (order>1)
xue@1 615 {
xue@1 616 //N=2^order, 4|N, 2|hN
xue@1 617 //keep subsequence 0, 2C, 4C, ... 2(N-1)C
xue@1 618 //process sequence C, 3C, ..., (2N-1)C only
xue@1 619 //RDCT4 routine
xue@1 620 int hN=N/2, N2=N*2;
xue@1 621 for (int i=0; i<hN; i++)
xue@1 622 {
xue@1 623 double b=Output[(2*(2*i)+1)*C], c=Output[(2*(N-1-2*i)+1)*C], theta=-(i+0.25)*M_PI/N;
xue@1 624 double co=cos(theta), si=sin(theta);
xue@1 625 qX[i].x=b*co-c*si, qX[i].y=c*co+b*si;
xue@1 626 }
xue@1 627 CFFTC(order-1, qW, qX, qbitinv);
xue@1 628 for (int i=0; i<hN; i++)
xue@1 629 {
xue@1 630 double gr=qX[i].x, gi=qX[i].y, theta=-i*M_PI/N;
xue@1 631 double co=cos(theta), si=sin(theta);
xue@1 632 Output[(4*i+1)*C]=gr*co-gi*si;
xue@1 633 Output[(N2-4*i-1)*C]=-gr*si-gi*co;
xue@1 634 }
xue@1 635 N=(N>>1);
xue@1 636 C=(C<<1);
xue@1 637 for (int i=1; i<N/4; i++)
xue@1 638 qW[i]=qW[i*2];
xue@1 639 for (int i=1; i<N/2; i++)
xue@1 640 qbitinv[i]=qbitinv[i*2];
xue@1 641
xue@1 642 //splitting subsequence 0, 2C, 4C, ..., 2(N-1)C
xue@1 643 for (int i=0, N2=N*2; i<N; i++)
xue@1 644 {
xue@1 645 even[i]=Output[i*C]+Output[(N2-1-i)*C];
xue@1 646 odd[i]=Output[i*C]-Output[(N2-1-i)*C];
xue@1 647 }
xue@1 648 for (int i=0; i<N; i++)
xue@1 649 {
xue@1 650 Output[2*i*C]=even[i];
xue@1 651 Output[(2*i+1)*C]=odd[i];
xue@1 652 }
xue@1 653 order--;
xue@1 654 }
xue@1 655 //order==1
xue@1 656 //use C and 3C in DCT4
xue@1 657 //use 0 and 2C in DCT1
xue@1 658 double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1 659 double tmp1=c1*Output[C]+c2*Output[3*C];
xue@1 660 Output[3*C]=c2*Output[C]-c1*Output[3*C];
xue@1 661 Output[C]=tmp1;
xue@1 662 tmp1=Output[0]+Output[2*C];
xue@1 663 Output[2*C]=(Output[0]-Output[2*C])*lmd0;
xue@1 664 Output[0]=tmp1*lmd0;
xue@1 665
xue@1 666 if (createbitinv) free(qbitinv);
xue@1 667 free8(even);
xue@1 668 }//RDCT1
xue@1 669
xue@1 670 //---------------------------------------------------------------------------
Chris@5 671 /**
xue@1 672 function RDCT4: fast DCT4 implemented using FFT. DCT-IV has both the time and frequency scales
xue@1 673 0.5-sample or 0.5-bin shifted from DFT.
xue@1 674
xue@1 675 In: Input[Wid]: real waveform
xue@1 676 Order: integer, equals log2(Wid)
xue@1 677 hW[Wid/4]: half table of twiddle factors
xue@1 678 hX[Wid/2]: hal FFT data buffer
xue@1 679 hbitinv[Wid/2]: half bit-inversion table
xue@1 680 Out:Output[Wid]: DCT-IV of Input.
xue@1 681
xue@1 682 No return value. Content of hW not affected. On calling the fft buffers should be allocated to size
xue@1 683 0.5*Wid.
xue@1 684 */
xue@1 685 void RDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv)
xue@1 686 {
xue@1 687 if (Order==0)
xue@1 688 {
xue@1 689 Output[0]=Input[0]/sqrt(2);
xue@1 690 return;
xue@1 691 }
xue@1 692 if (Order==1)
xue@1 693 {
xue@1 694 double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1 695 double tmp1=c1*Input[0]+c2*Input[1];
xue@1 696 Output[1]=c2*Input[0]-c1*Input[1];
xue@1 697 Output[0]=tmp1;
xue@1 698 return;
xue@1 699 }
xue@1 700 int N=1<<Order, hN=N/2;
xue@1 701 for (int i=0; i<hN; i++)
xue@1 702 {
xue@1 703 double b=Input[2*i], c=Input[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1 704 double co=cos(theta), si=sin(theta);
xue@1 705 hX[i].x=b*co-c*si, hX[i].y=c*co+b*si;
xue@1 706 }
xue@1 707 CFFTC(Order-1, hW, hX, hbitinv);
xue@1 708 for (int i=0; i<hN; i++)
xue@1 709 {
xue@1 710 double gr=hX[i].x, gi=hX[i].y, theta=-i*M_PI/N;
xue@1 711 double co=cos(theta), si=sin(theta);
xue@1 712 Output[2*i]=gr*co-gi*si;
xue@1 713 Output[N-1-2*i]=-gr*si-gi*co;
xue@1 714 }
xue@1 715 }//RDCT4
xue@1 716
xue@1 717 //---------------------------------------------------------------------------
Chris@5 718 /**
xue@1 719 function RIDCT1: fast IDCT1 implemented using FFT.
xue@1 720
xue@1 721 In: Input[Wid]: DCT-I of some real waveform.
xue@1 722 Order: integer, equals log2(Wid)
xue@1 723 qW[Wid/8]: quarter table of twiddle factors
xue@1 724 qX[Wid/4]: quarter FFT data buffer
xue@1 725 qbitinv[Wid/4]: quarter bit-inversion table
xue@1 726 Out:Output[Wid]: IDCT-I of Input.
xue@1 727
xue@1 728 No return value. Content of qW is destroyed on return. On calling the fft buffers should be allocated
xue@1 729 to size 0.25*Wid.
xue@1 730 */
xue@1 731 void RIDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv)
xue@1 732 {
xue@1 733 const double lmd0=sqrt(0.5);
xue@1 734 if (Order==0)
xue@1 735 {
xue@1 736 Output[0]=Input[0]/lmd0;
xue@1 737 return;
xue@1 738 }
xue@1 739 if (Order==1)
xue@1 740 {
xue@1 741 double tmp1=(Input[0]+Input[1])*lmd0;
xue@1 742 Output[1]=(Input[0]-Input[1])*lmd0;
xue@1 743 Output[0]=tmp1;
xue@1 744 return;
xue@1 745 }
xue@1 746 int order=Order-1, N=1<<(Order-1), C=1;
xue@1 747 bool createbitinv=!qbitinv;
xue@1 748 if (createbitinv) qbitinv=CreateBitInvTable(Order-2);
xue@1 749 double *even=(double*)malloc8(sizeof(double)*N*2);
xue@1 750 double *odd=&even[N];
xue@1 751
xue@1 752 while (order>0)
xue@1 753 {
xue@1 754 //N=2^order, 4|N, 2|hN
xue@1 755 //keep subsequence 0, 2C, 4C, ... 2(N-1)C
xue@1 756 //process sequence C, 3C, ..., (2N-1)C only
xue@1 757 //data pass from Input
xue@1 758 for (int i=0; i<N; i++)
xue@1 759 {
xue@1 760 odd[i]=Input[(i*2+1)*C];
xue@1 761 }
xue@1 762
xue@1 763 //IDCT4 routine
xue@1 764 //RIDCT4(odd, odd, order, qW, qX, qbitinv);
xue@1 765
xue@1 766 if (order==1)
xue@1 767 {
xue@1 768 double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1 769 double tmp1=c1*odd[0]+c2*odd[1];
xue@1 770 odd[1]=c2*odd[0]-c1*odd[1];
xue@1 771 odd[0]=tmp1;
xue@1 772 }
xue@1 773 else
xue@1 774 {
xue@1 775 int hN=N/2;
xue@1 776 for (int i=0; i<hN; i++)
xue@1 777 {
xue@1 778 double b=odd[2*i], c=odd[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1 779 double co=cos(theta), si=sin(theta);
xue@1 780 qX[i].x=b*co-c*si, qX[i].y=c*co+b*si;
xue@1 781 }
xue@1 782 CFFTC(order-1, qW, qX, qbitinv);
xue@1 783 double i2N=2.0/N;
xue@1 784 for (int i=0; i<hN; i++)
xue@1 785 {
xue@1 786 double gr=qX[i].x, gi=qX[i].y, theta=-i*M_PI/N;
xue@1 787 double co=cos(theta), si=sin(theta);
xue@1 788 odd[2*i]=(gr*co-gi*si)*i2N;
xue@1 789 odd[N-1-2*i]=(-gr*si-gi*co)*i2N;
xue@1 790 }
xue@1 791 }
xue@1 792
xue@1 793 order--;
xue@1 794 N=(N>>1);
xue@1 795 C=(C<<1);
xue@1 796 for (int i=1; i<N/4; i++)
xue@1 797 qW[i]=qW[i*2];
xue@1 798 for (int i=1; i<N/2; i++)
xue@1 799 qbitinv[i]=qbitinv[i*2];
xue@1 800
xue@1 801 odd=&even[N];
xue@1 802 }
xue@1 803 //order==0
xue@1 804 even[0]=Input[0];
xue@1 805 even[1]=Input[C];
xue@1 806 double tmp1=(even[0]+even[1])*lmd0;
xue@1 807 Output[1]=(even[0]-even[1])*lmd0;
xue@1 808 Output[0]=tmp1;
xue@1 809 order++;
xue@1 810
xue@1 811 while (order<Order)
xue@1 812 {
xue@1 813 N=(N<<1);
xue@1 814 odd=&even[N];
xue@1 815 for (int i=0; i<N; i++)
xue@1 816 {
xue@1 817 Output[N*2-1-i]=(Output[i]-odd[i])/2;
xue@1 818 Output[i]=(Output[i]+odd[i])/2;
xue@1 819 }
xue@1 820 order++;
xue@1 821 }
xue@1 822
xue@1 823 if (createbitinv) free(qbitinv);
xue@1 824 free8(even);
xue@1 825 }//RIDCT1
xue@1 826
xue@1 827 //---------------------------------------------------------------------------
Chris@5 828 /**
xue@1 829 function RIDCT4: fast IDCT4 implemented using FFT.
xue@1 830
xue@1 831 In: Input[Wid]: DCT-IV of some real waveform
xue@1 832 Order: integer, equals log2(Wid)
xue@1 833 hW[Wid/4]: half table of twiddle factors
xue@1 834 hX[Wid/2]: hal FFT data buffer
xue@1 835 hbitinv[Wid/2]: half bit-inversion table
xue@1 836 Out:Output[Wid]: IDCT-IV of Input.
xue@1 837
xue@1 838 No return value. Content of hW not affected. On calling the fft buffers should be allocated to size
xue@1 839 0.5*Wid.
xue@1 840 */
xue@1 841 void RIDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv)
xue@1 842 {
xue@1 843 if (Order==0)
xue@1 844 {
xue@1 845 Output[0]=Input[0]*sqrt(2);
xue@1 846 return;
xue@1 847 }
xue@1 848 if (Order==1)
xue@1 849 {
xue@1 850 double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1 851 double tmp1=c1*Input[0]+c2*Input[1];
xue@1 852 Output[1]=c2*Input[0]-c1*Input[1];
xue@1 853 Output[0]=tmp1;
xue@1 854 return;
xue@1 855 }
xue@1 856 int N=1<<Order, hN=N/2;
xue@1 857 for (int i=0; i<hN; i++)
xue@1 858 {
xue@1 859 double b=Input[2*i], c=Input[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1 860 double co=cos(theta), si=sin(theta);
xue@1 861 hX[i].x=b*co-c*si, hX[i].y=c*co+b*si;
xue@1 862 }
xue@1 863 CFFTC(Order-1, hW, hX, hbitinv);
xue@1 864 double i2N=2.0/N;
xue@1 865 for (int i=0; i<hN; i++)
xue@1 866 {
xue@1 867 double gr=hX[i].x, gi=hX[i].y, theta=-i*M_PI/N;
xue@1 868 double co=cos(theta), si=sin(theta);
xue@1 869 Output[2*i]=(gr*co-gi*si)*i2N;
xue@1 870 Output[N-1-2*i]=(-gr*si-gi*co)*i2N;
xue@1 871 }
xue@1 872 }//RIDCT4
xue@1 873
xue@1 874 //---------------------------------------------------------------------------
Chris@5 875 /**
xue@1 876 function RLCT: real local cosine transform of equal frame widths Wid=2^Order
xue@1 877
xue@1 878 In: data[Count]: real waveform
xue@1 879 Order: integer, equals log2(Wid), Wid being the cosine atom size
xue@1 880 g[wid]: lap window, designed so that g[k] increases from 0 to 1 and g[k]^2+g[wid-1-k]^2=1
xue@1 881 example: wid=4, g[k]=sin(pi*(k+0.5)/8).
xue@1 882 Out:spec[Fr][Wid]: the local cosine transform
xue@1 883
xue@1 884 No return value.
xue@1 885 */
xue@1 886 void RLCT(double** spec, double* data, int Count, int Order, int wid, double* g)
xue@1 887 {
xue@1 888 int Wid=1<<Order, Fr=Count/Wid, hwid=wid/2;
xue@1 889 int* hbitinv=CreateBitInvTable(Order-1);
xue@1 890 cdouble *hx=(cdouble*)malloc8(sizeof(cdouble)*Wid*3/4), *hw=(cdouble*)&hx[Wid/2];
xue@1 891 double norm=sqrt(2.0/Wid);
xue@1 892 SetTwiddleFactors(Wid/2, hw);
xue@1 893
xue@1 894 for (int fr=0; fr<Fr; fr++)
xue@1 895 {
xue@1 896 if (fr==0)
xue@1 897 {
xue@1 898 memcpy(spec[fr], data, sizeof(double)*(Wid-hwid));
xue@1 899 for (int i=0, k=Wid+hwid-1, l=Wid-hwid; i<hwid; i++, k--, l++)
xue@1 900 spec[fr][l]=data[l]*g[wid-1-i]-data[k]*g[i];
xue@1 901 }
xue@1 902 else if (fr==Fr-1)
xue@1 903 {
xue@1 904 for (int i=hwid, j=fr*Wid, k=fr*Wid-1, l=0; i<wid; i++, j++, k--, l++)
xue@1 905 spec[fr][l]=data[k]*g[wid-1-i]+data[j]*g[i];
xue@1 906 memcpy(&spec[fr][hwid], &data[fr*Wid+hwid], sizeof(double)*(Wid-hwid));
xue@1 907 }
xue@1 908 else
xue@1 909 {
xue@1 910 for (int i=hwid, j=fr*Wid, k=fr*Wid-1, l=0; i<wid; i++, j++, k--, l++)
xue@1 911 spec[fr][l]=data[k]*g[wid-1-i]+data[j]*g[i];
xue@1 912 if (Wid>wid) memcpy(&spec[fr][hwid], &data[fr*Wid+hwid], sizeof(double)*(Wid-wid));
xue@1 913 for (int i=0, j=(fr+1)*Wid-hwid, k=(fr+1)*Wid+hwid-1, l=Wid-hwid; i<hwid; i++, j++, k--, l++)
xue@1 914 spec[fr][l]=data[j]*g[wid-1-i]-data[k]*g[i];
xue@1 915 }
xue@1 916 }
xue@1 917 for (int fr=0; fr<Fr; fr++)
xue@1 918 {
xue@1 919 if (fr==Fr-1)
xue@1 920 {
xue@1 921 for (int i=1; i<Wid/4; i++) hw[i]=hw[2*i], hbitinv[i]=hbitinv[2*i];
xue@1 922 RDCT1(spec[fr], spec[fr], Order, hw, hx, hbitinv);
xue@1 923 }
xue@1 924 else
xue@1 925 RDCT4(spec[fr], spec[fr], Order, hw, hx, hbitinv);
xue@1 926
xue@1 927 ////The following line can be removed if the corresponding line in RILCT(...) is removed
xue@1 928 for (int i=0; i<Wid; i++) spec[fr][i]*=norm;
xue@1 929 }
xue@1 930 free(hbitinv);
xue@1 931 free8(hx);
xue@1 932 }//RLCT
xue@1 933
xue@1 934 //---------------------------------------------------------------------------
Chris@5 935 /**
xue@1 936 function RILCT: inverse local cosine transform of equal frame widths Wid=2^Order
xue@1 937
xue@1 938 In: spec[Fr][Wid]: the local cosine transform
xue@1 939 Order: integer, equals log2(Wid), Wid being the cosine atom size
xue@1 940 g[wid]: lap window, designed so that g[k] increases from 0 to 1 and g[k]^2+g[wid-1-k]^2=1.
xue@1 941 example: wid=4, g[k]=sin(pi*(k+0.5)/8).
xue@1 942 Out:data[Count]: real waveform
xue@1 943
xue@1 944 No return value.
xue@1 945 */
xue@1 946 void RILCT(double* data, double** spec, int Fr, int Order, int wid, double* g)
xue@1 947 {
xue@1 948 int Wid=1<<Order, Count=Fr*Wid, hwid=wid/2, *hbitinv=CreateBitInvTable(Order-1);
xue@1 949 cdouble *hx=(cdouble*)malloc8(sizeof(cdouble)*Wid*3/4), *hw=&hx[Wid/2];
xue@1 950 double norm=sqrt(Wid/2.0);
xue@1 951 SetTwiddleFactors(Wid/2, hw);
xue@1 952
xue@1 953 for (int fr=0; fr<Fr; fr++)
xue@1 954 {
xue@1 955 if (fr==Fr-1)
xue@1 956 {
xue@1 957 for (int i=1; i<Wid/4; i++) hw[i]=hw[2*i], hbitinv[i]=hbitinv[i*2];
xue@1 958 RIDCT1(spec[fr], &data[fr*Wid], Order, hw, hx, hbitinv);
xue@1 959 }
xue@1 960 else
xue@1 961 RIDCT4(spec[fr], &data[fr*Wid], Order, hw, hx, hbitinv);
xue@1 962 }
xue@1 963 //unfolding
xue@1 964 for (int fr=1; fr<Fr; fr++)
xue@1 965 {
xue@1 966 double* h=&data[fr*Wid];
xue@1 967 for (int i=0; i<hwid; i++)
xue@1 968 {
xue@1 969 double a=h[i], b=h[-1-i], c=g[i+hwid], d=g[hwid-1-i];
xue@1 970 h[i]=a*c-b*d, h[-i-1]=b*c+a*d;
xue@1 971 }
xue@1 972 }
xue@1 973
xue@1 974 free8(hx);
xue@1 975 ////The following line can be removed if the corresponding line in RLCT(...) is removed
xue@1 976 for (int i=0; i<Count; i++) data[i]*=norm;
xue@1 977 }//RILCT
xue@1 978
xue@1 979 //---------------------------------------------------------------------------
Chris@5 980 /**
xue@1 981 function CMFTC: radix-2 complex multiresolution Fourier transform
xue@1 982
xue@1 983 In: x[Wid]: complex waveform
xue@1 984 Order: integer, equals log2(Wid)
xue@1 985 W[Wid/2]: twiddle factors
xue@1 986 Out:X[Order+1][Wid]: multiresolution FT of x. X[0] is the same as x itself.
xue@1 987
xue@1 988 No return value.
xue@1 989
xue@1 990 Further reading: Wen X. and M. Sandler, "Calculation of radix-2 discrete multiresolution Fourier
xue@1 991 transform," Signal Processing, vol.87 no.10, 2007, pp.2455-2460.
xue@1 992 */
xue@1 993 void CMFTC(cdouble* x, int Order, cdouble** X, cdouble* W)
xue@1 994 {
xue@1 995 X[0]=x;
xue@1 996 for (int n=1, L=1<<(Order-1), M=2; n<=Order; n++, L>>=1, M<<=1)
xue@1 997 {
xue@1 998 cdouble *Xn=X[n], *Xp=X[n-1];
xue@1 999 for (int l=0; l<L; l++)
xue@1 1000 {
xue@1 1001 cdouble* lX=&Xn[l*M];
xue@1 1002 for (int m=0; m<M/2; m++)
xue@1 1003 {
xue@1 1004 lX[m].x=Xp[l*M+m].x+Xp[l*M+M/2+m].x;
xue@1 1005 lX[m].y=Xp[l*M+m].y+Xp[l*M+M/2+m].y;
xue@1 1006 double tmpr=x[l*M+m].x-x[l*M+M/2+m].x, tmpi=x[l*M+m].y-x[l*M+M/2+m].y;
xue@1 1007 int iw=m*L;
xue@1 1008 double wr=W[iw].x, wi=W[iw].y;
xue@1 1009 lX[M/2+m].x=tmpr*wr-tmpi*wi;
xue@1 1010 lX[M/2+m].y=tmpr*wi+tmpi*wr;
xue@1 1011 }
xue@1 1012 if (n==1) {}
xue@1 1013 else if (n==2) //two-point DFT
xue@1 1014 {
xue@1 1015 cdouble *aX=&X[n][l*M+M/2];
xue@1 1016 double tmp;
xue@1 1017 tmp=aX[0].x+aX[1].x; aX[1].x=aX[0].x-aX[1].x; aX[0].x=tmp;
xue@1 1018 tmp=aX[0].y+aX[1].y; aX[1].y=aX[0].y-aX[1].y; aX[0].y=tmp;
xue@1 1019 }
xue@1 1020 else if (n==3) //4-point DFT
xue@1 1021 {
xue@1 1022 cdouble *aX=&X[n][l*M+M/2];
xue@1 1023 double tmp, tmp2;
xue@1 1024 tmp=aX[0].x+aX[2].x; aX[2].x=aX[0].x-aX[2].x; aX[0].x=tmp;
xue@1 1025 tmp=aX[0].y+aX[2].y; aX[2].y=aX[0].y-aX[2].y; aX[0].y=tmp;
xue@1 1026 tmp=aX[1].y+aX[3].y; tmp2=aX[1].y-aX[3].y; aX[1].y=tmp;
xue@1 1027 tmp=aX[3].x-aX[1].x; aX[1].x+=aX[3].x; aX[3].x=tmp2; aX[3].y=tmp;
xue@1 1028 tmp=aX[0].x+aX[1].x; aX[1].x=aX[0].x-aX[1].x; aX[0].x=tmp;
xue@1 1029 tmp=aX[0].y+aX[1].y; aX[1].y=aX[0].y-aX[1].y; aX[0].y=tmp;
xue@1 1030 tmp=aX[2].x+aX[3].x; aX[3].x=aX[2].x-aX[3].x; aX[2].x=tmp;
xue@1 1031 tmp=aX[2].y+aX[3].y; aX[3].y=aX[2].y-aX[3].y; aX[2].y=tmp;
xue@1 1032 }
xue@1 1033 else //n>3
xue@1 1034 {
xue@1 1035 cdouble *aX=&X[n][l*M+M/2];
xue@1 1036 for (int an=1, aL=1, aM=M/2; an<n; aL*=2, aM/=2, an++)
xue@1 1037 {
xue@1 1038 for (int al=0; al<aL; al++)
xue@1 1039 for (int am=0; am<aM/2; am++)
xue@1 1040 {
xue@1 1041 int iw=am*2*aL*L;
xue@1 1042 cdouble *lX=&aX[al*aM];
xue@1 1043 double x1r=lX[am].x, x1i=lX[am].y,
xue@1 1044 x2r=lX[aM/2+am].x, x2i=lX[aM/2+am].y;
xue@1 1045 lX[am].x=x1r+x2r, lX[am].y=x1i+x2i;
xue@1 1046 x1r=x1r-x2r, x1i=x1i-x2i;
xue@1 1047 lX[aM/2+am].x=x1r*W[iw].x-x1i*W[iw].y,
xue@1 1048 lX[aM/2+am].y=x1r*W[iw].y+x1i*W[iw].x;
xue@1 1049 }
xue@1 1050 }
xue@1 1051 }
xue@1 1052 }
xue@1 1053 }
xue@1 1054 }//CMFTC
xue@1 1055
xue@1 1056
xue@1 1057 //---------------------------------------------------------------------------
xue@1 1058 /*
xue@1 1059 Old versions no longer in use. For reference only.
xue@1 1060 */
xue@1 1061 void RFFTC_ual_old(double* Input, double *Amp, double *Arg, int Order, cdouble* W, double* XR, double* XI, int* bitinv)
xue@1 1062 {
xue@1 1063 int N=1<<Order, i, j, jj, k, *bitinv1=bitinv, Groups, ElemsPerGroup, X0, X1, X2;
xue@1 1064 cdouble Temp, zp, zn;
xue@1 1065
xue@1 1066 if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1 1067 if (XR!=Input) for (i=0; i<N; i++) XR[i]=Input[bitinv[i]];
xue@1 1068 else for (i=0; i<N; i++) {jj=bitinv[i]; if (i<jj) {Temp.x=XR[i]; XR[i]=XR[jj]; XR[jj]=Temp.x;}}
xue@1 1069 N/=2;
xue@1 1070 double* XII=&XR[N];
xue@1 1071 Order--;
xue@1 1072 if (!bitinv1) free(bitinv);
xue@1 1073 for (i=0; i<Order; i++)
xue@1 1074 {
xue@1 1075 ElemsPerGroup=1<<i;
xue@1 1076 Groups=1<<(Order-i-1);
xue@1 1077 X0=0;
xue@1 1078 for (j=0; j<Groups; j++)
xue@1 1079 {
xue@1 1080 for (k=0; k<ElemsPerGroup; k++)
xue@1 1081 {
xue@1 1082 X1=X0+k;
xue@1 1083 X2=X1+ElemsPerGroup;
xue@1 1084 int kGroups=k*2*Groups;
xue@1 1085 Temp.x=XR[X2]*W[kGroups].x-XII[X2]*W[kGroups].y,
xue@1 1086 XII[X2]=XR[X2]*W[kGroups].y+XII[X2]*W[kGroups].x;
xue@1 1087 XR[X2]=Temp.x;
xue@1 1088 Temp.x=XR[X1]+XR[X2], Temp.y=XII[X1]+XII[X2];
xue@1 1089 XR[X2]=XR[X1]-XR[X2], XII[X2]=XII[X1]-XII[X2];
xue@1 1090 XR[X1]=Temp.x, XII[X1]=Temp.y;
xue@1 1091 }
xue@1 1092 X0=X0+(ElemsPerGroup<<1);
xue@1 1093 }
xue@1 1094 }
xue@1 1095 zp.x=XR[0]+XII[0], zn.x=XR[0]-XII[0];
xue@1 1096 XR[0]=zp.x;
xue@1 1097 XI[0]=0;
xue@1 1098 XR[N]=zn.x;
xue@1 1099 XI[N]=0;
xue@1 1100 for (int k=1; k<N/2; k++)
xue@1 1101 {
xue@1 1102 zp.x=XR[k]+XR[N-k], zn.x=XR[k]-XR[N-k],
xue@1 1103 zp.y=XII[k]+XII[N-k], zn.y=XII[k]-XII[N-k];
xue@1 1104 XR[k]=0.5*(zp.x+W[k].y*zn.x+W[k].x*zp.y);
xue@1 1105 XI[k]=0.5*(zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1 1106 XR[N-k]=0.5*(zp.x-W[k].y*zn.x-W[k].x*zp.y);
xue@1 1107 XI[N-k]=0.5*(-zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1 1108 }
xue@1 1109 XR[N/2]=XR[N/2];
xue@1 1110 XI[N/2]=-XII[N/2];
xue@1 1111 N*=2;
xue@1 1112
xue@1 1113 for (int k=N/2+1; k<N; k++) XR[k]=XR[N-k], XI[k]=-XI[N-k];
xue@1 1114 if (Amp) for (i=0; i<N; i++) Amp[i]=sqrt(XR[i]*XR[i]+XI[i]*XI[i]);
xue@1 1115 if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(XI[i], XR[i]);
xue@1 1116 }//RFFTC_ual_old
xue@1 1117
xue@1 1118 void CIFFTR_old(cdouble* Input, int Order, cdouble* W, double* X, int* bitinv)
xue@1 1119 {
xue@1 1120 int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2, *bitinv1=bitinv;
xue@1 1121 cdouble Temp;
xue@1 1122 if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1 1123
xue@1 1124 Order--;
xue@1 1125 N/=2;
xue@1 1126 double* XII=&X[N];
xue@1 1127
xue@1 1128 X[0]=0.5*(Input[0].x+Input[N].x);
xue@1 1129 XII[0]=0.5*(Input[0].x-Input[N].x);
xue@1 1130 for (int i=1; i<N/2; i++)
xue@1 1131 {
xue@1 1132 double frp=Input[i].x+Input[N-i].x, frn=Input[i].x-Input[N-i].x,
xue@1 1133 fip=Input[i].y+Input[N-i].y, fin=Input[i].y-Input[N-i].y;
xue@1 1134 X[i]=0.5*(frp+frn*W[i].y-fip*W[i].x);
xue@1 1135 XII[i]=0.5*(fin+frn*W[i].x+fip*W[i].y);
xue@1 1136 X[N-i]=0.5*(frp-frn*W[i].y+fip*W[i].x);
xue@1 1137 XII[N-i]=0.5*(-fin+frn*W[i].x+fip*W[i].y);
xue@1 1138 }
xue@1 1139 X[N/2]=Input[N/2].x;
xue@1 1140 XII[N/2]=-Input[N/2].y;
xue@1 1141
xue@1 1142 ElemsPerGroup=1<<Order;
xue@1 1143 Groups=1;
xue@1 1144
xue@1 1145 for (i=0; i<Order; i++)
xue@1 1146 {
xue@1 1147 ElemsPerGroup/=2;
xue@1 1148 X0=0;
xue@1 1149 for (j=0; j<Groups; j++)
xue@1 1150 {
xue@1 1151 int kGroups=bitinv[j]/2;
xue@1 1152 for (k=0; k<ElemsPerGroup; k++)
xue@1 1153 {
xue@1 1154 X1=X0+k;
xue@1 1155 X2=X1+ElemsPerGroup;
xue@1 1156 Temp.x=X[X2]*W[kGroups].x+XII[X2]*W[kGroups].y,
xue@1 1157 XII[X2]=-X[X2]*W[kGroups].y+XII[X2]*W[kGroups].x;
xue@1 1158 X[X2]=Temp.x;
xue@1 1159 Temp.x=X[X1]+X[X2], Temp.y=XII[X1]+XII[X2];
xue@1 1160 X[X2]=X[X1]-X[X2], XII[X2]=XII[X1]-XII[X2];
xue@1 1161 X[X1]=Temp.x, XII[X1]=Temp.y;
xue@1 1162 }
xue@1 1163 X0=X0+(ElemsPerGroup<<1);
xue@1 1164 }
xue@1 1165 Groups*=2;
xue@1 1166 }
xue@1 1167
xue@1 1168 N*=2;
xue@1 1169 Order++;
xue@1 1170 for (i=0; i<N; i++)
xue@1 1171 {
xue@1 1172 int jj=bitinv[i];
xue@1 1173 if (i<jj)
xue@1 1174 {
xue@1 1175 Temp.x=X[i];
xue@1 1176 X[i]=X[jj];
xue@1 1177 X[jj]=Temp.x;
xue@1 1178 }
xue@1 1179 }
xue@1 1180 for (int i=0; i<N; i++) X[i]/=(N/2);
xue@1 1181 if (!bitinv1) free(bitinv);
xue@1 1182 }//RFFTC_ual_old
xue@1 1183