annotate fft.cpp @ 13:de3961f74f30 tip

Add Linux/gcc Makefile; build fix

/*
    Harmonic sinusoidal modelling and tools

    C++ code package for harmonic sinusoidal modelling and relevant signal processing.
    Centre for Digital Music, Queen Mary, University of London.
    This file copyright 2011 Wen Xue.

    This program is free software; you can redistribute it and/or
    modify it under the terms of the GNU General Public License as
    published by the Free Software Foundation; either version 2 of the
    License, or (at your option) any later version.
*/
//---------------------------------------------------------------------------

#include <string.h>
#include <stdlib.h>
#include "fft.h"

/** \file 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 Log2: Log2

  Returns x for Log2(2^x) if x is integer.
*/
int Log2(int x)
{
  return floor(log(x)/log(2)+0.5);
}//Log2

/**
  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