kfr

Fast, modern C++ DSP framework, FFT, Sample Rate Conversion, FIR/IIR/Biquad Filters (SSE, AVX, AVX-512, ARM NEON)
Log | Files | Refs | README

fir.cpp (6933B)


      1 /**
      2  * KFR (https://www.kfrlib.com)
      3  * Copyright (C) 2016-2023 Dan Cazarin
      4  * See LICENSE.txt for details
      5  */
      6 
      7 #include <kfr/base.hpp>
      8 #ifdef HAVE_DFT
      9 #include <kfr/dft.hpp>
     10 #endif
     11 #include <kfr/dsp.hpp>
     12 #include <kfr/io.hpp>
     13 
     14 using namespace kfr;
     15 
     16 // Define a macro to check if Python is installed
     17 #ifndef PYTHON_IS_INSTALLED
     18 #define PYTHON_IS_INSTALLED 1
     19 #endif
     20 
     21 int main()
     22 {
     23     // Print the version of the KFR library being used
     24     println(library_version());
     25 
     26     // --------------------------------------------------------------------------------------
     27     // --------------------- FIR filter design (using window functions) ---------------------
     28     // --------------------------------------------------------------------------------------
     29 
     30     // Define options for plotting DSP data
     31     const std::string options = "phaseresp=True";
     32 
     33     // Define filter coefficient buffers of different sizes
     34     univector<fbase, 15> taps15;
     35     univector<fbase, 127> taps127;
     36     univector<fbase, 8191> taps8191;
     37 
     38     // Prepare window function expressions, data is not processed at this point
     39     // Hann window for 15-point filter
     40     expression_handle<fbase> hann = to_handle(window_hann(taps15.size()));
     41     // Kaiser window with alpha = 3.0 for 127-point filter
     42     expression_handle<fbase> kaiser = to_handle(window_kaiser(taps127.size(), 3.0));
     43     // Blackman-Harris window for 8191-point filter
     44     expression_handle<fbase> blackman_harris = to_handle(window_blackman_harris(taps8191.size()));
     45 
     46     // Design a 15-point lowpass FIR filter with a cutoff frequency of 0.15 using Hann window
     47     fir_lowpass(taps15, 0.15, hann, true);
     48 #if PYTHON_IS_INSTALLED
     49     // Plot the filter, frequency, and impulse response
     50     // internally plot_save calls Python (matplotlib and numpy packages used) to save SVG files
     51     plot_save("fir_lowpass_hann", taps15,
     52               options + ", phasearg='auto', title='15-point lowpass FIR, Hann window'");
     53 #endif
     54 
     55     // Design a 127-point lowpass FIR filter with a cutoff frequency of 0.2 using the Kaiser window
     56     // The resulting coefficients are in taps127
     57     fir_lowpass(taps127, 0.2, kaiser, true);
     58 #if PYTHON_IS_INSTALLED
     59     // Plot the filter, frequency, and impulse response
     60     plot_save("fir_lowpass_kaiser", taps127,
     61               options + ", phasearg='auto', title='127-point lowpass FIR, Kaiser window (\\alpha=3.0)'");
     62 #endif
     63 
     64     // Design a 127-point highpass FIR filter with a cutoff frequency of 0.2 using the Kaiser window
     65     // The resulting coefficients are in taps127
     66     fir_highpass(taps127, 0.2, kaiser, true);
     67 #if PYTHON_IS_INSTALLED
     68     // Plot the filter, frequency, and impulse response
     69     plot_save("fir_highpass_kaiser", taps127,
     70               options + ", phasearg='auto', title='127-point highpass FIR, Kaiser window (\\alpha=3.0)'");
     71 #endif
     72 
     73     // Design a 127-point bandpass FIR filter with cutoff frequencies of 0.2 and 0.4 using the Kaiser window
     74     fir_bandpass(taps127, 0.2, 0.4, kaiser, true);
     75 #if PYTHON_IS_INSTALLED
     76     // Plot the filter, frequency, and impulse response
     77     plot_save("fir_bandpass_kaiser", taps127,
     78               options + ", phasearg='auto', title='127-point bandpass FIR, Kaiser window (\\alpha=3.0)'");
     79 #endif
     80 
     81     // Design a 127-point bandstop FIR filter with cutoff frequencies of 0.2 and 0.4 using the Kaiser window
     82     fir_bandstop(taps127, 0.2, 0.4, kaiser, true);
     83 #if PYTHON_IS_INSTALLED
     84     // Plot the filter, frequency, and impulse response
     85     plot_save("fir_bandstop_kaiser", taps127,
     86               options + ", phasearg='auto', title='127-point bandstop FIR, Kaiser window (\\alpha=3.0)'");
     87 #endif
     88 
     89     // Design an 8191-point lowpass FIR filter with a cutoff frequency of 0.15 using the Blackman-Harris
     90     // window
     91     fir_lowpass(taps8191, 0.15, blackman_harris, true);
     92     // Initialize a buffer that is 8191+150 samples long and set it to zero
     93     univector<fbase, 8191 + 150> taps8191_150 = scalar(0);
     94 
     95     // Shift the filter coefficients by 150 samples
     96     taps8191_150.slice(150) = taps8191;
     97 
     98 #if PYTHON_IS_INSTALLED
     99     // Plot the filter, frequency, and impulse response
    100     // phasearg is used to get the correct phase shift (offset by 150 samples)
    101     plot_save(
    102         "fir_lowpass_blackman", taps8191_150,
    103         options +
    104             ", phasearg=4095+150, title='8191-point lowpass FIR, Blackman-Harris window', padwidth=16384");
    105 #endif
    106 
    107     // --------------------------------------------------------------------------------------
    108     // -------------------------- Using FIR filter as an expression -------------------------
    109     // --------------------------------------------------------------------------------------
    110 
    111     // Generate 10,000 samples of white noise
    112     univector<float> noise = truncate(gen_random_range(random_init(1, 2, 3, 4), -1.f, +1.f), 10000);
    113 
    114     // Apply the bandstop filter designed earlier to the noise
    115     univector<float> filtered_noise = fir(noise, fir_params{ taps127 });
    116 
    117 #if PYTHON_IS_INSTALLED
    118     // Plot the original noise and the filtered noise
    119     plot_save("noise", noise, "title='Original noise', div_by_N=True");
    120     plot_save("filtered_noise", filtered_noise, "title='Filtered noise', div_by_N=True");
    121 #endif
    122 
    123     // --------------------------------------------------------------------------------------
    124     // ------------------------------- FIR filter as a class --------------------------------
    125     // --------------------------------------------------------------------------------------
    126 
    127     // Redesign the 127-point bandpass FIR filter with cutoff frequencies of 0.2 and 0.4 using the Kaiser
    128     // window
    129     fir_bandpass(taps127, 0.2, 0.4, kaiser, true);
    130     // Initialize an FIR filter class with float input/output and fbase-typed taps
    131     filter_fir<fbase, float> fir_filter(taps127);
    132 
    133     // Apply the FIR filter to the noise data
    134     univector<float> filtered_noise2;
    135     fir_filter.apply(filtered_noise2, noise);
    136 
    137 #if PYTHON_IS_INSTALLED
    138     // Plot the results of the filtered noise
    139     plot_save("filtered_noise2", filtered_noise2, "title='Filtered noise 2', div_by_N=True");
    140 #endif
    141 
    142 #ifdef HAVE_DFT
    143     // --------------------------------------------------------------------------------------
    144     // ---------------------- Convolution filter (optimized using DFT) ----------------------
    145     // --------------------------------------------------------------------------------------
    146 
    147     // Initialize a convolution filter with float input/output and fbase-typed taps
    148     convolve_filter<fbase> conv_filter(taps127);
    149 
    150     // Apply the convolution filter to the noise data
    151     univector<fbase> filtered_noise3;
    152     conv_filter.apply(filtered_noise3, noise * fbase(1.0));
    153 
    154 #if PYTHON_IS_INSTALLED
    155     // Plot the results, same as filtered_noise2
    156     plot_save("filtered_noise3", filtered_noise3, "title='Filtered noise 3', div_by_N=True");
    157 #endif
    158 #endif
    159 
    160     println("SVG plots have been saved to svg directory");
    161 
    162     return 0;
    163 }