diff fft.cpp @ 1:6422640a802f

first upload

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