FFT filter

ArdeII · 03-03-2014, 06:15 AM

Hi,
Time to give something back to the community

This is simple unoptimised version of FFT lowpass filter using WDL fft.
If you want to optimise it please send me the updated version

Code:

/*
 * FFTFilter Arto Vaarala 2014 www.kirnuarp.com
 * An example how to use Cockos WDL fft class
 * The code is not optimised to keep it as much readable as possible
 *
 * FFTFilter* fftFilter=new FFTFilter();
 * fftFilter->init(1024,GetSampleRate());
 * ...
 * fftFilter->process(input, output, nFrames);
 */
#ifndef FFTFILTERREALTYPE
#define FFTFILTERREALTYPE double
#endif
class FFTFilter{
    float cutFrequency;
    float sampleRate;
    int fftSize;
    int channels;
    int inBufferSize;
    int outBufferSize;
    FFTFILTERREALTYPE* inBuffer;
    FFTFILTERREALTYPE* outBuffer;
    FFTFILTERREALTYPE* window;
    int inBufPos;
    int outBufPos;
    WDL_FFT_COMPLEX* unsortedbuffer;
    WDL_FFT_COMPLEX*sortedbuffer;
    WDL_FFT_COMPLEX*currentFFTBuffer;
public:
    FFTFilter():inBufPos(0),outBufPos(0){}
    
    // freq 0...1
    void setCutFrequency(float freq){
        cutFrequency=freq*(sampleRate/2);
    }
    
    void reset(){
        WDL_fft_init();
        inBufPos=512;
        outBufPos=0;
        inBufferSize=0;
        outBufferSize=0;
        memset(inBuffer,0,sizeof(FFTFILTERREALTYPE)*channels*fftSize);
        memset(outBuffer,0,sizeof(FFTFILTERREALTYPE)*channels*fftSize);
    }
    
    void init(int size,float sr,int chs=1){
        cutFrequency=sr/2;
        sampleRate=sr;
        fftSize=size;
        channels=chs;
        unsortedbuffer=new WDL_FFT_COMPLEX[fftSize];
        sortedbuffer=new WDL_FFT_COMPLEX[fftSize];
        currentFFTBuffer=new WDL_FFT_COMPLEX[fftSize];

        inBuffer=new FFTFILTERREALTYPE[channels*fftSize];
        outBuffer=new FFTFILTERREALTYPE[channels*fftSize];
        reset();
        window=new FFTFILTERREALTYPE[fftSize];
        
        memset(inBuffer,0,sizeof(FFTFILTERREALTYPE)*channels*fftSize);
        memset(outBuffer,0,sizeof(FFTFILTERREALTYPE)*channels*fftSize);
        memset(window,0,sizeof(FFTFILTERREALTYPE)*fftSize);
        
        float tsc=2*M_PI/(fftSize-1);
        for (int i=0; i<fftSize; i++) {
            if(i<fftSize/2){
                window[i]=0.42323-0.49755*cos(i*tsc)+0.07922*cos(2*i*tsc);
            }
            // Make window symmetric !!!
            else{
                window[i]= 1-window[i-(fftSize/2)];
            }
        }
    }
    
    ~FFTFilter(){
        delete window;
        delete inBuffer;
        delete outBuffer;
        delete unsortedbuffer;
        delete sortedbuffer;
        delete currentFFTBuffer;
    }

