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: #ifndef fftH
xue@1: #define fftH
xue@1: 
Chris@5: /**
Chris@5:   \file fft.h - fast Fourier and cosine transforms
xue@1: 
xue@1:   Arguments of the function in this unit roughly follow the following rules:
xue@1:     Wid   : size of transform
xue@1:     Order : log2(Wid)
xue@1:     W     : FFT buffer: twiddle factors, standard size Wid/2
xue@1:     X     : FFT buffer: data buffer, standard size Wid
xue@1:     bitinv: FFT buffer: table of bit-inversed integers, standard size Wid
xue@1:     Window: window function, standard size Wid
xue@1:     spec  : complex spectrogram
xue@1:     Amp   : amplitude spectrogram
xue@1:     Arg   : phase spectrogram
xue@1:     prefix h- : half-size
xue@1:     prefix q- : quarter-size
xue@1: */
xue@1: 
xue@1: 
xue@1: #include <math.h>
xue@1: #include "xcomplex.h"
xue@1: #include "align8.h"
xue@1: 
xue@1: //---------------------------------------------------------------------------
Chris@5: /**
xue@1:   macro AllocateFFTBuffer: allocates FFT buffers X and W of given length with an extra LDATA buffer the
xue@1:   right size for hosting such arrays as the standard-size window function or output amplitude spectrum.
xue@1: 
xue@1:   This macro allocates buffer with 64-bit alignment, which must be freed with free8() or FreeFFTBuffer
xue@1:   macro.
xue@1: */
xue@1: #define AllocateFFTBuffer(_WID, _LDATA, _W, _X) \
xue@1:   double *_LDATA=(double*)malloc8(sizeof(double)*_WID*4); \
xue@1:   cdouble *_X=(cdouble*)&_LDATA[_WID]; \
xue@1:   cdouble *_W=(cdouble*)&_X[_WID]; \
xue@1:   SetTwiddleFactors(_WID, _W);
xue@1: #define FreeFFTBuffer(_LDATA) free8(_LDATA);
xue@1: 
xue@1: //--compelx FFT and IFFT-----------------------------------------------------
xue@1: void CFFTCbii(int Order, cdouble* W, cdouble* X); //complex FFT subroutine after bit-inversed ordering
xue@1: void CFFTC(int Order, cdouble* W, cdouble* X, int* bitinv=0); //complex FFT, compact in-place
xue@1: void CFFTC(cdouble* Input, double* Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv=0); //complex FFT
xue@1: void CFFTCW(double* Window, int Order, cdouble* W, cdouble* X, int* bitinv=0);  //complex FFT with window, compact in-place
xue@1: void CFFTCW(cdouble* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* bitinv=0);  //complex FFT with window
xue@1: void CIFFTC(int Order, cdouble* W, cdouble* X, int* bitinv=0);  //complex IFFT, compact in-place
xue@1: void CIFFTC(cdouble* Input, int Order, cdouble* W, cdouble* X, int* bitinv=0);  //complex IFFT
xue@1: 
xue@1: //--real FFT and IFFT--------------------------------------------------------
xue@1: void RFFTC_ual(double* Input, double* Amp, double* Arg, int Order, cdouble* W, cdouble* X, int* hbitinv=0); //real-to-complex FFT
xue@1: void RFFTC_ual_old(double* Input, double *Amp, double *Arg, int Order, cdouble* W, double* XR, double* XI, int* bitinv=0);  //depreciated
xue@1: void RFFTCW(double* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv=0);  //real-to-complex FFT with window, double-type input
xue@1: void RFFTCW(__int16* Input, double* Window, double *Amp, double *Arg, int Order, cdouble* W, cdouble* X, int* hbitinv=0); //real-to-complex FFT with window, __int16 input
xue@1: void CIFFTR(cdouble* Input, int Order, cdouble* W, double* X, int* hbitinv=0);  //complex-to-real IFFT
xue@1: 
xue@1: /*
xue@1:   macro RFFTC: was defined as RFFTC_ual (literally "unaligned RFFTC") enveloped by ALIGN8(). However, as
xue@1:   ALIGN8() is currently disabled for lack of compiler support, RFFTC is equivalent to RFFTC_ual.
xue@1: */
xue@1: #define RFFTC(Input, Amp, Arg, Order, W, X, hbitinv) RFFTC_ual(Input, Amp, Arg, Order, W, X, hbitinv);
xue@1: 
xue@1: //--MFT----------------------------------------------------------------------
xue@1: void CMFTC(double* xR, double* xI, int Order, cdouble** X, cdouble* W); //complex fast multiresolution FT
xue@1: 
xue@1: //--DCT and IDCT-------------------------------------------------------------
xue@1: void RDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv=0);   //DCT1
xue@1: void RDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv=0);   //DCT4
xue@1: void RIDCT1(double* Input, double* Output, int Order, cdouble* qW, cdouble* qX, int* qbitinv=0);  //IDCT1
xue@1: void RIDCT4(double* Input, double* Output, int Order, cdouble* hW, cdouble* hX, int* hbitinv=0);  //IDCT4
xue@1: 
xue@1: //--LCT----------------------------------------------------------------------
xue@1: void RLCT(double** spec, double* data, int Count, int Order, int wid, double* g); //local cosine transform
xue@1: void RILCT(double* data, double** spec, int Fr, int Order, int wid, double* g);   //inverse local cosine transform
xue@1: 
xue@1: //--tools--------------------------------------------------------------------
xue@1: double Atan2(double, double);
xue@11: int Log2(int);
xue@1: int* CreateBitInvTable(int Order);  //creates table of bit-inversed integers
xue@1: void SetTwiddleFactors(int N, cdouble* w);  //set twiddle factors
xue@1: 
xue@1: 
xue@1: #endif