C / C++ FFT (Fast Fourier Transform) Based Signal Filter

The following design is a FFT (Fast Fourier Transform) based signal filter developed in C / C++.
You can use this type of filter to amplify or dampen very specific bands.
Alternatively you could use it as a band pass, low pass, or high pass filter by simply setting coefficient ranges to zero.
And yet another great application is to use it as a pitch shifter, for example, with digital audio, and keep the timing the same. I.e. change the pitch without playing the sound faster or slower.

//**************************** class DSP_FFT_Filter ***************************
// Defines are filter based on Fast Fourier Transform
// transferFunctionCoef and numberOfElements define the transfer function.
// these coefficients are multiplied with the complex vectors returned by FFT
// Note: coeff vectors appear in reverse order, ie. 
// transferFunctionCoef[numberOfElements-1] is the base frequency,
// transferFunctionCoef[numberOfElements-2] is first harmonic.
class DSP_FFT_Filter : public DSP_ProcessorBase<double, double>
{
public:
    DSP_FFT_Filter(float *transferFunctionCoef, int numberOfElements, int blockSize);
    virtual ~DSP_FFT_Filter()
    {
        if (m_pTransferFunctionCoef) delete[] m_pTransferFunctionCoef;
#ifdef _DEBUG
        if (m_pFFTCoefficients) delete[] m_pFFTCoefficients;
#endif
    }
    virtual char* GetProcessorName(void){return "DSP_FFT_Filter";}
    virtual bool ProcessSamples(const DSP_Samples &inputSamples);
    virtual bool ProcessFFTCoeff(float *data, unsigned long nn);
                                                                 
    float *m_pTransferFunctionCoef;
    int m_NumberOfElementsInTransferFunction;

#ifdef _DEBUG
    float *m_pFFTCoefficients;
    int m_NumberOfCoeffients;
#endif
};

The header above is simple, implementation follows below:

//************************ class DSP_FFT_Filter *******************************

// returns true if number is a power of 2
bool IsPowerOfTwo(unsigned number)
{
    bool result = false;
    unsigned mask = 1;
    for (int i = 0; i < 32; i++)
    {
        if ((number & mask) == number) // if one bit is set only, return true
        {
            result = true;
            break;
        }
        else
            mask <<= 1; // shift mask and try again...
    }
    return result;
}

// returns bit number of MSB, ex: number = 8 -> retVal = 3
int GetHighestBit(unsigned number)
{
    unsigned mask = 0x80000000;
    for (int i = 31; i >= 0; i--)
    {
        if (number & mask) 
            return i;
        else 
            mask >>= 1;
    }
    return 0;
}

// bool four1(float data[], unsigned long nn, int isign);
// Replaces data[1..2*nn] by its discrete Fourier transform, if isign is input as 1; or replaces
// data[1..2*nn] by nn times its inverse discrete Fourier transform, if isign is input as - 1.
// data is a complex array of length nn or, equivalently, a real array of length 2*nn. nn MUST
// be an integer power of 2 (this is not checked for!).
 //       data[0] is NOT used
 //       data[1] contains the real part of element #1
 //       data[2] contains the imaginary part of element #1
 //       and so on ...

// in case of real data, nn MUST BE A POWER OF TWO!!!
// returns true if check was ok
#define SWAP(a,b) tempr=(a);(a)=(b);(b)=tempr
bool four1(float data[], unsigned long nn, int isign)
{
    bool result = false;
    assert(IsPowerOfTwo(nn));
    IF_NLOG(IsPowerOfTwo(nn), DSPLIB_ET_FFT_INPUT_MUST_BE_POWER_OF_TWO, "four1")
    {
        unsigned long n, mmax, m, j, istep, i;
        double wtemp, wr, wpr, wpi, wi, theta; 
        float tempr, tempi;
    
        n = nn << 1;
        j = 1;
        for (i = 1; i < n; i += 2) 
        {                               // This is the bit-reversal section of the routine. 
            if (j > i) 
            {
                SWAP(data[j],data[i]); //Exchange the two complex numbers.
                SWAP(data[j+1],data[i+1]);
            }
            m=nn;
            while (m >= 2 && j > m) 
            {
                j -=m;
                m >>= 1;
            }
            j +=m;
        }
        //Here begins the Danielson-Lanczos section of the routine.
                                     
        mmax = 2;
        while (n > mmax) 
        {                   // Outer loop executed log 2 nn times.
            istep = mmax << 1;
            theta = isign * (6.28318530717959/mmax); //Initialize the trigonometric recurrence.
            wtemp = sin (0.5 * theta);
            wpr = -2.0 * wtemp * wtemp;

            wpi=sin(theta);
            wr=1.0;
            wi=0.0;
            for (m=1;m<mmax;m+=2)
            {                           //Here are the two nested inner loops.
                for (i=m;i<=n;i+=istep) 
                {
                    j = i + mmax; //This is the Danielson-Lanczos for-mula:
                    tempr = float(wr * double(data[j]) - wi * double(data[j+1]));
                    tempi = float(wr * double(data[j+1]) + wi * double(data[j]));
                    data[j] = data[i] - tempr;
                    data[j+1] = data[i+1] - tempi;
                    data[i] += tempr;
                    data[i+1] += tempi;
                }
                wr = (wtemp = wr) * wpr - wi * wpi + wr;// Trigonometric recurrence.
                wi = wi * wpr + wtemp * wpi + wi;
            }
            mmax=istep;
        }
        result = true;
    }
    return result;
}
              