    // If flush == true out buffer will be filled when in buffer reaches its end. Used when processing offline samples.
    // numSamples must be dividable by fftSize
    // All the complexity comes from poor ovarlap-add implementation :)
    int process(FFTFILTERREALTYPE* inSamples,FFTFILTERREALTYPE* outSamples,int numSamples,bool flush=false){
        int processed=0;
        int origNumSamples=numSamples;
        while(numSamples>0){
            int size=inBufPos+numSamples>fftSize/2?fftSize/2:numSamples;
            memcpy(inBuffer+inBufPos,inSamples,sizeof(FFTFILTERREALTYPE)*size);
            inSamples+=size;
            inBufPos+=size;
            numSamples-=size;

            if(inBufPos==fftSize){

                if(outBufferSize>0){
                    // move data to the beginning of the buffer
                    int outRemaining=fftSize-outBufferSize;
                    memmove(outBuffer,outBuffer+outBufferSize,sizeof(FFTFILTERREALTYPE)*outRemaining);
                    // zero the end part of the buffer
                    memset(outBuffer+outRemaining,0,sizeof(FFTFILTERREALTYPE)*(fftSize-outRemaining));
                    outBufPos=0;
                }
                
                // process fft
                fft(inBuffer,outBuffer);
                outBufferSize+=fftSize/2;
                inBufPos-=fftSize/2;
                // move data to the beginning of the buffer
                memmove(inBuffer,inBuffer+inBufPos,sizeof(FFTFILTERREALTYPE)*inBufPos);
            }
            
            if(outBufferSize>=fftSize){
                memcpy(outSamples,outBuffer,sizeof(FFTFILTERREALTYPE)*size);
                outBufferSize-=size;
                outBufPos+=size;
                outSamples+=size;
                processed+=size;
                
                // flush the remainings before finishing
                if(processed==origNumSamples-(fftSize/2)&&flush){
                    for (int i = 0; i < fftSize/2; i++,outBufPos++,inBufPos++,outSamples++){
                        *outSamples=*(outBuffer+outBufPos)+(*(inBuffer+inBufPos)*window[i]);
                    }
                    processed+=size;
                }
            }
        }
        
        return processed;
    }
    
    void fft(FFTFILTERREALTYPE *samples,FFTFILTERREALTYPE *outBuffer){
       // Bypass fft
    /*   for (int i = 0; i < fftSize; i++,outBuffer++,samples++){
            *outBuffer+=window[i]*(*samples);
        }
        return;
      */ 
        memset(sortedbuffer,0,sizeof(WDL_FFT_COMPLEX)*fftSize);
        memset(unsortedbuffer,0,sizeof(WDL_FFT_COMPLEX)*fftSize);
        memset(currentFFTBuffer,0,sizeof(WDL_FFT_COMPLEX)*fftSize);
        
        for (int i = 0; i < fftSize; i++,samples++){
            currentFFTBuffer[i].re = (WDL_FFT_REAL) *samples;
            currentFFTBuffer[i].im = 0;
        }
        
        WDL_fft(currentFFTBuffer, fftSize,false);
        
        // sort
        for (int i = 0; i < fftSize; i++){
            int j = WDL_fft_permute(fftSize, i);
            
            sortedbuffer[i].re = currentFFTBuffer[j].re;
            sortedbuffer[i].im = currentFFTBuffer[j].im;
        }
        
        // http://forum.cockos.com/showthread.php?t=56163
        /*
         Alright, I think I actually got it. From what I have gathered, the pos freqs go from 0 to N/2 where N = fftsize. That means that the pos freqs have 2 more bins than the neg freqs, which go from N/2 + 1 to N - 1.
         
         Here's an example I did of a hard high-pass, which is located where the "//do processing on sortedbuffer here" is in the previous code segment:
         */
        // Positive freqs
        for(int x = 0; x <= fftSize/2; x++){
             float F = (sampleRate * (x)) / (fftSize);
             // Lowpass filter
	     if(F > cutFrequency|| x==0){
		sortedbuffer[x].re = 0;
		sortedbuffer[x].im = 0;
            }
        }
        
		//Negative freqs
	for(int x = fftSize/2 + 1; x < fftSize; x++) {
		sortedbuffer[x].re = sortedbuffer[fftSize - x].re;
		sortedbuffer[x].im = -1 * sortedbuffer[fftSize - x].im;
        }
        
        for (int i = 0; i < fftSize; i++){
            int j = WDL_fft_permute(fftSize, i);
            unsortedbuffer[j].re = sortedbuffer[i].re;
            unsortedbuffer[j].im = sortedbuffer[i].im;
        }
        
        WDL_fft(unsortedbuffer, fftSize, true);
        
        for (int i = 0; i < fftSize; i++,outBuffer++){
            *outBuffer+=window[i]*(unsortedbuffer[i].re/fftSize);
        }
    }
};