Frequency-domain filtering with FFT

Signal filtering can be expressed in two (equivalent) approaches.

  1. convolution of the signal and impulse response of the filter - in time domain
  2. multiplication of signal spectrum and filter's transfer function (filter spectrum) - in frequency domain

This example illustrates the second approach. Signal spectrum is calculated using FFT. The filter is represented as it's transfer function - it is a pretty simple low-pass filter. The two complex arrays are multiplied and the result is plotted. Inverse FFT is used to return to a time-domain signal.

#include "aquila/global.h"
#include "aquila/transform/OouraFft.h"
#include "aquila/tools/TextPlot.h"
#include <algorithm>
#include <cmath>
#include <functional>

int main()
{
    // input signal parameters
    const std::size_t SIZE = 64;
    const Aquila::FrequencyType sampleFreq = 2000;
    const double dt = 1.0/sampleFreq;
    const Aquila::FrequencyType f1 = 96, f2 = 813;
    const Aquila::FrequencyType f_lp = 500;

    double x[SIZE];
    for (std::size_t i = 0; i < SIZE; ++i)
    {
        x[i] = 32*std::sin(2*M_PI*f1*i*dt) + 8*std::sin(2*M_PI*f2*i*dt + 0.75*M_PI);
    }
    Aquila::TextPlot plt("Signal waveform before filtration");
    plt.plot(x, SIZE);

    // calculate the FFT
    Aquila::OouraFft fft(SIZE);
    Aquila::ComplexType spectrum[SIZE];
    fft.fft(x, spectrum);
    plt.setTitle("Signal spectrum before filtration");
    plt.plotSpectrum(spectrum, SIZE);

    // generate a low-pass filter spectrum
    Aquila::ComplexType filterSpectrum[SIZE];
    for (std::size_t i = 0; i < SIZE; ++i)
    {
        if (i < (SIZE * f_lp / sampleFreq))
        {
            // passband
            filterSpectrum[i] = 1.0;
        }
        else
        {
            // stopband
            filterSpectrum[i] = 0.0;
        }
    }
    plt.setTitle("Filter spectrum");
    plt.plotSpectrum(filterSpectrum, SIZE);

    // the following line does the multiplication of two spectra
    // (which is complementary to convolution in time domain)
    std::transform(spectrum, spectrum + SIZE, filterSpectrum, spectrum,
                   std::multiplies<Aquila::ComplexType>());
    plt.setTitle("Signal spectrum after filtration");
    plt.plotSpectrum(spectrum, SIZE);

    // Inverse FFT moves us back to time domain
    double x1[SIZE];
    fft.ifft(spectrum, x1);
    plt.setTitle("Signal waveform after filtration");
    plt.plot(x1, SIZE);

    return 0;
}

The output should resemble something like this:

Signal waveform before filtration
     * *                   *                  ** *              
    *                   ***                 *                   
  *   *                      *                  *               
   *                        *                *                  
        **            *                            *            
                       *                          *             
*         *                   *           **                    
 *                             *                    *          *
                     *          *        *                      
            *       *                                        ** 
           *       *                                 **         
                                  *    *                        
                 *                      *                  *    
              *                  *   *                  *       
             * *  *                 *                  *    *   
                *                  *  *                  **     

Signal spectrum before filtration
   *                            











                          *     

    *                           
***  ********************* *****

Filter spectrum
****************                














                ****************

Signal spectrum after filtration
   *                            













    *                           
***  ***************************

Signal waveform after filtration
    ***                  ***                  ***               
   *   *                *   *                *   *              
        *                    *              *                   
  *                    *                          *             
         *                                 *                    
 *                    *       *                    *            

*         *          *         *          *         *          *

           *        *           *        *                    * 
                                                     *          
            *                           *                    *  
                   *             *                    *         
             *                         *                    *   
                  *               *                    *        
              ****                 ****                 ****
« Return to list

Comments

blog comments powered by Disqus