DSP_FFT_Filter::DSP_FFT_Filter(float *transferFunctionCoef, int numberOfElements, int blockSize)
        : DSP_ProcessorBase<double, double>(blockSize)
{
    m_pTransferFunctionCoef = NULL;
    m_NumberOfElementsInTransferFunction = 0;
    
    // copy over transferFunction coefficients
    IF_NLOG(transferFunctionCoef, DSPLIB_ET_FFT_FILTER_TRANSFER_FUNC_COEFS_MISSING, "")
    {
        m_pTransferFunctionCoef = new float[numberOfElements];
        m_NumberOfElementsInTransferFunction = numberOfElements;
        for (int i = 0; i < numberOfElements; i++)
        {
            m_pTransferFunctionCoef[i] = transferFunctionCoef[i];
        }
    }
#ifdef _DEBUG
    m_pFFTCoefficients = NULL;
    m_NumberOfCoeffients = 0;
#endif

}


bool DSP_FFT_Filter::ProcessSamples(const DSP_Samples &inputSamples)
{
    bool result = false;

    // ensure input is power of two...if not we will need to add some extra
    // padding
    // so figure out if input length is power of two, if not compute extra space needed
    int extraSpaceRequired = 0;
    if (IsPowerOfTwo(inputSamples.m_NumberOfSamples) == false)
    {
        extraSpaceRequired = (1 << (GetHighestBit(inputSamples.m_NumberOfSamples) + 1)) - inputSamples.m_NumberOfSamples;
    }

    DSP_Samples *outputSamples = (DSP_Samples*)m_pOutputSamples;
    IF_NLOG(outputSamples, DSPLIB_ET_OUTPUTSAMPLE_OBJ_MISSING, "")
    {
        outputSamples->ReallocWithoutCopy(inputSamples.m_NumberOfSamples + extraSpaceRequired);
        assert(outputSamples->m_BlockSize >= inputSamples.m_NumberOfSamples + extraSpaceRequired);
        IF_NLOG(outputSamples->m_BlockSize >= inputSamples.m_NumberOfSamples + extraSpaceRequired, DSPLIB_ET_OUTPUTSAMPLE_OBJ_BLOCKSIZE_TOO_SMALL, "")
        {
            outputSamples->m_NumberOfSamples = inputSamples.m_NumberOfSamples + extraSpaceRequired;    
            outputSamples->m_SampleRate = inputSamples.m_SampleRate;
                                                                           
            assert(IsPowerOfTwo(inputSamples.m_NumberOfSamples + extraSpaceRequired));
            IF_NLOG(IsPowerOfTwo(inputSamples.m_NumberOfSamples + extraSpaceRequired), DSPLIB_ET_FFT_INPUT_MUST_BE_POWER_OF_TWO, "")
            {
                int originalNumberOfSamples = inputSamples.m_NumberOfSamples + extraSpaceRequired;
                int numberOfSamples = inputSamples.m_NumberOfSamples;
    
                // allocate temp buffer
                float *tmpData = new float[(numberOfSamples + extraSpaceRequired) * 2 + 2];
                memset(tmpData, 0, sizeof(float) * ((numberOfSamples + extraSpaceRequired) * 2 + 2));

                IF_NLOG(tmpData, ET_OUT_OF_MEMORY, "")
                {
                    // convert samples into floats first...
                    for (int i = 0; i < (numberOfSamples + extraSpaceRequired) && i < inputSamples.m_NumberOfSamples; i++) 
                    {               
                        tmpData[1 + i * 2] = float(inputSamples.m_pSamples[i]);
                    }

                    // run FFT
                    four1(tmpData, numberOfSamples + extraSpaceRequired, 1);

        #ifdef _DEBUG 
                    // if in debug mode, store FFT results for later review
                    if (m_pFFTCoefficients) delete[] m_pFFTCoefficients;
                    m_NumberOfCoeffients = numberOfSamples + extraSpaceRequired;
                    m_pFFTCoefficients = new float[(numberOfSamples + extraSpaceRequired) * 2];
                    for(int jj = 0; jj < numberOfSamples + extraSpaceRequired; jj++)
                    {
                        m_pFFTCoefficients[jj * 2] = tmpData[1 + jj*2];
                        m_pFFTCoefficients[jj * 2 + 1] = tmpData[1 + jj*2 + 1];
                    }
        #endif
                    IF_NLOG(ProcessFFTCoeff(tmpData, numberOfSamples + extraSpaceRequired), DSPLIB_ET_FFT_FILTER_PROCESS_COEFS_FAILED, "")
                    {           
                        // run IFFT
                        four1(tmpData, numberOfSamples + extraSpaceRequired, -1);

                        // convert samples back to doubles and normalize...and we're done!
                        for (i = 0; i < numberOfSamples + extraSpaceRequired; i++)
                        {
                            outputSamples->m_pSamples[i] = double(tmpData[1+i*2]) / double(originalNumberOfSamples);
                        }
                        outputSamples->m_NumberOfSamples = originalNumberOfSamples;
                        result = true;
                    }
                    delete[] tmpData;
                }    
            }
        }
    }
    return result;
}

bool DSP_FFT_Filter::ProcessFFTCoeff(float *tmpData, unsigned long nn)
{
    bool result = true;

    // apply transfer function coef by MULTIPLYING (modify filter attenuation)
    for (int i = 0; i < m_NumberOfElementsInTransferFunction && i < int(nn); i++)
    {
        // and write back result
        tmpData[1 + i*2] = tmpData[1 + i*2] * m_pTransferFunctionCoef[i];
        tmpData[1 + i*2 + 1] = tmpData[1 + i*2 + 1] * m_pTransferFunctionCoef[i];
    }
    return result;
}