Multiply results of FFT class? (convolution)


#1

Hey all,

I’m trying to quickly implement convolution to test a filter design I’ve made. Here’s my excerpt of code that would hypothetically perform a convolution (notes: I’d of course pre-allocate my FFTs in a non testbed setting, and in this test case windowedFilter is already in the frequency domain and guaranteed to match the buffer in doubled size):

void FIRFilter::processFilter(AudioSampleBuffer &buffer)
{
    FFT fourierTransform(log2(buffer.getNumSamples()), false);
    FFT inverseFourierTransform(log2(buffer.getNumSamples()), true);
    buffer.setSize(buffer.getNumChannels(), buffer.getNumSamples()*2); //double size for FFT
    
    for(int channel = 0; channel < buffer.getNumChannels(); channel++)
    {
        float* bufferArray = buffer.getWritePointer(channel);
        fourierTransform.performRealOnlyForwardTransform(bufferArray);
        
        for(int n = 0; n < buffer.getNumSamples(); n++)
        {
            bufferArray[n] *= windowedFilter.getSample(0, n);
        }
        
        inverseFourierTransform.performRealOnlyInverseTransform(bufferArray);
    }
    
    buffer.setSize(buffer.getNumChannels(), buffer.getNumSamples()/2); //undo size change
}

The result of this is… nothing. The input buffer is passed out of this method exactly as it came in.

I tried changing the line within the inner for loop to be ‘bufferArray[n] *= 0’ just to see if it’d do anything- but no, the same result (input = output) occurs.

I figure there’s not an override for * in the Complex class that FFT uses. Is there any way that I can go about multiplying two buffers that were processed with FFT’s performRealOnlyForwardTransform, or no? If not, what other libraries out there might work well for my purpose (performing convolution efficiently)?


DSP: Designing a FIR Lowpass Filter / Calculating Coefficients?
#2

The problem is that setSize may delete your content. You should be using buffer.setSize (buffer.getNumChannels(), buffer.getNumSamples()/2, true). See the documentation here for more details.


#3

I’m just getting back to this now, thank you for that tip. I’ll be sure to use setSize that way, though it did not fix my problem.

Here’s the code I have currently, where the issue persists of the output signal being unaffected:

void FIRFilter::processFilter(AudioSampleBuffer &buffer)
{
    FFT fourierTransform(log2(buffer.getNumSamples()), false);
    FFT inverseFourierTransform(log2(buffer.getNumSamples()), true);
    
    float* windowArray = windowedFilter.getWritePointer(0);
    FFT::Complex* windowComplexArray = (FFT::Complex*) windowArray;
    
    buffer.setSize(buffer.getNumChannels(), buffer.getNumSamples()*2, true, true); //double size for FFT
    
    for(int channel = 0; channel < buffer.getNumChannels(); channel++)
    {
        float* bufferArray = buffer.getWritePointer(channel);
        fourierTransform.performRealOnlyForwardTransform(bufferArray);
        FFT::Complex* bufferComplexArray = (FFT::Complex*) bufferArray;
        
        for(int point=0; point < buffer.getNumSamples()/2; point++)
        {
            double newReal = (bufferComplexArray[point].r * windowComplexArray[point].r - bufferComplexArray[point].i*windowComplexArray[point].i);
            double newImag = (bufferComplexArray[point].r * windowComplexArray[point].i + bufferComplexArray[point].i * windowComplexArray[point].r);
            bufferComplexArray[point].r = newReal;
            bufferComplexArray[point].i = newImag;
        }
        
        inverseFourierTransform.performRealOnlyInverseTransform(bufferArray);
    }
    
    buffer.setSize(buffer.getNumChannels(), buffer.getNumSamples()/2, true);
}

As I said before, nothing seems to happen here, though I did confirm it reaches and runs through the inner-most for loop. setting windowedFilter to all zeros (and thusly windowArray and windowComplexArray) hypothetically should result in silence, but it does nothing (buffer before method=buffer after method).

Replacing the inner lines with

bufferComplexArray[point].r = 0;
bufferComplexArray[point].i = 0;

also does nothing.

Instead of using FFT::Complex, replacing the innermost for loop with one that just goes through all of bufferArray and replaces it with zeros also does… nothing.

It gets weirder… if I stick buffer.clear(); as the first line within this processFilter method I do get silence (as expected), but if I stick buffer.clear(); as the last line, after everything, it does nothing. Wtf?


#4

Oh… is setSize not lockfree or causing some other issues? >.>

Writing the code as follows does perform an operation on the audio:

void FIRFilter::processFilter(AudioSampleBuffer &buffer)
{
    FFT fourierTransform(log2(buffer.getNumSamples()), false);
    FFT inverseFourierTransform(log2(buffer.getNumSamples()), true);
    
    float* windowArray = windowedFilter.getWritePointer(0);
    FFT::Complex* windowComplexArray = (FFT::Complex*) windowArray;
    
    float bufferArray[buffer.getNumSamples()*2];
    
    for(int channel = 0; channel < buffer.getNumChannels(); channel++)
    {
        for(int sample=0; sample < buffer.getNumSamples(); sample++)
        {
            bufferArray[sample] = buffer.getSample(channel, sample);
        }
        
        fourierTransform.performRealOnlyForwardTransform(bufferArray);
        FFT::Complex* bufferComplexArray = (FFT::Complex*) bufferArray;
        
        for(int point=0; point < buffer.getNumSamples()/2; point++)
        {
            double newReal = (bufferComplexArray[point].r * windowComplexArray[point].r - bufferComplexArray[point].i*windowComplexArray[point].i);
            double newImag = (bufferComplexArray[point].r * windowComplexArray[point].i + bufferComplexArray[point].i * windowComplexArray[point].r);
            bufferComplexArray[point].r = newReal;
            bufferComplexArray[point].i = newImag;
        }
        
        inverseFourierTransform.performRealOnlyInverseTransform(bufferArray);
        
        for(int sample=0; sample < buffer.getNumSamples(); sample++)
        {
            buffer.setSample(channel, sample, bufferArray[sample]);
        }
    }
}

#5

Hi TristineWild,

Can you tell me what values does your windowArray has? I’m trying to multiply the values of the FFT (an then do the IFFT) but i am only able to multiply it by a constant. For example, if I want to multiply bin 10 by 0.4 and all the others by 0.3 I get no sound as a result.

Best regards,
JP


#6

Well, that’s because you can’t do filtering this way…
Convolution in spectral domain has the same constraints as in time domain. If you convolve a signal of M samples with another of N samples, you get a result of M + N samples.

When you do an FFT of M samples of you signal and then apply something different than a multiplication, you are actually doing the same as convolving in time domain a signal of N samples with an infinite impulse. But as you do that only on M samples, you end up with aliasing the infinite impulse to N samples.
And then everything does bad.

That’s why we have FIR design functions like Remez, because you can’t design a perfect filter in spectral domain and just multiply your FFT result with it (even in FFT convolution, you use zero-padding).

Yet something else I have to put in my book on audio digital signal processing because this is fundamental in the way you do filter design.


#7

FFT convolution is a very sophisticated technique, with lots of particularities and constrains.

Your approach, although logical, is a bit simplistic and overlooks many theoretical points.

There are many pappers dedicated to the subject, I particularly recommend this chapter:
http://www.dspguide.com/ch18/2.htm

It doesn’t explore all the subject, but will give you a better picture.