Harmonic Notch Filter

• Tips and Tricks

My basement is covered with power lines and florescent lights which makes collecting ECG and EEG data  rather difficult due to the 60 cycle hum.  I found the following notch filter to work very well at eliminating the background signal without effecting the highly amplified signals I was looking for. 

The notch filter is based on the a transfer function with the form $$H(z)=\frac{1}{2}(1+A(z))$$ where A(z) is an all pass filter. The original paper [1] describes a method to combine all the notch locations to get one transfer function.  For use with a DSP,  I prefer to have biquad coefficients and this is easy to do if each notch is taken as a second order IIR.

A second order all pass filter has the form$$A(z) = \frac{a_2 + a_1 z^{-1} + z^{-2}}{1 + a_1 z^{-1} + a_2 z^{-2}}$$ Note the symmetry of the coefficients.  The idea is to put the notches where the all pass filter phase goes through $\theta(\omega_n)=-(2 n - 1)\pi$.  The width of the notch is dealt with by setting the all pass filter to $\theta(\omega_n-\frac{BW_n}{2})=-(2 n - 1)\pi + \frac{\pi}{2}$.  For the simple second order notch we can define the following variables:$$\omega_n = 2\pi j\frac{f_n}{ f_s}$$ $$\omega_{BW}=2\pi\frac{f_{BW}}{f_s}$$ $$\omega_1=\omega_n - \frac{\omega_{BW}}{2}$$ $$\beta_1 = \omega_1 - \frac{\pi}{4}$$ $$\beta_2 = \omega_n - \frac{\pi}{2}$$ $$p_0 = \tan(\beta_1)$$ $$p_1 = \tan(\beta_2)$$ where $f_s$ is the sample rate, $j f_n$ is the $j$th harmonic of the notch frequency $f_n$ and $f_{BW}$ is the desired spread between 3 dB down amplitudes at the notch.  

A matrix of the form $$q = \left[\matrix{\sin(\omega_1)-p_0 \cos(\omega_1) & \sin(2 \omega_1) - p_0\cos(2 \omega_1) \\ \sin(\omega_n) - p_1 \cos(\omega_n) & \sin(2 \omega_n) - p_1\cos(2 \omega_n)}\right]$$ is then created and the inverse is computed so that the coefficients are found from $$c = q^{-1} p$$ or explicitly $$c_0=q^{-1}_{0 0} p_0 + q^{-1}_{0 1} p_1$$ $$c_1 = q^{-1}_{1 0} p_0 + q^{-1}_{1 1} p_1$$ 

The final filter is then given with $$H(z) = \frac{1 + c_1/2 + c_0 z^{-1} + (1 + c_1/2)z^{-2}}{1 + c_0 z^{-1} + c_1 z^{-2}}$$

Here is an example which uses this form to generate the first 8 odd harmonic notches of 60 Hz.   Obviously there are better ways to output the coefficients, but for my purposes cutting and pasting was simple enough. Also note that creating a 50 Hz harmonic notch filter only requires change a 6 to a 5 on one line of code. In this example, the sample rate was 2kHz.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

main()
{
  int i, j;
  double omega, acoef[3], bcoef[3], hnr, hni, hdr, hdi, mag;
  double c1, c2, f, s1, s2;
  double omegan1, omegabw1;
  double omega1, omega2, theta1, theta2, beta1, beta2;
  double p[2], q[2][2], det, qinv[2][2];
  double am[2], tmp;

  for(j=1; j<16; j+=2)
  {
    omegan1 = 2.0*M_PI*60.0*j/2000.0;
    omegabw1 = 2.0*M_PI*15.0/2000.0;
    omega1 = omegan1 - omegabw1/2.0;
    omega2 = omegan1;
    theta1 = -M_PI/2.0;
    theta2 = -M_PI;
    beta1 = theta1/2.0 + omega1;
    beta2 = theta2/2.0 + omega2;
    p[0] = tan(beta1);
    p[1] = tan(beta2);
    q[0][0] = sin(omega1) - p[0]*cos(omega1);
    q[0][1] = sin(2.0*omega1) - p[0]*cos(2.0*omega1);
    q[1][0] = sin(omega2) - p[1]*cos(omega2);
    q[1][1] = sin(2.0*omega2) - p[1]*cos(2.0*omega2);
    det = q[0][0]*q[1][1] - q[0][1]*q[1][0];
    qinv[0][0] = q[1][1]/det;
    qinv[0][1] = -q[0][1]/det;
    qinv[1][0] = -q[1][0]/det;
    qinv[1][1] = q[0][0]/det;
    am[0] = qinv[0][0]*p[0] + qinv[0][1]*p[1];
    am[1] = qinv[1][0]*p[0] + qinv[1][1]*p[1];

    printf("%d coefficients are %13.10lf %13.10lf\n", j, am[0], am[1]);
  }
}

I then used the output of this program with the following to create the biquad filter coefficients:

