How to implement overlap add method?

Question

I implemented my filter, where overlap add method to prevent circular convultion is used.
input - file with noise, output should be filtered file.
My result: out is slightly modified, frequencies aren`t cut

My guess is that I wrongly multiply in the frequency domain input signal on the filter kernel (My intention is to cut off frequencies that aren't in range [300,3700]). How multiplication should be done?

I construct kernel using blackmanwindow - is my understanding correct? ( I compute amount of frequency per one sample of filter, then go through samples and see if it is in range I want to cut off I calculate frequency using formula for Blackman window.)

I just started learning DSP.

Here is my implementation (what is wrong with it???):

void DeleteFrequencies(char* fileWithNoise, char* resultFile, const int bufferSize, int lowestFrequency, int highestFrequency, int sampleRate )
{


    // |1|. files
        std::fstream in;
        std::fstream out;

in.open (fileWithNoise, std::ios::in | std::ios::binary);
out.open(resultFile, std::ios::out | std::ios::binary);

// |2|. Filter kernel design. I shall use blackman window
    // fundamental params
const int filterKernelLength = 200; // 512
const int filterMaxFrequency = sampleRate / 2; // 8000 / 2
const int frequencyPerSamle = filterMaxFrequency / filterKernelLength;
double *RealFilterResp = new double [bufferSize / 2];
double *ImmFilterResp = new double [bufferSize / 2];
    // coefficients for Blackman window
const double a0 = 0.42659;
const double a1 = 0.49656;
const double a2 = 0.076849;
    // construct filter kernel
for (int i = 0 ; i < bufferSize / 2; ++i)
{
    if ( i >= filterKernelLength ) // padd filter kernel with zeroes
    {
        RealFilterResp[i] = 0;
        ImmFilterResp[i] = 0;
    }
    else if (i * frequencyPerSamle < lowestFrequency || i * frequencyPerSamle > highestFrequency)
    {
        // apply blackman window (to eleminate frequencies < 300 hz and > 3700 hz)
        RealFilterResp[i] = a0 - a1 * cos (2 * M_PI * i / (bufferSize / 2 - 1)) + a2 * cos (4 * M_PI / (bufferSize / 2 - 1));
        ImmFilterResp[i] = a0 - a1 * cos (2 * M_PI * i / (bufferSize / 2 - 1)) + a2 * cos (4 * M_PI / (bufferSize / 2 - 1));
    }
    else
    {
        RealFilterResp[i] = 1;
        ImmFilterResp[i] = 1;
    }
}

// |3|. overlap add method
    // calculate parameters for overlap add method (we use it to prevent circular convultion)
const int FFT_length = pow (2.0 ,(int)(log(bufferSize + filterKernelLength - 1.0)/log(2.0)) + 1.0);
double *OLAP = new double[bufferSize / 2 ]; // holds the overlapping samples from segment to segment
memset(OLAP,0, bufferSize / 2 * sizeof (double));

double *RealX = new  double[bufferSize];
memset(RealX, 0, bufferSize * sizeof(double)); 
double *ImmX = new  double[bufferSize];
memset(ImmX, 0, bufferSize * sizeof(double)); 


short* audioDataBuffer = new short[bufferSize];
memset(audioDataBuffer, 0 , sizeof(short) * bufferSize);

    // start reading from file by chunks of bufferSize
while (in.good())
{
    // get proper chunk of data
    FillBufferFromFile(audioDataBuffer, bufferSize, in); // read chunk from file
    ShortArrayToDoubleArray(audioDataBuffer, RealX, bufferSize); // fill RealPart


    ForwardRealFFT(RealX, ImmX, bufferSize); // go to frequency domain

    // perform convultion as multiplication in frequency domain
    for (int j = 0; j < bufferSize / 2; ++j)
    {
        double tmp = RealX[j] * RealFilterResp[j] - ImmX[j] * ImmFilterResp[j];
        ImmX[j] = RealX[j] * ImmFilterResp[j] + ImmX[j] * RealFilterResp[j];
        RealX[j] = tmp;
    }

    // Inverse FFT
    ReverseRealFFT(RealX, ImmX, bufferSize); // go to time domain

    // add last segment overlap to this segment
    for (int j = 0; j < filterKernelLength - 2; ++j )
    {
        RealX[j] += OLAP[j];
    }

    // save samples that will overlap the next segment
    for (int j = bufferSize/2 + 1; j < bufferSize; ++j )
    {
        OLAP[j - bufferSize/2 - 1] = audioDataBuffer[j];
    }

    // write results

    DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);
    FillFileFromBuffer(audioDataBuffer, bufferSize, out);
}

/*ReverseRealFFT(RealX, ImmX, bufferSize
);
DoubleArrayToShortArray(RealX, audioDataBuffer, bufferSize);*/
delete [] audioDataBuffer;
delete [] RealFilterResp;
delete [] ImmFilterResp;
delete [] OLAP;
delete [] RealX;
delete [] ImmX;
in.close();
out.close();

}

'My guess is that I wrongly multiply ...' Can't you extinct your doubts using a debugger? — πάντα ῥεῖ
@ πάντα ῥεῖ - I think that there is logical error, which debugger coudn`t help me to spot, I have been debugging this code for 11 hours in a row and I have check it thousand times before posting to SO. — spin_eight

SleuthEye SleuthEye · Accepted Answer · 2014-06-09T15:16:00

If your intention is to use the window method to implement the filter, the window should multiply the time-domain sequence corresponding to the infinite impulse response of the ideal bandpass filter.

Specifically, for a bandpass filter of bandwidth w₀=2*pi*(3700-300)/8000 centered at w_c=2*pi*(300+3700)/8000, the ideal impulse response would be (for -infinity < n < infinity):

w0*sinc(0.5*w0*n/pi) * cos(wc*n) / pi

Which you would shift to the interval [0,N-1], and then apply the window that you computed:

double sinc(double x) {
  if (fabs(x)<1e-6) return 1.0;
  return sin(M_PI * x)/(M_PI * x);
}

void bandpassDesign(int N, double* filterImpulseResponse) {
  double w0 = 2*(3700-300)*M_PI/8000;
  double wc = 2*(300+3700)*M_PI/8000;
  double shift = 0.5*N;

  for (int i = 0; i < bufferSize; ++i) {
    double truncatedIdealResponse = w0*sinc(0.5*w0*(i-shift)/M_PI) * cos(wc*i) / M_PI;
    double window = a0 - a1 * cos (2 * M_PI * i / (N- 1)) + a2 * cos (4 * M_PI * i / (N- 1));
    filterImpulseResponse[i] = truncatedIdealResponse * window;
  }
}

You can then take the FFT to obtain the frequency-domain coefficients. Remember that if you intend on filtering data using this filter, the time sequence has to be zero padded. For example, if you wish to use a 1024-point FFT with the overlap-add method, and assuming a 128-point filter kernel meets your filter design specifications, you would call bandpassDesign with N=128, pad with 1024-128=896 zeros, then take the FFT.

How to implement overlap add method?

2 Answers