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@0: #include <cmath>
c@0: #include <vamp-sdk/FFT.h>
c@0: #include <assert.h>
c@0: #include <iostream>
c@0: using namespace std;
c@0: 
c@0: using Vamp::FFT;
c@0: 
c@0: FIRFilter::FIRFilter(const unsigned int lengthInput, const unsigned int numberOfCoefficients) :
c@0:     _lengthInput(lengthInput),
c@0:     _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@0: void
c@0: FIRFilter::initialise()
c@0: {
c@0:     _lengthFIRFFT = pow(2,(ceil(log2(_lengthInput+_numberOfCoefficients-1))));
c@0:     
c@0:     fftInput = new double[_lengthFIRFFT];
c@0:     fftCoefficients = new double[_lengthFIRFFT];
c@0:     fftReal1 = new double[_lengthFIRFFT];
c@0:     fftImag1 = new double[_lengthFIRFFT];
c@0:     fftReal2 = new double[_lengthFIRFFT];
c@0:     fftImag2 = new double[_lengthFIRFFT];
c@0:     fftFilteredReal = new double[_lengthFIRFFT];
c@0:     fftFilteredImag = new double[_lengthFIRFFT];
c@0:     fftOutputReal = new double[_lengthFIRFFT];
c@0:     fftOutputImag = new double[_lengthFIRFFT];
c@0:     
c@0:     for(int i = 0; i < _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@0:     float max = 0;
c@0:     for(int i = 0; i < _lengthInput; i++){
c@0:         fftInput[i] = input[i];
c@0:         max = max > fftInput[i] ? max : fftInput[i];
c@0:         //cout << fftInput[i] << endl;
c@0:     }
c@0:     //cout << max << endl;
c@0:     for(int i = 0; i < _numberOfCoefficients; i++){
c@0:         fftCoefficients[i] = coefficients[i];
c@0:         //cout << fftCoefficients[i] << endl;
c@0:     }
c@0:     
c@0:     FFT::forward(_lengthFIRFFT, fftInput, NULL, fftReal1, fftImag1);
c@0:     FFT::forward(_lengthFIRFFT, fftCoefficients, NULL, fftReal2, fftImag2);
c@0:     for (int i = 0; i < _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@0:     FFT::inverse(_lengthFIRFFT, fftFilteredReal, fftFilteredImag, fftOutputReal, fftOutputImag);
c@0:     
c@0:     max = 0;
c@0:     for(int i = 0; i < _lengthInput; i++){
c@0:         output[i] = fftOutputReal[i];
c@0:         max = max > output[i] ? max : output[i];
c@0:     }
c@0: }
c@0: 
c@0: void
c@0: FIRFilter::cleanup()
c@0: {
c@0:     delete []fftInput;
c@0:     fftInput = NULL;
c@0:     delete []fftCoefficients;
c@0:     fftCoefficients = NULL;
c@0:     delete []fftReal1;
c@0:     fftReal1 = NULL;
c@0:     delete []fftImag1;
c@0:     fftImag1 = NULL;
c@0:     delete []fftReal2;
c@0:     fftReal2 = NULL;
c@0:     delete []fftImag2;
c@0:     fftImag2 = NULL;
c@0:     delete []fftFilteredReal;
c@0:     fftFilteredReal = NULL;
c@0:     delete []fftFilteredImag;
c@0:     fftFilteredImag = NULL;
c@0:     delete []fftOutputReal;
c@0:     fftOutputReal = NULL;
c@0:     delete []fftOutputImag;
c@0:     fftOutputImag = NULL;
c@0: }