#define A11 (-1.9162329361)
#define A12 (0.9507867324)
#define A31 (-1.6490444725)
#define A32 (0.9530853152)
#define A51 (-1.1482741905)
#define A52 (0.9535607368)
#define A71 (-0.4858941845)
#define A72 (0.9538156135)
#define A91 (0.2449035617)
#define A92 (0.9540193348)
#define A111 (0.9414624597)
#define A112 (0.9542403314)
#define A131 (1.5060265873)
#define A132 (0.9545758641)
#define A151 (1.8597673176)
#define A152 (0.9554750804)
:
:
// notch filter

  bcoef[0][0] = (1.0 + A12)/2.0;
  bcoef[0][1] = A11;
  bcoef[0][2] = (1.0 + A12)/2.0;
  acoef[0][0] = 1.0;
  acoef[0][1] = A11;
  acoef[0][2] = A12;

  bcoef[1][0] = (1.0 + A32)/2.0;
  bcoef[1][1] = A31;
  bcoef[1][2] = (1.0 + A32)/2.0;
  acoef[1][0] = 1.0;
  acoef[1][1] = A31;
  acoef[1][2] = A32;

  bcoef[2][0] = (1.0 + A52)/2.0;
  bcoef[2][1] = A51;
  bcoef[2][2] = (1.0 + A52)/2.0;
  acoef[2][0] = 1.0;
  acoef[2][1] = A51;
  acoef[2][2] = A52;

  bcoef[3][0] = (1.0 + A72)/2.0;
  bcoef[3][1] = A71;
  bcoef[3][2] = (1.0 + A72)/2.0;
  acoef[3][0] = 1.0;
  acoef[3][1] = A71;
  acoef[3][2] = A72;

  bcoef[4][0] = (1.0 + A92)/2.0;
  bcoef[4][1] = A91;
  bcoef[4][2] = (1.0 + A92)/2.0;
  acoef[4][0] = 1.0;
  acoef[4][1] = A91;
  acoef[4][2] = A92;

  bcoef[5][0] = (1.0 + A112)/2.0;
  bcoef[5][1] = A111;
  bcoef[5][2] = (1.0 + A112)/2.0;
  acoef[5][0] = 1.0;
  acoef[5][1] = A111;
  acoef[5][2] = A12;

  bcoef[6][0] = (1.0 + A132)/2.0;
  bcoef[6][1] = A131;
  bcoef[6][2] = (1.0 + A132)/2.0;
  acoef[6][0] = 1.0;
  acoef[6][1] = A131;
  acoef[6][2] = A132;

  bcoef[7][0] = (1.0 + A152)/2.0;
  bcoef[7][1] = A151;
  bcoef[7][2] = (1.0 + A152)/2.0;
  acoef[7][0] = 1.0;
  acoef[7][1] = A151;
  acoef[7][2] = A152;

This "brute force" method of programming is not efficient, but for what I was doing at the time being quick and dirty got the job done with the least thinking.  I then used to coefficients on 3 sets of data collected simultaneously and filtered all the data at once:

  for(k=0; k<3; k++)
    for(i=0; i<8; i++)
      for(j=0; j<3; j++)
        stage[k][i][j] = 0.0;
  for(j=0; j<3; j++)  // loop over each curve
  {

      stage[j][0][2] = 0.0;
      for(i=0; i<3; i++)
          stage[j][0][2] += rawin[-3*i + j]*bcoef[0][i];
      stage[j][0][2] -= acoef[0][1]*stage[j][0][1] + acoef[0][2]*stage[j][0][0];
      for(k=1; k<7; k++) // loop over each harmonic
      {
         stage[j][k][2] = 0.0;
         for(i=0; i<3; i++)
             stage[j][k][2] += stage[j][k-1][2-i]*bcoef[k][i];
          stage[j][k][2] -= acoef[k][1]*stage[j][k][1] + acoef[k][2]*stage[j][k][0];
      }
      passptr[j] = 0.0;
      for(i=0; i<3; i++)
          passptr[j] += stage[j][6][2-i]*bcoef[7][i];
      passptr[j] -= acoef[7][1]*passptr[j-3] + acoef[7][2]*passptr[j-6];
      for(k=0; k<7; k++)
      {
           stage[j][k][0] = stage[j][k][1];
           stage[j][k][1] = stage[j][k][2];
      }
  }

The "stage" variable contains the state of each biquad block as the signal passes though.  These are zeroed out before any processing happens.  "rawin" is the source data and this enters the first biquad.  All the subsequent biquads get data from the previous biquad so that the final signal has been filtered with multiple notches.  The output is saved at "passptr", which actually gets incremented by 3 in this case because there are 3 samples for every time stamp.  

At the very end, each biquad is advanced in time by shifting the internal storage one step.  The biquads are then ready for the next input sample.

This is only one way to approach a set of notch filters.  Using the method of the original paper, a complete $q$ matrix for all the notches can be created at once, and the transfer function can be computed.  One can then use standard methods to break the resulting transfer function into biquads.  In any event, this method of computing a notch filter is pretty slick.

[1] Soo-Chang Pei and Chien-Cheng Tseng, "IIR Multiple Notch Filter Design Based on Allpass Filter", 1996 IEEE TENCON. Digital Signal Processing Applications , pg 267-272

Sign in