firfilt: adding option to generate notch filter

examples/firfilt_cccf_notch_example.c

@@ -0,0 +1,101 @@
+//
+// firfilt_cccf_notch_example.c
+//
+// This example demonstrates how to create a notching non-recursive
+// (finite impulse response) filter.
+//
+
+#include <stdio.h>
+#include <math.h>
+#include <complex.h>
+
+#include "liquid.h"
+
+#define OUTPUT_FILENAME "firfilt_cccf_notch_example.m"
+
+int main() {
+    // options
+    unsigned int num_samples = 600;     // number of samples
+    unsigned int m           = 25;      // prototype filter semi-length
+    float        As          = 30.0f;   // prototype filter stop-band suppression
+    float        f0          = 0.15f;   // notch frequency
+
+    // design filter from prototype
+    firfilt_cccf q = firfilt_cccf_create_notch(m,As,f0);
+    firfilt_cccf_print(q);
+
+    // allocate memory for data arrays
+    float complex x[num_samples];   // original input
+    float complex y[num_samples];   // filtered signal
+
+    // generate input signal
+    unsigned int i;
+    for (i=0; i<num_samples; i++) {
+        x[i] = cexpf(_Complex_I*2*M_PI*0.037*i) +
+               cexpf(_Complex_I*2*M_PI*f0*   i);
+    }
+
+    // run notch filter
+    firfilt_cccf_execute_block(q, x, num_samples, y);
+
+    // destroy filter object
+    firfilt_cccf_destroy(q);
+
+    // 
+    // plot results to output file
+    //
+    FILE * fid = fopen(OUTPUT_FILENAME,"w");
+    fprintf(fid,"%% %s : auto-generated file\n", OUTPUT_FILENAME);
+    fprintf(fid,"clear all;\n");
+    fprintf(fid,"close all;\n");
+    fprintf(fid,"\n");
+    fprintf(fid,"num_samples=%u;\n",num_samples);
+    fprintf(fid,"x=zeros(1,num_samples);\n");
+    fprintf(fid,"y=zeros(1,num_samples);\n");
+
+    // save input, output arrays
+    for (i=0; i<num_samples; i++) {
+        fprintf(fid,"x(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(x[i]), cimagf(x[i]));
+        fprintf(fid,"y(%4u) = %12.4e + j*%12.4e;\n", i+1, crealf(y[i]), cimagf(y[i]));
+    }
+
+    // plot output
+    fprintf(fid,"t=0:(num_samples-1);\n");
+    fprintf(fid,"figure;\n");
+    fprintf(fid,"subplot(3,1,1);\n");
+    fprintf(fid,"  plot(t,real(x),'-','Color',[1 1 1]*0.5,'LineWidth',1,...\n");
+    fprintf(fid,"       t,imag(x),'-','Color',[0 0.2 0.5],'LineWidth',2);\n");
+    fprintf(fid,"  xlabel('time');\n");
+    fprintf(fid,"  ylabel('original input');\n");
+    fprintf(fid,"  legend('real','imag','location','northeast');\n");
+    fprintf(fid,"  axis([0 num_samples -2.2 2.2]);\n");
+    fprintf(fid,"  grid on;\n");
+    fprintf(fid,"subplot(3,1,2);\n");
+    fprintf(fid,"  plot(t,real(y),'-','Color',[1 1 1]*0.5,'LineWidth',1,...\n");
+    fprintf(fid,"       t,imag(y),'-','Color',[0 0.5 0.2],'LineWidth',2);\n");
+    fprintf(fid,"  xlabel('time');\n");
+    fprintf(fid,"  ylabel('notch filter output');\n");
+    fprintf(fid,"  legend('real','imag','location','northeast');\n");
+    fprintf(fid,"  axis([0 num_samples -2.2 2.2]);\n");
+    fprintf(fid,"  grid on;\n");
+    fprintf(fid,"subplot(3,1,3);\n");
+    fprintf(fid,"  nfft = 2^nextpow2(num_samples);\n");
+    fprintf(fid,"  f    = [0:(nfft-1)]/nfft-0.5;\n");
+    fprintf(fid,"  X    = 20*log10(abs(fftshift(fft(x,nfft))));\n");
+    fprintf(fid,"  Y    = 20*log10(abs(fftshift(fft(y,nfft))));\n");
+    fprintf(fid,"  plot(f,X,'-','Color',[1 1 1]*0.5,'LineWidth',1,...\n");
+    fprintf(fid,"       f,Y,'-','Color',[0.5 0 0],  'LineWidth',2);\n");
+    fprintf(fid,"  xlabel('time');\n");
+    fprintf(fid,"  ylabel('PSD (dB)');\n");
+    fprintf(fid,"  legend('original','filtered','location','northeast');\n");
+    fprintf(fid,"  axis([-0.5 0.5 -20 60]);\n");
+    fprintf(fid,"  grid on;\n");
+
+    // close output file
+    fclose(fid);
+    printf("results written to '%s'\n", OUTPUT_FILENAME);
+
+    printf("done.\n");
+    return 0;
+}
+

include/liquid.h

@@ -2501,6 +2501,14 @@ FIRFILT() FIRFILT(_create_rect)(unsigned int _n);                           \
 FIRFILT() FIRFILT(_create_dc_blocker)(unsigned int _m,                      \
                                       float        _As);                    \
                                                                             \
+/* Create notch filter from prototype                                   */  \
+/*  _m  : prototype filter semi-length such that filter length is 2*m+1 */  \
+/*  _As : prototype filter stop-band attenuation [dB], _As > 0          */  \
+/*  _f0 : center frequency for notch, _fc in [-0.5, 0.5]                */  \
+FIRFILT() FIRFILT(_create_notch)(unsigned int _m,                           \
+                                 float        _As,                          \
+                                 float        _f0);                         \
+                                                                            \
 /* Re-create filter object of potentially a different length with       */  \
 /* different coefficients. If the length of the filter does not change, */  \
 /* not memory reallocation is invoked.                                  */  \

makefile.in

@@ -1420,6 +1420,7 @@ example_programs :=						\
 	examples/firdecim_crcf_example				\
 	examples/firfarrow_rrrf_example				\
 	examples/firfilt_cccf_example				\
+	examples/firfilt_cccf_notch_example			\
 	examples/firfilt_crcf_example				\
 	examples/firfilt_crcf_dcblocker_example			\
 	examples/firfilt_rrrf_example				\

src/filter/src/firfilt.c

@@ -224,6 +224,47 @@ FIRFILT() FIRFILT(_create_dc_blocker)(unsigned int _m,
     return FIRFILT(_create)(h, h_len);
 }
 
+// create notch filter
+FIRFILT() FIRFILT(_create_notch)(unsigned int _m,
+                                 float        _As,
+                                 float        _f0)
+{
+    // validate input
+    if (_m < 1 || _m > 1000) {
+        fprintf(stderr,"error: %s:%u, firfilt_%s_create_notch(), filter semi-length (%u) must be in [1,1000]\n",
+                __FILE__, __LINE__, EXTENSION_FULL, _m);
+        exit(1);
+    } else if (_As <= 0.0f) {
+        fprintf(stderr,"error: %s:%u, firfilt_%s_create_notch(), prototype stop-band suppression (%12.4e) must be greater than zero\n",
+                __FILE__, __LINE__, EXTENSION_FULL, _As);
+        exit(1);
+    } else if (_f0 < -0.5f || _f0 > 0.5f) {
+        fprintf(stderr,"error: %s:%u, firfilt_%s_create_notch(), notch frequency must be in [-0.5,0.5]\n",
+                __FILE__, __LINE__, EXTENSION_FULL, _As);
+        exit(1);
+    }
+
+    // create float array coefficients and design filter
+    unsigned int h_len = 2*_m+1;
+    float hf[h_len];
+    liquid_firdes_dcblocker(_m, _As, hf);
+
+    // copy coefficients to type-specific array and mix to appropriate frequency
+    TC h[h_len];
+    unsigned int i;
+    for (i=0; i<h_len; i++) {
+        float phi = 2.0f * M_PI * _f0 * ((float)i - (float)_m);
+#if TC_COMPLEX
+        h[i] = cexpf(_Complex_I*phi) * (TC) hf[i];
+#else
+        h[i] = cosf(phi) * (TC) hf[i];
+#endif
+    }
+
+    // return filter object and return
+    return FIRFILT(_create)(h, h_len);
+}
+
 // re-create firfilt object
 //  _q      :   original firfilt object
 //  _h      :   new coefficients [size: _n x 1]