c@0: //
c@0: //  FIRFilter.cpp
c@0: //  Tempogram
c@0: //
c@0: //  Created by Carl Bussey on 25/06/2014.
c@0: //  Copyright (c) 2014 Carl Bussey. All rights reserved.
c@0: //
c@0: 
c@0: #include "FIRFilter.h"
c@7: 
c@0: using namespace std;
c@0: using Vamp::FFT;
c@0: 
c@0: FIRFilter::FIRFilter(const unsigned int lengthInput, const unsigned int numberOfCoefficients) :
c@8:     m_lengthInput(lengthInput),
c@8:     m_numberOfCoefficients(numberOfCoefficients),
c@0:     fftInput(NULL),
c@0:     fftCoefficients(NULL),
c@0:     fftReal1(NULL),
c@0:     fftImag1(NULL),
c@0:     fftReal2(NULL),
c@0:     fftImag2(NULL),
c@0:     fftFilteredReal(NULL),
c@0:     fftFilteredImag(NULL),
c@0:     fftOutputReal(NULL),
c@0:     fftOutputImag(NULL)
c@0: {
c@0:     initialise();
c@0: }
c@0: 
c@0: FIRFilter::~FIRFilter()
c@0: {
c@0:     cleanup();
c@0: }
c@0: 
c@9: //allocate memory
c@0: void
c@0: FIRFilter::initialise()
c@0: {
c@9:     //next power of 2
c@8:     m_lengthFIRFFT = pow(2,(ceil(log2(m_lengthInput+m_numberOfCoefficients-1))));
c@0:     
c@8:     fftInput = new double[m_lengthFIRFFT];
c@8:     fftCoefficients = new double[m_lengthFIRFFT];
c@8:     fftReal1 = new double[m_lengthFIRFFT];
c@8:     fftImag1 = new double[m_lengthFIRFFT];
c@8:     fftReal2 = new double[m_lengthFIRFFT];
c@8:     fftImag2 = new double[m_lengthFIRFFT];
c@8:     fftFilteredReal = new double[m_lengthFIRFFT];
c@8:     fftFilteredImag = new double[m_lengthFIRFFT];
c@8:     fftOutputReal = new double[m_lengthFIRFFT];
c@8:     fftOutputImag = new double[m_lengthFIRFFT];
c@0:     
c@8:     for(int i = 0; i < m_lengthFIRFFT; i++){
c@0:         fftInput[i] = fftCoefficients[i] = fftReal1[i] = fftImag1[i] = fftReal2[i] = fftImag2[i] = fftFilteredReal[i] = fftFilteredImag[i] = fftOutputReal[i] = fftOutputImag[i] = 0.0;
c@0:     }
c@0: }
c@0: 
c@0: void
c@0: FIRFilter::process(const float* input, const float* coefficients, float* output)
c@0: {
c@9:     //Copy to same length FFT buffers
c@8:     for(int i = 0; i < m_lengthFIRFFT; i++){
c@8:         fftInput[i] = i < m_lengthInput ? input[i] : 0.0;
c@8:         fftCoefficients[i] = i < m_numberOfCoefficients ? coefficients[i] : 0.0;
c@0:     }
c@0:     
c@8:     FFT::forward(m_lengthFIRFFT, fftInput, NULL, fftReal1, fftImag1);
c@8:     FFT::forward(m_lengthFIRFFT, fftCoefficients, NULL, fftReal2, fftImag2);
c@7:     
c@9:     //Multiply FFT coefficients. Multiplication in freq domain is convolution in time domain.
c@8:     for (int i = 0; i < m_lengthFIRFFT; i++){
c@0:         fftFilteredReal[i] = (fftReal1[i] * fftReal2[i]) - (fftImag1[i] * fftImag2[i]);
c@0:         fftFilteredImag[i] = (fftReal1[i] * fftImag2[i]) + (fftReal2[i] * fftImag1[i]);
c@0:     }
c@8:     FFT::inverse(m_lengthFIRFFT, fftFilteredReal, fftFilteredImag, fftOutputReal, fftOutputImag);
c@0:     
c@9:     //copy to output
c@8:     for (int i = 0; i < m_lengthInput; i++){
c@0:         output[i] = fftOutputReal[i];
c@0:     }
c@0: }
c@0: 
c@9: //remove memory allocations
c@0: void
c@0: FIRFilter::cleanup()
c@0: {
c@0:     delete []fftInput;
c@0:     delete []fftCoefficients;
c@0:     delete []fftReal1;
c@0:     delete []fftImag1;
c@0:     delete []fftReal2;
c@0:     delete []fftImag2;
c@0:     delete []fftFilteredReal;
c@0:     delete []fftFilteredImag;
c@0:     delete []fftOutputReal;
c@0:     delete []fftOutputImag;
c@7:     fftInput = fftCoefficients = fftReal1 = fftImag1 = fftReal2 = fftImag2 = fftFilteredReal = fftFilteredImag = fftOutputReal = fftOutputImag = NULL;
c@0: }