xue@11: /*
xue@11:     Harmonic sinusoidal modelling and tools
xue@11: 
xue@11:     C++ code package for harmonic sinusoidal modelling and relevant signal processing.
xue@11:     Centre for Digital Music, Queen Mary, University of London.
xue@11:     This file copyright 2011 Wen Xue.
xue@11: 
xue@11:     This program is free software; you can redistribute it and/or
xue@11:     modify it under the terms of the GNU General Public License as
xue@11:     published by the Free Software Foundation; either version 2 of the
xue@11:     License, or (at your option) any later version.
xue@11: */
xue@1: //---------------------------------------------------------------------------
xue@1: 
Chris@2: #include <string.h>
xue@1: #include <stdlib.h>
xue@1: #include "fft.h"
xue@1: 
Chris@5: /** \file fft.h */
Chris@5: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function Atan2: (0, 0)-safe atan2
xue@1: 
xue@1:   Returns 0 is x=y=0, atan2(x,y) otherwise.
xue@1: */
xue@1: double Atan2(double y, double x)
xue@1: {
xue@1:   if (x==0 && y==0) return 0;
xue@1:   else return atan2(y, x);
xue@1: }//Atan2
xue@1: 
Chris@5: /**
xue@11:   function Log2: Log2
xue@11: 
xue@11:   Returns x for Log2(2^x) if x is integer.
xue@11: */
xue@11: int Log2(int x)
xue@11: {
xue@11:   return floor(log(x)/log(2)+0.5);
xue@11: }//Log2
xue@11: 
xue@11: /**
xue@1:   function BitInv: inverse bit order of Value within an $Order-bit expression.
xue@1: 
xue@1:   In: integer Value smaller than 2^Order
xue@1: 
xue@1:   Returns an integer whose lowest Order bits are the lowest Order bits of Value in reverse order.
xue@1: */
xue@1: int BitInv(int Value, int Order)
xue@1: {
xue@1:   int Result;
xue@1:   Result=0;
xue@1:   for (int i=0;i<Order;i++)
xue@1:   {
xue@1:     Result=(Result<<1)+(Value&0x00000001);
xue@1:     Value=Value>>1;
xue@1:   }
xue@1:   return Result;
xue@1: }//BitInv
xue@1: 
Chris@5: /**
xue@1:   function SetTwiddleFactors: fill w[N/2] with twiddle factors used in N-point complex FFT.
xue@1: 
xue@1:   In: N
xue@1:   Out: array w[N/2] containing twiddle factors
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void SetTwiddleFactors(int N, cdouble* w)
xue@1: {
xue@1:   double ep=-M_PI*2/N;
xue@1:   for (int i=0; i<N/2; i++)
xue@1:   {
xue@1:     double tmp=ep*i;
xue@1:     w[i].x=cos(tmp), w[i].y=sin(tmp);
xue@1:   }
xue@1: }//SetTwiddleFactors
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CFFTCbii: basic complex DIF-FFT module, applied after bit-inversed ordering of inputs
xue@1: 
xue@1:   In: Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex waveform
xue@1:   Out: X[Wid]: complex spectrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CFFTCbii(int Order, cdouble* W, cdouble* X)
xue@1: {
xue@1:   int i, j, k, ElemsPerGroup, Groups, X0, X1, X2;
xue@1:   cdouble Temp;
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup=1<<i;
xue@1:     Groups=1<<(Order-i-1);
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         int kGroups=k*Groups;
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         //X(X2)<-X(X2)*W
xue@1:         Temp.x=X[X2].x*W[kGroups].x-X[X2].y*W[kGroups].y,
xue@1:         X[X2].y=X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1:         X[X2].x=Temp.x;
xue@1:         Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1:         X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1:         X[X1]=Temp;
xue@1:       }
xue@1:       X0+=ElemsPerGroup*2;
xue@1:     }
xue@1:   }
xue@1: }//CFFTCbii
xue@1: 
Chris@5: /**
xue@1:   function CFFTC: in-place complex FFT
xue@1: 
xue@1:   In: Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex waveform
xue@1:       bitinv[Wid]:  bit-inversion table
xue@1:   Out: X[Wid]: complex spectrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CFFTC(int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order, i, jj;
xue@1:   cdouble Temp;
xue@1:   int* bitinv1=bitinv;
xue@1:   if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1:   for (i=1; i<N-1; i++)
xue@1:   {
xue@1:     jj=bitinv[i];
xue@1:     if (i<jj)
xue@1:     {
xue@1:       Temp=X[i];
xue@1:       X[i]=X[jj];
xue@1:       X[jj]=Temp;
xue@1:     }
xue@1:   }
xue@1:   if (!bitinv1) free(bitinv);
xue@1:   CFFTCbii(Order, W, X);
xue@1: }//CFFTC
xue@1: 
Chris@5: /**
xue@1:   function CFFTC: complex FFT
xue@1: 
xue@1:   In: Input[Wid]: complex waveform
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       bitinv[Wid]:  bit-inversion table
xue@1:   Out:X[Wid]: complex spectrum
xue@1:       Amp[Wid]: amplitude spectrum
xue@1:       Arg[Wid]: phase spectrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CFFTC(cdouble* Input, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int i, N=1<<Order;
xue@1: 
xue@1:   if (X!=Input) memcpy(X, Input, sizeof(cdouble)*N);
xue@1:   CFFTC(Order, W, X, bitinv);
xue@1: 
xue@1:   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:   if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(X[i].y, X[i].x);
xue@1: }//CFFTC
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CIFFTCbii: basic complex IFFT module, applied after bit-inversed ordering of inputs
xue@1: 
xue@1:   In: Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex spectrum
xue@1:   Out: X[Wid]: complex waveform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CIFFTCbii(int Order, cdouble* W, cdouble* X)
xue@1: {
xue@1:   int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2;
xue@1:   cdouble Temp;
xue@1: 
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup=1<<i;
xue@1:     Groups=1<<(Order-i-1);
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         int kGroups=k*Groups;
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         Temp.x=X[X2].x*W[kGroups].x+X[X2].y*W[kGroups].y,
xue@1:         X[X2].y=-X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1:         X[X2].x=Temp.x;
xue@1:         Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1:         X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1:         X[X1]=Temp;
xue@1:       }
xue@1:       X0=X0+ElemsPerGroup*2;
xue@1:     }
xue@1:   }
xue@1:   for (i=0; i<N; i++)
xue@1:   {
xue@1:     X[i].x/=N;
xue@1:     X[i].y/=N;
xue@1:   }
xue@1: }//CIFFTCbii
xue@1: 
Chris@5: /**
xue@1:   function CIFFTC: in-place complex IFFT
xue@1: 
xue@1:   In: Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex spectrum
xue@1:       bitinv[Wid]:  bit-inversion table
xue@1:   Out: X[Wid]: complex waveform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CIFFTC(int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order, i, jj, *bitinv1=bitinv;
xue@1:   cdouble Temp;
xue@1:   if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1:   for (i=1; i<N-1; i++)
xue@1:   {
xue@1:     jj=bitinv[i];
xue@1:     if (i<jj)
xue@1:     {
xue@1:       Temp=X[i];
xue@1:       X[i]=X[jj];
xue@1:       X[jj]=Temp;
xue@1:     }
xue@1:   }
xue@1:   if (!bitinv1) free(bitinv);
xue@1:   CIFFTCbii(Order, W, X);
xue@1: }//CIFFTC
xue@1: 
Chris@5: /**
xue@1:   function CIFFTC: complex IFFT
xue@1: 
xue@1:   In: Input[Wid]: complex spectrum
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       bitinv[Wid]:  bit-inversion table
xue@1:   Out:X[Wid]: complex waveform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CIFFTC(cdouble* Input, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order;
xue@1:   if (X!=Input) memcpy(X, Input, sizeof(cdouble)*N);
xue@1:   if (bitinv) CIFFTC(Order, W, X, bitinv);
xue@1:   else CIFFTC(Order, W, X);
xue@1: }//CIFFTC
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CIFFTR: complex-to-real IFFT
xue@1: 
xue@1:   In: Input[Wid/2+1]: complex spectrum, imaginary parts of Input[0] and Input[Wid/2] are ignored
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:X[Wid]: real waveform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CIFFTR(cdouble* Input, int Order, cdouble* W, double* X, int* hbitinv)
xue@1: {
xue@1:   int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2, *hbitinv1=hbitinv;
xue@1:   cdouble Temp;
xue@1: 
xue@1:   Order--; N/=2;
xue@1:   if (!hbitinv) hbitinv=CreateBitInvTable(Order);
xue@1: 
xue@1:   cdouble* Xc=(cdouble*)X;
xue@1: 
xue@1:   Xc[0].y=0.5*(Input[0].x-Input[N].x);
xue@1:   Xc[0].x=0.5*(Input[0].x+Input[N].x);
xue@1:   for (int i=1; i<N/2; i++)
xue@1:   {
xue@1:     double frp=Input[i].x+Input[N-i].x, frn=Input[i].x-Input[N-i].x,
xue@1:            fip=Input[i].y+Input[N-i].y, fin=Input[i].y-Input[N-i].y;
xue@1:     Xc[i].x=0.5*(frp+frn*W[i].y-fip*W[i].x);
xue@1:     Xc[i].y=0.5*(fin+frn*W[i].x+fip*W[i].y);
xue@1:     Xc[N-i].x=0.5*(frp-frn*W[i].y+fip*W[i].x);
xue@1:     Xc[N-i].y=0.5*(-fin+frn*W[i].x+fip*W[i].y);
xue@1:   }
xue@1:   Xc[N/2].x=Input[N/2].x;
xue@1:   Xc[N/2].y=-Input[N/2].y;
xue@1: 
xue@1:   ElemsPerGroup=1<<Order;
xue@1:   Groups=1;
xue@1: 
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup/=2;
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       int kGroups=hbitinv[j];
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         Temp.x=Xc[X2].x*W[kGroups].x+Xc[X2].y*W[kGroups].y,
xue@1:         Xc[X2].y=-Xc[X2].x*W[kGroups].y+Xc[X2].y*W[kGroups].x;
xue@1:         Xc[X2].x=Temp.x;
xue@1:         Temp.x=Xc[X1].x+Xc[X2].x, Temp.y=Xc[X1].y+Xc[X2].y;
xue@1:         Xc[X2].x=Xc[X1].x-Xc[X2].x, Xc[X2].y=Xc[X1].y-Xc[X2].y;
xue@1:         Xc[X1].x=Temp.x, Xc[X1].y=Temp.y;
xue@1:       }
xue@1:       X0=X0+(ElemsPerGroup<<1);
xue@1:     }
xue@1:     Groups*=2;
xue@1:   }
xue@1: 
xue@1:   for (i=0; i<N; i++)
xue@1:   {
xue@1:     int jj=hbitinv[i];
xue@1:     if (i<jj)
xue@1:     {
xue@1:       Temp=Xc[i];
xue@1:       Xc[i]=Xc[jj];
xue@1:       Xc[jj]=Temp;
xue@1:     }
xue@1:   }
xue@1:   double norm=1.0/N;;
xue@1:   N*=2;
xue@1:   for (int i=0; i<N; i++) X[i]*=norm;
xue@1:   if (!hbitinv1) free(hbitinv);
xue@1: }//CIFFTR
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CreateBitInvTable: creates a table of bit-inversed integers
xue@1: 
xue@1:   In: Order: interger, equals log2(size of table)
xue@1: 
xue@1:   Returns a pointer to a newly allocated array containing the table. The returned pointer must be freed
xue@1:   using free(), NOT delete[].
xue@1: */
xue@1: int* CreateBitInvTable(int Order)
xue@1: {
xue@1:   int N=1<<Order;
xue@1:   int* result=(int*)malloc(sizeof(int)*N);
xue@1:   for (int i=0; i<N; i++) result[i]=BitInv(i, Order);
xue@1:   return result;
xue@1: }//CreateBitInvTable
xue@1: 
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RFFTC_ual: unaligned real-to-complex FFT
xue@1: 
xue@1:   In: Input[Wid]: real waveform
xue@1:       Order; integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:X[Wid}: complex spectrum
xue@1:       Amp[Wid]: amplitude spectrum
xue@1:       Arg[Wid]: phase spetrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RFFTC_ual(double* Input, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1: {
xue@1:   int N=1<<Order, i, j, k, *hbitinv1=hbitinv, Groups, ElemsPerGroup, X0, X1, X2;
xue@1:   cdouble Temp, zp, zn;
xue@1: 
xue@1:   N/=2; Order--;
xue@1: 
xue@1:   //Input being NULL implies external initialization of X. This is used by RFFTCW but is not
xue@1:   //recommended for external use.
xue@1:   if (Input)
xue@1:   {
xue@1:     if (!hbitinv) hbitinv=CreateBitInvTable(Order);
xue@1: 
xue@1:     if (Input==(double*)X)
xue@1:     {
xue@1:       //Input being identical to X is not recommended for external use.
xue@1:       for (int i=0; i<N; i++)
xue@1:       {
xue@1:         int bi=hbitinv[i];
xue@1:         if (i<bi) {cdouble tmp=X[i]; X[i]=X[bi]; X[bi]=tmp;}
xue@1:       }
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       for (i=0; i<N; i++) X[i]=((cdouble*)Input)[hbitinv[i]];
xue@1:     }
xue@1:     if (!hbitinv1) free(hbitinv);
xue@1:   }
xue@1: 
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup=1<<i;
xue@1:     Groups=1<<(Order-i-1);
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         int kGroups=k*2*Groups;
xue@1:         Temp.x=X[X2].x*W[kGroups].x-X[X2].y*W[kGroups].y,
xue@1:         X[X2].y=X[X2].x*W[kGroups].y+X[X2].y*W[kGroups].x;
xue@1:         X[X2].x=Temp.x;
xue@1:         Temp.x=X[X1].x+X[X2].x, Temp.y=X[X1].y+X[X2].y;
xue@1:         X[X2].x=X[X1].x-X[X2].x, X[X2].y=X[X1].y-X[X2].y;
xue@1:         X[X1]=Temp;
xue@1:       }
xue@1:       X0=X0+(ElemsPerGroup<<1);
xue@1:     }
xue@1:   }
xue@1:   zp.x=X[0].x+X[0].y, zn.x=X[0].x-X[0].y;
xue@1:   X[0].x=zp.x;
xue@1:   X[0].y=0;
xue@1:   X[N].x=zn.x;
xue@1:   X[N].y=0;
xue@1:   for (int k=1; k<N/2; k++)
xue@1:   {
xue@1:     zp.x=X[k].x+X[N-k].x, zn.x=X[k].x-X[N-k].x,
xue@1:     zp.y=X[k].y+X[N-k].y, zn.y=X[k].y-X[N-k].y;
xue@1:     X[k].x=0.5*(zp.x+W[k].y*zn.x+W[k].x*zp.y);
xue@1:     X[k].y=0.5*(zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1:     X[N-k].x=0.5*(zp.x-W[k].y*zn.x-W[k].x*zp.y);
xue@1:     X[N-k].y=0.5*(-zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1:   }
xue@1:   //X[N/2].x=X[N/2].x;
xue@1:   X[N/2].y=-X[N/2].y;
xue@1:   N*=2;
xue@1: 
xue@1:   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:   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:   if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(X[i].x, X[i].y);
xue@1: }//RFFTC_ual
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RFFTCW: real-to-complex FFT with window
xue@1: 
xue@1:   In: Input[Wid]: real waveform
xue@1:       Window[Wid]:  window function
xue@1:       Order; integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:X[Wid}: complex spectrum
xue@1:       Amp[Wid]: amplitude spectrum
xue@1:       Arg[Wid]: phase spetrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RFFTCW(double* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1: {
xue@1:   int N=1<<Order, *hbitinv1=hbitinv;
xue@1:   if (!hbitinv) hbitinv=CreateBitInvTable(Order-1);
xue@1:   N/=2; 
xue@1: 
xue@1:   if (Input==(double*)X)
xue@1:   { //so that X[n].x IS Input[2n], X[n].y IS Input[2n+1]
xue@1:     for (int n=0; n<N; n++)
xue@1:     {
xue@1:       int bi=hbitinv[n], n2=n+n, bi2=bi+bi;
xue@1:       if (n<bi)
xue@1:       {
xue@1:         double tmp=X[n].x*Window[n2]; X[n].x=X[bi].x*Window[bi2]; X[bi].x=tmp;
xue@1:         tmp=X[n].y*Window[n2+1]; X[n].y=X[bi].y*Window[bi2+1]; X[bi].y=tmp;
xue@1:       }
xue@1:       else if (n==bi)
xue@1:       { //so that X[n].x IS Input[bi]
xue@1:         X[n].x*=Window[bi2], X[n].y*=Window[bi2+1];
xue@1:       }
xue@1:     }
xue@1:   }
xue@1:   else
xue@1:   {
xue@1:     for (int n=0; n<N; n++)
xue@1:     {
xue@1:       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:     }
xue@1:   }
xue@1: 
xue@1:   RFFTC_ual(0, Amp, Arg, Order, W, X, hbitinv);
xue@1:   if (!hbitinv1) free(hbitinv);
xue@1: }//RFFTCW
xue@1: 
Chris@5: /**
xue@1:   function RFFTCW: real-to-complex FFT with window
xue@1: 
xue@1:   In: Input[Wid]: real waveform as _int16 array
xue@1:       Window[Wid]:  window function
xue@1:       Order; integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:X[Wid}: complex spectrum
xue@1:       Amp[Wid]: amplitude spectrum
xue@1:       Arg[Wid]: phase spetrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RFFTCW(__int16* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv)
xue@1: {
xue@1:   int N=1<<Order, *hbitinv1=hbitinv;
xue@1: 
xue@1:   N/=2;
xue@1:   if (!hbitinv) hbitinv=CreateBitInvTable(Order-1);
xue@1:   for (int n=0; n<N; n++)
xue@1:   {
xue@1:     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:   }
xue@1: 
xue@1:   RFFTC_ual(0, Amp, Arg, Order, W, X, hbitinv);
xue@1:   if (!hbitinv1) free(hbitinv);
xue@1: }//RFFTCW
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CFFTCW: complex FFT with window
xue@1: 
xue@1:   In: Window[Wid]: window function
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex waveform
xue@1:       bitinv[Wid]: bit-inversion table
xue@1:   Out:X[Wid], complex spectrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CFFTCW(double* Window, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order;
xue@1:   for (int i=0; i<N; i++) X[i].x*=Window[i], X[i].y*=Window[i];
xue@1:   CFFTC(Order, W, X, bitinv);
xue@1: }//CFFTCW
xue@1: 
Chris@5: /**
xue@1:   function CFFTCW: complex FFT with window
xue@1: 
xue@1:   In: Input[Wid]: complex waveform
xue@1:       Window[Wid]: window function
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:       X[Wid]: complex waveform
xue@1:       bitinv[Wid]: bit-inversion table
xue@1:   Out:X[Wid], complex spectrum
xue@1:       Amp[Wid], amplitude spectrum
xue@1:       Arg[Wid], phase spectrum
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void CFFTCW(cdouble* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order;
xue@1:   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:   CFFTC(X, Amp, Arg, Order, W, X, bitinv);
xue@1: }//CFFTCW
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RDCT1: fast DCT1 implemented using FFT. DCT-I has the time scale 0.5-sample shifted from the DFT.
xue@1: 
xue@1:   In: Input[Wid]: real waveform
xue@1:       Order: integer, equals log2(Wid)
xue@1:       qW[Wid/8]: quarter table of twiddle factors
xue@1:       qX[Wid/4]: quarter FFT data buffer
xue@1:       qbitinv[Wid/4]: quarter bit-inversion table
xue@1:   Out:Output[Wid]: DCT-I of Input.
xue@1: 
xue@1:   No return value. Content of qW is destroyed on return. On calling the fft buffers should be allocated
xue@1:   to size 0.25*Wid.
xue@1: */
xue@1: void RDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv)
xue@1: {
xue@1:   const double lmd0=sqrt(0.5);
xue@1:   if (Order==0)
xue@1:   {
xue@1:     Output[0]=Input[0]*lmd0;
xue@1:     return;
xue@1:   }
xue@1:   if (Order==1)
xue@1:   {
xue@1:     double tmp1=(Input[0]+Input[1])*lmd0;
xue@1:     Output[1]=(Input[0]-Input[1])*lmd0;
xue@1:     Output[0]=tmp1;
xue@1:     return;
xue@1:   }
xue@1:   int order=Order-1, N=1<<(Order-1), C=1;
xue@1:   bool createbitinv=!qbitinv;
xue@1:   if (createbitinv) qbitinv=CreateBitInvTable(Order-2);
xue@1:   double *even=(double*)malloc8(sizeof(double)*N*2);
xue@1:   double *odd=&even[N];
xue@1:   //data pass from Input to Output, combined with the first sequence split
xue@1:   for (int i=0, N2=N*2; i<N; i++)
xue@1:   {
xue@1:     even[i]=Input[i]+Input[N2-1-i];
xue@1:     odd[i]=Input[i]-Input[N2-1-i];
xue@1:   }
xue@1:   for (int i=0; i<N; i++) Output[i*2]=even[i], Output[i*2+1]=odd[i];
xue@1:   while (order>1)
xue@1:   {
xue@1:   //N=2^order, 4|N, 2|hN
xue@1:   //keep subsequence 0, 2C, 4C, ... 2(N-1)C
xue@1:   //process sequence C, 3C, ..., (2N-1)C only
xue@1:     //RDCT4 routine
xue@1:     int hN=N/2, N2=N*2;
xue@1:     for (int i=0; i<hN; i++)
xue@1:     {
xue@1:       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:       double co=cos(theta), si=sin(theta);
xue@1:       qX[i].x=b*co-c*si, qX[i].y=c*co+b*si;
xue@1:     }
xue@1:     CFFTC(order-1, qW, qX, qbitinv);
xue@1:     for (int i=0; i<hN; i++)
xue@1:     {
xue@1:       double gr=qX[i].x, gi=qX[i].y, theta=-i*M_PI/N;
xue@1:       double co=cos(theta), si=sin(theta);
xue@1:       Output[(4*i+1)*C]=gr*co-gi*si;
xue@1:       Output[(N2-4*i-1)*C]=-gr*si-gi*co;
xue@1:     }
xue@1:     N=(N>>1);
xue@1:     C=(C<<1);
xue@1:     for (int i=1; i<N/4; i++)
xue@1:       qW[i]=qW[i*2];
xue@1:     for (int i=1; i<N/2; i++)
xue@1:       qbitinv[i]=qbitinv[i*2];
xue@1: 
xue@1:     //splitting subsequence 0, 2C, 4C, ..., 2(N-1)C
xue@1:     for (int i=0, N2=N*2; i<N; i++)
xue@1:     {
xue@1:       even[i]=Output[i*C]+Output[(N2-1-i)*C];
xue@1:       odd[i]=Output[i*C]-Output[(N2-1-i)*C];
xue@1:     }
xue@1:     for (int i=0; i<N; i++)
xue@1:     {
xue@1:       Output[2*i*C]=even[i];
xue@1:       Output[(2*i+1)*C]=odd[i];
xue@1:     }
xue@1:     order--;
xue@1:   }
xue@1:   //order==1
xue@1:   //use C and 3C in DCT4
xue@1:   //use 0 and 2C in DCT1
xue@1:   double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1:   double tmp1=c1*Output[C]+c2*Output[3*C];
xue@1:   Output[3*C]=c2*Output[C]-c1*Output[3*C];
xue@1:   Output[C]=tmp1;
xue@1:   tmp1=Output[0]+Output[2*C];
xue@1:   Output[2*C]=(Output[0]-Output[2*C])*lmd0;
xue@1:   Output[0]=tmp1*lmd0;
xue@1: 
xue@1:   if (createbitinv) free(qbitinv);
xue@1:   free8(even);
xue@1: }//RDCT1
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RDCT4: fast DCT4 implemented using FFT. DCT-IV has both the time and frequency scales
xue@1:   0.5-sample or 0.5-bin shifted from DFT.
xue@1: 
xue@1:   In: Input[Wid]: real waveform
xue@1:       Order: integer, equals log2(Wid)
xue@1:       hW[Wid/4]: half table of twiddle factors
xue@1:       hX[Wid/2]: hal FFT data buffer
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:Output[Wid]: DCT-IV of Input.
xue@1: 
xue@1:   No return value. Content of hW not affected. On calling the fft buffers should be allocated to size
xue@1:   0.5*Wid.
xue@1: */
xue@1: void RDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv)
xue@1: {
xue@1:   if (Order==0)
xue@1:   {
xue@1:     Output[0]=Input[0]/sqrt(2);
xue@1:     return;
xue@1:   }
xue@1:   if (Order==1)
xue@1:   {
xue@1:     double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1:     double tmp1=c1*Input[0]+c2*Input[1];
xue@1:     Output[1]=c2*Input[0]-c1*Input[1];
xue@1:     Output[0]=tmp1;
xue@1:     return;
xue@1:   }
xue@1:   int N=1<<Order, hN=N/2;
xue@1:   for (int i=0; i<hN; i++)
xue@1:   {
xue@1:     double b=Input[2*i], c=Input[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1:     double co=cos(theta), si=sin(theta);
xue@1:     hX[i].x=b*co-c*si, hX[i].y=c*co+b*si;
xue@1:   }
xue@1:   CFFTC(Order-1, hW, hX, hbitinv);
xue@1:   for (int i=0; i<hN; i++)
xue@1:   {
xue@1:     double gr=hX[i].x, gi=hX[i].y, theta=-i*M_PI/N;
xue@1:     double co=cos(theta), si=sin(theta);
xue@1:     Output[2*i]=gr*co-gi*si;
xue@1:     Output[N-1-2*i]=-gr*si-gi*co;
xue@1:   }
xue@1: }//RDCT4
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RIDCT1: fast IDCT1 implemented using FFT.
xue@1: 
xue@1:   In: Input[Wid]: DCT-I of some real waveform.
xue@1:       Order: integer, equals log2(Wid)
xue@1:       qW[Wid/8]: quarter table of twiddle factors
xue@1:       qX[Wid/4]: quarter FFT data buffer
xue@1:       qbitinv[Wid/4]: quarter bit-inversion table
xue@1:   Out:Output[Wid]: IDCT-I of Input.
xue@1: 
xue@1:   No return value. Content of qW is destroyed on return. On calling the fft buffers should be allocated
xue@1:   to size 0.25*Wid.
xue@1: */
xue@1: void RIDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv)
xue@1: {
xue@1:   const double lmd0=sqrt(0.5);
xue@1:   if (Order==0)
xue@1:   {
xue@1:     Output[0]=Input[0]/lmd0;
xue@1:     return;
xue@1:   }
xue@1:   if (Order==1)
xue@1:   {
xue@1:     double tmp1=(Input[0]+Input[1])*lmd0;
xue@1:     Output[1]=(Input[0]-Input[1])*lmd0;
xue@1:     Output[0]=tmp1;
xue@1:     return;
xue@1:   }
xue@1:   int order=Order-1, N=1<<(Order-1), C=1;
xue@1:   bool createbitinv=!qbitinv;
xue@1:   if (createbitinv) qbitinv=CreateBitInvTable(Order-2);
xue@1:   double *even=(double*)malloc8(sizeof(double)*N*2);
xue@1:   double *odd=&even[N];
xue@1: 
xue@1:   while (order>0)
xue@1:   {
xue@1:     //N=2^order, 4|N, 2|hN
xue@1:     //keep subsequence 0, 2C, 4C, ... 2(N-1)C
xue@1:     //process sequence C, 3C, ..., (2N-1)C only
xue@1:     //data pass from Input
xue@1:     for (int i=0; i<N; i++)
xue@1:     {
xue@1:       odd[i]=Input[(i*2+1)*C];
xue@1:     }
xue@1: 
xue@1:     //IDCT4 routine
xue@1:     //RIDCT4(odd, odd, order, qW, qX, qbitinv);
xue@1: 
xue@1:       if (order==1)
xue@1:       {
xue@1:         double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1:         double tmp1=c1*odd[0]+c2*odd[1];
xue@1:         odd[1]=c2*odd[0]-c1*odd[1];
xue@1:         odd[0]=tmp1;
xue@1:       }
xue@1:       else
xue@1:       {
xue@1:         int hN=N/2;
xue@1:         for (int i=0; i<hN; i++)
xue@1:         {
xue@1:           double b=odd[2*i], c=odd[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1:           double co=cos(theta), si=sin(theta);
xue@1:           qX[i].x=b*co-c*si, qX[i].y=c*co+b*si;
xue@1:         }
xue@1:         CFFTC(order-1, qW, qX, qbitinv);
xue@1:         double i2N=2.0/N;
xue@1:         for (int i=0; i<hN; i++)
xue@1:         {
xue@1:           double gr=qX[i].x, gi=qX[i].y, theta=-i*M_PI/N;
xue@1:           double co=cos(theta), si=sin(theta);
xue@1:           odd[2*i]=(gr*co-gi*si)*i2N;
xue@1:           odd[N-1-2*i]=(-gr*si-gi*co)*i2N;
xue@1:         }
xue@1:       }
xue@1: 
xue@1:     order--;
xue@1:     N=(N>>1);
xue@1:     C=(C<<1);
xue@1:     for (int i=1; i<N/4; i++)
xue@1:       qW[i]=qW[i*2];
xue@1:     for (int i=1; i<N/2; i++)
xue@1:       qbitinv[i]=qbitinv[i*2];
xue@1: 
xue@1:     odd=&even[N];
xue@1:   }
xue@1:   //order==0
xue@1:   even[0]=Input[0];
xue@1:   even[1]=Input[C];
xue@1:   double tmp1=(even[0]+even[1])*lmd0;
xue@1:   Output[1]=(even[0]-even[1])*lmd0;
xue@1:   Output[0]=tmp1;
xue@1:   order++;
xue@1: 
xue@1:   while (order<Order)
xue@1:   {
xue@1:     N=(N<<1);
xue@1:     odd=&even[N];
xue@1:     for (int i=0; i<N; i++)
xue@1:     {
xue@1:       Output[N*2-1-i]=(Output[i]-odd[i])/2;
xue@1:       Output[i]=(Output[i]+odd[i])/2;
xue@1:     }
xue@1:     order++;
xue@1:   }
xue@1:   
xue@1:   if (createbitinv) free(qbitinv);
xue@1:   free8(even);
xue@1: }//RIDCT1
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RIDCT4: fast IDCT4 implemented using FFT.
xue@1: 
xue@1:   In: Input[Wid]: DCT-IV of some real waveform
xue@1:       Order: integer, equals log2(Wid)
xue@1:       hW[Wid/4]: half table of twiddle factors
xue@1:       hX[Wid/2]: hal FFT data buffer
xue@1:       hbitinv[Wid/2]: half bit-inversion table
xue@1:   Out:Output[Wid]: IDCT-IV of Input.
xue@1: 
xue@1:   No return value. Content of hW not affected. On calling the fft buffers should be allocated to size
xue@1:   0.5*Wid.
xue@1: */
xue@1: void RIDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv)
xue@1: {
xue@1:   if (Order==0)
xue@1:   {
xue@1:     Output[0]=Input[0]*sqrt(2);
xue@1:     return;
xue@1:   }
xue@1:   if (Order==1)
xue@1:   {
xue@1:     double c1=cos(M_PI/8), c2=cos(3*M_PI/8);
xue@1:     double tmp1=c1*Input[0]+c2*Input[1];
xue@1:     Output[1]=c2*Input[0]-c1*Input[1];
xue@1:     Output[0]=tmp1;
xue@1:     return;
xue@1:   }
xue@1:   int N=1<<Order, hN=N/2;
xue@1:   for (int i=0; i<hN; i++)
xue@1:   {
xue@1:     double b=Input[2*i], c=Input[N-1-i*2], theta=-(i+0.25)*M_PI/N;
xue@1:     double co=cos(theta), si=sin(theta);
xue@1:     hX[i].x=b*co-c*si, hX[i].y=c*co+b*si;
xue@1:   }
xue@1:   CFFTC(Order-1, hW, hX, hbitinv);
xue@1:   double i2N=2.0/N;
xue@1:   for (int i=0; i<hN; i++)
xue@1:   {
xue@1:     double gr=hX[i].x, gi=hX[i].y, theta=-i*M_PI/N;
xue@1:     double co=cos(theta), si=sin(theta);
xue@1:     Output[2*i]=(gr*co-gi*si)*i2N;
xue@1:     Output[N-1-2*i]=(-gr*si-gi*co)*i2N;
xue@1:   }
xue@1: }//RIDCT4
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RLCT: real local cosine transform of equal frame widths Wid=2^Order
xue@1: 
xue@1:   In: data[Count]: real waveform
xue@1:       Order: integer, equals log2(Wid), Wid being the cosine atom size
xue@1:       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:         example: wid=4, g[k]=sin(pi*(k+0.5)/8).
xue@1:   Out:spec[Fr][Wid]: the local cosine transform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RLCT(double** spec, double* data, int Count, int Order, int wid, double* g)
xue@1: {
xue@1:   int Wid=1<<Order, Fr=Count/Wid, hwid=wid/2;
xue@1:   int* hbitinv=CreateBitInvTable(Order-1);
xue@1:   cdouble *hx=(cdouble*)malloc8(sizeof(cdouble)*Wid*3/4), *hw=(cdouble*)&hx[Wid/2];
xue@1:   double norm=sqrt(2.0/Wid);
xue@1:   SetTwiddleFactors(Wid/2, hw);
xue@1: 
xue@1:   for (int fr=0; fr<Fr; fr++)
xue@1:   {
xue@1:     if (fr==0)
xue@1:     {
xue@1:       memcpy(spec[fr], data, sizeof(double)*(Wid-hwid));
xue@1:       for (int i=0, k=Wid+hwid-1, l=Wid-hwid; i<hwid; i++, k--, l++)
xue@1:         spec[fr][l]=data[l]*g[wid-1-i]-data[k]*g[i];
xue@1:     }
xue@1:     else if (fr==Fr-1)
xue@1:     {
xue@1:       for (int i=hwid, j=fr*Wid, k=fr*Wid-1, l=0; i<wid; i++, j++, k--, l++)
xue@1:         spec[fr][l]=data[k]*g[wid-1-i]+data[j]*g[i];
xue@1:       memcpy(&spec[fr][hwid], &data[fr*Wid+hwid], sizeof(double)*(Wid-hwid));
xue@1:     }
xue@1:     else
xue@1:     {
xue@1:       for (int i=hwid, j=fr*Wid, k=fr*Wid-1, l=0; i<wid; i++, j++, k--, l++)
xue@1:         spec[fr][l]=data[k]*g[wid-1-i]+data[j]*g[i];
xue@1:       if (Wid>wid) memcpy(&spec[fr][hwid], &data[fr*Wid+hwid], sizeof(double)*(Wid-wid));
xue@1:       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:         spec[fr][l]=data[j]*g[wid-1-i]-data[k]*g[i];
xue@1:     }
xue@1:   }
xue@1:   for (int fr=0; fr<Fr; fr++)
xue@1:   {
xue@1:     if (fr==Fr-1)
xue@1:     {
xue@1:       for (int i=1; i<Wid/4; i++) hw[i]=hw[2*i], hbitinv[i]=hbitinv[2*i];
xue@1:       RDCT1(spec[fr], spec[fr], Order, hw, hx, hbitinv);
xue@1:     }
xue@1:     else
xue@1:       RDCT4(spec[fr], spec[fr], Order, hw, hx, hbitinv);
xue@1: 
xue@1: ////The following line can be removed if the corresponding line in RILCT(...) is removed
xue@1:     for (int i=0; i<Wid; i++) spec[fr][i]*=norm;
xue@1:   }
xue@1:   free(hbitinv);
xue@1:   free8(hx);
xue@1: }//RLCT
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function RILCT: inverse local cosine transform of equal frame widths Wid=2^Order
xue@1: 
xue@1:   In: spec[Fr][Wid]: the local cosine transform
xue@1:       Order: integer, equals log2(Wid), Wid being the cosine atom size
xue@1:       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:         example: wid=4, g[k]=sin(pi*(k+0.5)/8).
xue@1:   Out:data[Count]: real waveform
xue@1: 
xue@1:   No return value.
xue@1: */
xue@1: void RILCT(double* data, double** spec, int Fr, int Order, int wid, double* g)
xue@1: {
xue@1:   int Wid=1<<Order, Count=Fr*Wid, hwid=wid/2, *hbitinv=CreateBitInvTable(Order-1);
xue@1:   cdouble *hx=(cdouble*)malloc8(sizeof(cdouble)*Wid*3/4), *hw=&hx[Wid/2];
xue@1:   double norm=sqrt(Wid/2.0);
xue@1:   SetTwiddleFactors(Wid/2, hw);
xue@1: 
xue@1:   for (int fr=0; fr<Fr; fr++)
xue@1:   {
xue@1:     if (fr==Fr-1)
xue@1:     {
xue@1:       for (int i=1; i<Wid/4; i++) hw[i]=hw[2*i], hbitinv[i]=hbitinv[i*2];
xue@1:       RIDCT1(spec[fr], &data[fr*Wid], Order, hw, hx, hbitinv);
xue@1:     }
xue@1:     else
xue@1:       RIDCT4(spec[fr], &data[fr*Wid], Order, hw, hx, hbitinv);
xue@1:   }
xue@1:   //unfolding
xue@1:   for (int fr=1; fr<Fr; fr++)
xue@1:   {
xue@1:     double* h=&data[fr*Wid];
xue@1:     for (int i=0; i<hwid; i++)
xue@1:     {
xue@1:       double a=h[i], b=h[-1-i], c=g[i+hwid], d=g[hwid-1-i];
xue@1:       h[i]=a*c-b*d, h[-i-1]=b*c+a*d;
xue@1:     }
xue@1:   }
xue@1: 
xue@1:   free8(hx);
xue@1: ////The following line can be removed if the corresponding line in RLCT(...) is removed
xue@1:   for (int i=0; i<Count; i++) data[i]*=norm;
xue@1: }//RILCT
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   function CMFTC: radix-2 complex multiresolution Fourier transform
xue@1: 
xue@1:   In: x[Wid]: complex waveform
xue@1:       Order: integer, equals log2(Wid)
xue@1:       W[Wid/2]: twiddle factors
xue@1:   Out:X[Order+1][Wid]: multiresolution FT of x. X[0] is the same as x itself.
xue@1: 
xue@1:   No return value.
xue@1: 
xue@1:   Further reading: Wen X. and M. Sandler, "Calculation of radix-2 discrete multiresolution Fourier
xue@1:   transform," Signal Processing, vol.87 no.10, 2007, pp.2455-2460.
xue@1: */
xue@1: void CMFTC(cdouble* x, int Order, cdouble** X, cdouble* W)
xue@1: {
xue@1:   X[0]=x;
xue@1:   for (int n=1, L=1<<(Order-1), M=2; n<=Order; n++, L>>=1, M<<=1)
xue@1:   {
xue@1:     cdouble *Xn=X[n], *Xp=X[n-1];
xue@1:     for (int l=0; l<L; l++)
xue@1:     {
xue@1:       cdouble* lX=&Xn[l*M];
xue@1:       for (int m=0; m<M/2; m++)
xue@1:       {
xue@1:         lX[m].x=Xp[l*M+m].x+Xp[l*M+M/2+m].x;
xue@1:         lX[m].y=Xp[l*M+m].y+Xp[l*M+M/2+m].y;
xue@1:         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:         int iw=m*L;
xue@1:         double wr=W[iw].x, wi=W[iw].y;
xue@1:         lX[M/2+m].x=tmpr*wr-tmpi*wi;
xue@1:         lX[M/2+m].y=tmpr*wi+tmpi*wr;
xue@1:       }
xue@1:       if (n==1) {}
xue@1:       else if (n==2) //two-point DFT
xue@1:       {
xue@1:         cdouble *aX=&X[n][l*M+M/2];
xue@1:         double tmp;
xue@1:         tmp=aX[0].x+aX[1].x; aX[1].x=aX[0].x-aX[1].x; aX[0].x=tmp;
xue@1:         tmp=aX[0].y+aX[1].y; aX[1].y=aX[0].y-aX[1].y; aX[0].y=tmp;
xue@1:       }
xue@1:       else if (n==3) //4-point DFT
xue@1:       {
xue@1:         cdouble *aX=&X[n][l*M+M/2];
xue@1:         double tmp, tmp2;
xue@1:         tmp=aX[0].x+aX[2].x; aX[2].x=aX[0].x-aX[2].x; aX[0].x=tmp;
xue@1:         tmp=aX[0].y+aX[2].y; aX[2].y=aX[0].y-aX[2].y; aX[0].y=tmp;
xue@1:         tmp=aX[1].y+aX[3].y; tmp2=aX[1].y-aX[3].y; aX[1].y=tmp;
xue@1:         tmp=aX[3].x-aX[1].x; aX[1].x+=aX[3].x; aX[3].x=tmp2; aX[3].y=tmp;
xue@1:         tmp=aX[0].x+aX[1].x; aX[1].x=aX[0].x-aX[1].x; aX[0].x=tmp;
xue@1:         tmp=aX[0].y+aX[1].y; aX[1].y=aX[0].y-aX[1].y; aX[0].y=tmp;
xue@1:         tmp=aX[2].x+aX[3].x; aX[3].x=aX[2].x-aX[3].x; aX[2].x=tmp;
xue@1:         tmp=aX[2].y+aX[3].y; aX[3].y=aX[2].y-aX[3].y; aX[2].y=tmp;
xue@1:       }
xue@1:       else //n>3
xue@1:       {
xue@1:         cdouble *aX=&X[n][l*M+M/2];
xue@1:         for (int an=1, aL=1, aM=M/2; an<n; aL*=2, aM/=2, an++)
xue@1:         {
xue@1:           for (int al=0; al<aL; al++)
xue@1:             for (int am=0; am<aM/2; am++)
xue@1:             {
xue@1:               int iw=am*2*aL*L;
xue@1:               cdouble *lX=&aX[al*aM];
xue@1:               double x1r=lX[am].x, x1i=lX[am].y,
xue@1:                      x2r=lX[aM/2+am].x, x2i=lX[aM/2+am].y;
xue@1:               lX[am].x=x1r+x2r, lX[am].y=x1i+x2i;
xue@1:               x1r=x1r-x2r, x1i=x1i-x2i;
xue@1:               lX[aM/2+am].x=x1r*W[iw].x-x1i*W[iw].y,
xue@1:               lX[aM/2+am].y=x1r*W[iw].y+x1i*W[iw].x;
xue@1:             }
xue@1:         }
xue@1:       }
xue@1:     }
xue@1:   }
xue@1: }//CMFTC
xue@1: 
xue@1: 
xue@1: //---------------------------------------------------------------------------
xue@1: /*
xue@1:   Old versions no longer in use. For reference only.
xue@1: */
xue@1: void RFFTC_ual_old(double* Input, double *Amp, double *Arg, int Order, cdouble* W, double* XR, double* XI, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order, i, j, jj, k, *bitinv1=bitinv, Groups, ElemsPerGroup, X0, X1, X2;
xue@1:   cdouble Temp, zp, zn;
xue@1: 
xue@1:   if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1:   if (XR!=Input) for (i=0; i<N; i++) XR[i]=Input[bitinv[i]];
xue@1:   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:   N/=2;
xue@1:   double* XII=&XR[N];
xue@1:   Order--;
xue@1:   if (!bitinv1) free(bitinv);
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup=1<<i;
xue@1:     Groups=1<<(Order-i-1);
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         int kGroups=k*2*Groups;
xue@1:         Temp.x=XR[X2]*W[kGroups].x-XII[X2]*W[kGroups].y,
xue@1:         XII[X2]=XR[X2]*W[kGroups].y+XII[X2]*W[kGroups].x;
xue@1:         XR[X2]=Temp.x;
xue@1:         Temp.x=XR[X1]+XR[X2], Temp.y=XII[X1]+XII[X2];
xue@1:         XR[X2]=XR[X1]-XR[X2], XII[X2]=XII[X1]-XII[X2];
xue@1:         XR[X1]=Temp.x, XII[X1]=Temp.y;
xue@1:       }
xue@1:       X0=X0+(ElemsPerGroup<<1);
xue@1:     }
xue@1:   }
xue@1:   zp.x=XR[0]+XII[0], zn.x=XR[0]-XII[0];
xue@1:   XR[0]=zp.x;
xue@1:   XI[0]=0;
xue@1:   XR[N]=zn.x;
xue@1:   XI[N]=0;
xue@1:   for (int k=1; k<N/2; k++)
xue@1:   {
xue@1:     zp.x=XR[k]+XR[N-k], zn.x=XR[k]-XR[N-k],
xue@1:     zp.y=XII[k]+XII[N-k], zn.y=XII[k]-XII[N-k];
xue@1:     XR[k]=0.5*(zp.x+W[k].y*zn.x+W[k].x*zp.y);
xue@1:     XI[k]=0.5*(zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1:     XR[N-k]=0.5*(zp.x-W[k].y*zn.x-W[k].x*zp.y);
xue@1:     XI[N-k]=0.5*(-zn.y-W[k].x*zn.x+W[k].y*zp.y);
xue@1:   }
xue@1:   XR[N/2]=XR[N/2];
xue@1:   XI[N/2]=-XII[N/2];
xue@1:   N*=2;
xue@1: 
xue@1:   for (int k=N/2+1; k<N; k++) XR[k]=XR[N-k], XI[k]=-XI[N-k];
xue@1:   if (Amp) for (i=0; i<N; i++) Amp[i]=sqrt(XR[i]*XR[i]+XI[i]*XI[i]);
xue@1:   if (Arg) for (i=0; i<N; i++) Arg[i]=Atan2(XI[i], XR[i]);
xue@1: }//RFFTC_ual_old
xue@1: 
xue@1: void CIFFTR_old(cdouble* Input, int Order, cdouble* W, double* X, int* bitinv)
xue@1: {
xue@1:   int N=1<<Order, i, j, k, Groups, ElemsPerGroup, X0, X1, X2, *bitinv1=bitinv;
xue@1:   cdouble Temp;
xue@1:   if (!bitinv) bitinv=CreateBitInvTable(Order);
xue@1: 
xue@1:   Order--;
xue@1:   N/=2;
xue@1:   double* XII=&X[N];
xue@1: 
xue@1:   X[0]=0.5*(Input[0].x+Input[N].x);
xue@1:   XII[0]=0.5*(Input[0].x-Input[N].x);
xue@1:   for (int i=1; i<N/2; i++)
xue@1:   {
xue@1:     double frp=Input[i].x+Input[N-i].x, frn=Input[i].x-Input[N-i].x,
xue@1:            fip=Input[i].y+Input[N-i].y, fin=Input[i].y-Input[N-i].y;
xue@1:     X[i]=0.5*(frp+frn*W[i].y-fip*W[i].x);
xue@1:     XII[i]=0.5*(fin+frn*W[i].x+fip*W[i].y);
xue@1:     X[N-i]=0.5*(frp-frn*W[i].y+fip*W[i].x);
xue@1:     XII[N-i]=0.5*(-fin+frn*W[i].x+fip*W[i].y);
xue@1:   }
xue@1:   X[N/2]=Input[N/2].x;
xue@1:   XII[N/2]=-Input[N/2].y;
xue@1:   
xue@1:   ElemsPerGroup=1<<Order;
xue@1:   Groups=1;
xue@1: 
xue@1:   for (i=0; i<Order; i++)
xue@1:   {
xue@1:     ElemsPerGroup/=2;
xue@1:     X0=0;
xue@1:     for (j=0; j<Groups; j++)
xue@1:     {
xue@1:       int kGroups=bitinv[j]/2;
xue@1:       for (k=0; k<ElemsPerGroup; k++)
xue@1:       {
xue@1:         X1=X0+k;
xue@1:         X2=X1+ElemsPerGroup;
xue@1:         Temp.x=X[X2]*W[kGroups].x+XII[X2]*W[kGroups].y,
xue@1:         XII[X2]=-X[X2]*W[kGroups].y+XII[X2]*W[kGroups].x;
xue@1:         X[X2]=Temp.x;
xue@1:         Temp.x=X[X1]+X[X2], Temp.y=XII[X1]+XII[X2];
xue@1:         X[X2]=X[X1]-X[X2], XII[X2]=XII[X1]-XII[X2];
xue@1:         X[X1]=Temp.x, XII[X1]=Temp.y;
xue@1:       }
xue@1:       X0=X0+(ElemsPerGroup<<1);
xue@1:     }
xue@1:     Groups*=2;
xue@1:   }
xue@1: 
xue@1:   N*=2;
xue@1:   Order++;
xue@1:   for (i=0; i<N; i++)
xue@1:   {
xue@1:     int jj=bitinv[i];
xue@1:     if (i<jj)
xue@1:     {
xue@1:       Temp.x=X[i];
xue@1:       X[i]=X[jj];
xue@1:       X[jj]=Temp.x;
xue@1:     }
xue@1:   }
xue@1:   for (int i=0; i<N; i++) X[i]/=(N/2);
xue@1:   if (!bitinv1) free(bitinv);
xue@1: }//RFFTC_ual_old
xue@1: