/Users/ljp/code/gamma/trunk/Gamma/DFT.h

#ifndef GAMMA_DFT_H_INC
#define GAMMA_DFT_H_INC

/*  Gamma - Generic processing library
    See COPYRIGHT file for authors and license information */

#include <math.h>
#include "Gamma/mem.h"
#include "Gamma/scl.h"
#include "Gamma/tbl.h"
#include "Gamma/Sync.h"
#include "Gamma/Constants.h"
#include "Gamma/Containers.h"
#include "Gamma/FFT.h"
#include "Gamma/Types.h"


namespace gam{


enum SpectralType{
    COMPLEX,    
    MAG_PHASE,  
    MAG_FREQ    
};



template <class T=gam::real>
class SlidingWindow{
public:
    SlidingWindow(uint32_t winSize, uint32_t hopSize);
    ~SlidingWindow();

    void resize(uint32_t winSize, uint32_t hopSize);
    void sizeHop(uint32_t size);
    void sizeWin(uint32_t size);
    
    uint32_t sizeHop() const;
    uint32_t sizeWin() const;
    
    
    const T * window();
    const T * operator()();
    
    bool operator()(T input);
    
    
    bool operator()(T * dst, T input);

protected:
    void slide();   // Slides samples in window left by hop num samples
    // After processing the window samples, a call to slide() must be made
    // to prepare for the next window.

protected:
    T * mBuf;
    uint32_t mSizeWin, mSizeHop;
    uint32_t mTapW; // current index to write to
    uint32_t mHopCnt;   // counts samples for hop
    
private:
    uint32_t hopStart() const;
};



template <class T=gam::real>
class DFTBase : public Synced{
public:
    DFTBase();
    virtual ~DFTBase();

    T * aux(uint32_t num);
    
    Complex<T> * bins() const { return mBins; }
    
    Complex<T>& bin(uint32_t k){ return mBins[k]; }
    
    const Complex<T>& bin(uint32_t k) const { return mBins[k]; }
    
    double binFreq() const;     
    uint32_t numBins() const;   
    uint32_t sizeDFT() const;   
    Sync& syncFreq();           
    
    void numAux(uint32_t num);  

    virtual void onResync(double r);

protected:
    uint32_t mSizeDFT, mNumAux;
    union{
        T * mBuf;       // FFT buffer
        Complex<T> * mBins;
    };
    T * mAux;       // aux buffers
    Sync mSyncFreq;

    T normForward() const;  // get norm factor for forward transform values
    T * bufPos(){ return mBuf+1; }
    T * bufFrq(){ return mBuf; }
};



class DFT : public DFTBase<float>{
public:

    DFT(
        uint32_t winSize=1024, uint32_t padSize=0,
        SpectralType specType=COMPLEX,
        uint32_t numAux=0
    );

    virtual ~DFT();


    DFT& spectrumType(SpectralType v);

    DFT& precise(bool whether);

    void resize(uint32_t windowSize, uint32_t padSize);

    float freqRes() const;              
    float overlap() const;              
    bool overlapping() const;           
    uint32_t sizeHop() const;           
    uint32_t sizePad() const;           
    uint32_t sizeWin() const;           

    Sync& syncHop();                    

    
    bool operator()(float input);

    
    float operator()();

    
    void forward(const float * src);    

    
    virtual void inverse(float * dst);
    
    
    bool inverseOnNext();

    void spctToRect();      // convert spectrum to rectangle format
    void spctToPolar();     // convert spectrum to polar format
    void zero();            
    void zeroEnds();        

    virtual void onResync(double r);
    virtual void print(FILE * fp=stdout, const char * append="\n");
    
protected:
    uint32_t mSizeWin;              // samples in analysis window
    uint32_t mSizeHop;              // samples between forward transforms (= winSize() for DFT)
    SpectralType mSpctFormat;       // format of spectrum
    RFFT<float> mFFT;
    Sync mSyncHop;
    
    // Buffers
    float * mPadOA;         // Overlap-add buffer (alloc'ed only if zero-padded)
    float * mBufInv;        // Pointer to inverse sample buffer
    uint32_t mTapW, mTapR;  // DFT i/o read/write taps
    bool mPrecise;
};





class STFT : public DFT {
public:

    STFT(uint32_t winSize=1024, uint32_t hopSize=256, uint32_t padSize=0,
        WindowType winType = RECTANGLE,
        SpectralType specType = COMPLEX,
        uint32_t numAux=0
    );
    
    virtual ~STFT();


    using DFT::operator();
    using DFT::sizeHop;


    
    bool operator()(float input);

    void forward(const float * src);
    
    virtual void inverse(float * dst);


    void resize(uint32_t winSize, uint32_t padSize);

    STFT& inverseWindowing(bool v){
        mWindowInverse=v; computeInvWinMul(); return *this; }

    STFT& rotateForward(bool v){ mRotateForward=v; return *this; }

    STFT& sizeHop(uint32_t size);
    
    STFT& windowType(WindowType type);


    float unitsHop();
    
    float * phases();
    
    virtual void print(FILE * fp=stdout, const char * append="\n"); 

protected:
    void computeInvWinMul();    // compute inverse normalization factor (due to overlap-add)
    
    SlidingWindow<float> mSlide;
    float * mFwdWin;            // forward transform window
    float * mPhases;            // copy of current phases (mag-freq mode)
    WindowType mWinType;        // type of analysis window used
    float mFwdWinMul, mInvWinMul;
    bool mWindowInverse;        // whether to window inverse samples
    bool mRotateForward;
};





template <class T>
class SDFT : public DFTBase<T> {
public:

    SDFT(uint32_t sizeDFT, uint32_t binLo, uint32_t binHi);
    
    void forward(T input);
    
    SDFT& interval(uint32_t binLo, uint32_t binHi);

    void resize(uint32_t sizeDFT, uint32_t binLo, uint32_t binHi);
        
protected:
    uint32_t mBinLo, mBinHi;
    DelayN<T> mDelay;
    //T * mBufI;                // alias to imaginaries
    Complex<T> mF1;
    Complex<T> mFL;
    //double mC0, mS0;      // phasor rotators
    //double mCL, mSL;      // low bin phasor
    T mNorm;                // fwd transform normalization
};



// Implementation_______________________________________________________________

//---- SlidingWindow
#define TEM template<class T>
TEM SlidingWindow<T>::SlidingWindow(uint32_t winSize, uint32_t hopSize)
    : mBuf(0), mSizeWin(0), mSizeHop(0), mHopCnt(0)
{
    resize(winSize, hopSize);
    mem::deepZero(mBuf, sizeWin());
}

TEM SlidingWindow<T>::~SlidingWindow(){
    if(mBuf){ free(mBuf); mBuf = 0; }
}

TEM void SlidingWindow<T>::resize(uint32_t winSize, uint32_t hopSize){
    sizeWin(winSize);
    sizeHop(hopSize);
    //mTapW = hopStart();   // for single-buffer slide mode
    mTapW = 0;              // for single-buffer rotate mode
}

TEM void SlidingWindow<T>::sizeWin(uint32_t size){
    if(mem::resize(mBuf, sizeWin(), size)) mSizeWin = size;
}

TEM void SlidingWindow<T>::sizeHop(uint32_t size){
    mSizeHop = scl::clip(size, sizeWin(), (uint32_t)1);
}

TEM inline uint32_t SlidingWindow<T>::sizeHop() const { return mSizeHop; }
TEM inline uint32_t SlidingWindow<T>::sizeWin() const { return mSizeWin; }
TEM inline const T * SlidingWindow<T>::window(){ return mBuf; }
TEM inline const T * SlidingWindow<T>::operator()(){ return mBuf; }

TEM inline bool SlidingWindow<T>::operator()(T input){

// faster, but requires explicit call to slide() by user which is error prone
//  mBuf[mTapW] = input;
//  
//  if(++mTapW == sizeWin()){
//      mTapW = hopStart();
//      return true;
//  }
//  return false;
    
    mBuf[mTapW] = input;
    
    if(++mTapW >= sizeHop()){
        mTapW = 0;
        mem::rotateLeft(sizeHop(), mBuf, sizeWin());
        return true;
    }
    return false;
}

TEM inline bool SlidingWindow<T>::operator()(T * output, T input){
    mBuf[mTapW] = input;
    if(++mTapW == sizeWin()) mTapW = 0;

    if(++mHopCnt == sizeHop()){
        mem::copyAllFromRing(mBuf, sizeWin(), mTapW, output);
        mHopCnt = 0;
        return true;
    }
    return false;
}

TEM inline void SlidingWindow<T>::slide(){
    mem::deepMove(mBuf, mBuf + sizeHop(), hopStart());
}

TEM inline uint32_t SlidingWindow<T>::hopStart() const { return sizeWin() - sizeHop(); }



//---- DFTBase

TEM DFTBase<T>::DFTBase() : mSizeDFT(0), mNumAux(0), mBuf(0), mAux(0){ initSynced(); }

TEM DFTBase<T>::~DFTBase(){ //printf("~DFTBase\n");
    mem::free(mBuf);
    mem::free(mAux);
}

TEM inline T *      DFTBase<T>::aux(uint32_t num){ return mAux + numBins() * num; }
TEM inline double   DFTBase<T>::binFreq() const { return spu() / (double)sizeDFT(); }
TEM inline uint32_t DFTBase<T>::numBins() const { return (sizeDFT() + 2)>>1; }
TEM inline uint32_t DFTBase<T>::sizeDFT() const { return mSizeDFT; }
TEM inline Sync&    DFTBase<T>::syncFreq(){ return mSyncFreq; }
TEM inline T        DFTBase<T>::normForward() const { return (T)2 / (T)sizeDFT(); }

TEM void DFTBase<T>::numAux(uint32_t num){
    if(mem::resize(mAux, mNumAux * numBins(), num * numBins())) mNumAux = num;
}

TEM void DFTBase<T>::onResync(double r){
    mSyncFreq.ups(binFreq());
}



//---- DFT

inline DFT& DFT::spectrumType(SpectralType v){ mSpctFormat=v; return *this; }
inline DFT& DFT::precise(bool w){ mPrecise=w; return *this; }

inline float DFT::freqRes() const { return spu() / sizeWin(); }
inline float DFT::overlap() const { return float(sizeWin()) / sizeHop(); }
inline bool DFT::overlapping() const { return sizeHop() < sizeWin(); }
inline uint32_t DFT::sizeHop() const { return mSizeHop; }
inline uint32_t DFT::sizePad() const { return mSizeDFT - mSizeWin; }
inline uint32_t DFT::sizeWin() const { return mSizeWin; }
inline Sync& DFT::syncHop(){ return mSyncHop; }

inline bool DFT::operator()(float input){
    bufPos()[mTapW] = input;

    if(++mTapW >= sizeHop()){
        forward(bufPos());
        mTapW = 0;
        return true;
    }
    return false;
}

inline float DFT::operator()(){
    if(++mTapR >= sizeHop()){
        inverse(0); // this is a virtual method
        mTapR = 0;
    }
    return mBufInv[mTapR];
}

inline bool DFT::inverseOnNext(){ return mTapR == (sizeHop() - 1); }
inline void DFT::zero(){ mem::deepZero(mBuf, sizeDFT() + 2); }
inline void DFT::zeroEnds(){ 
    //bins0()[0]=0; bins0()[numBins()-1]=0; 
    //bins1()[0]=0; bins1()[numBins()-1]=0;
    bins()[0](0,0);
    bins()[numBins()-1](0,0);
}


//---- STFT

inline bool STFT::operator()(float input){
    if(mSlide(bufPos(), input)){
        forward(bufPos());
        return true;
    }
    return false;
}

inline float STFT::unitsHop(){ return (float)DFT::sizeHop() * ups(); }
inline float * STFT::phases(){ return mPhases; }



//---- SDFT

TEM SDFT<T>::SDFT(uint32_t sizeDFT, uint32_t binLo, uint32_t binHi)
    : DFTBase<T>(), mBinLo(0), mBinHi(0), mDelay(0)
{
    resize(sizeDFT, binLo, binHi);
}

TEM void SDFT<T>::resize(uint32_t sizeDFT, uint32_t binLo, uint32_t binHi){
    // may be able to keep these smaller?
    mem::resize(this->mBuf, this->mSizeDFT + 2, sizeDFT + 2);
    mem::deepZero(this->mBuf, sizeDFT + 2);

    mDelay.resize(sizeDFT);
    mDelay.assign(T(0));

    this->mSizeDFT = sizeDFT;

    interval(binLo, binHi);

    //this->onSyncChange();
}

TEM SDFT<T>& SDFT<T>::interval(uint32_t binLo, uint32_t binHi){
    mBinLo = binLo;
    mBinHi = binHi;
    
    double theta = M_2PI / this->sizeDFT();

    mF1.fromPhase(theta);
    mFL.fromPhase(theta*mBinLo);
    mNorm = T(2) / T(this->sizeDFT());
    return *this;
}

//
TEM inline void SDFT<T>::forward(T input){
    T dif = (input - mDelay(input)) * mNorm;    // ffd comb zeroes
                                                // difference between temporal 'frames'
    Complex<T> c = mFL;                         // phasor at low bin

    // apply complex resonators:
    // multiply freq samples by 1st harmonic (shift time signal)
    // add time sample to all bins (set time sample at n=0)
    
    for(uint32_t k=mBinLo; k<mBinHi; ++k){
        Complex<T>& b = this->bins(k);
        b = b*c + dif;
        c *= mF1;
    }
}

//TEM inline T SDFT<T>::inverse(){}



/*
TEM class Counter{
public:

    bool operator()(){
        if(++mVal == mMax){
            onWrap();
            mVal = (T)0;
            return true;
        }
        return false;
    }

    virtual void onWrap();

protected:
    T mVal, mMax;
};
*/

#undef TEM

} // gam::

#endif