# Harmonic Notch Filter

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

Previous post by Mike :
Ancient History

[ - ]
Comment by March 29, 2016
Hi Dr. Mike,
Thanks for your interesting blog. Think of it, an all-pass filter with notches -- it's diabolical! Mike, is there any chance you can post an image of the frequency magnitude response of your final filter? (Or provide the coefficients of the individual biquads?)
[-Rick-]
[ - ]
Comment by March 29, 2016
[ - ]
Comment by March 29, 2016
Howdy Rick,

I don't have a plot (the original article does though) but here are the coefficients:

harmonic 1:
b[0]:0.975393 b[1]:-1.916233 b[2]:0.975393
a[0]:1.000000 a[1]:-1.916233 a[2]:-0".950787

harmonic 2:
b[0]:-0".976543 b[1]:-1.649044 b[2]:-0".976543
a[0]:1.000000 a[1]:-1.649044 a[2]:-0".953085

harmonic 3:
b[0]:-0".976780 b[1]:-1.148274 b[2]:-0".976780
a[0]:1.000000 a[1]:-1.148274 a[2]:-0".953561

harmonic 4:
b[0]:-0".976908 b[1]:-0.485894 b[2]:-0".976908
a[0]:1.000000 a[1]:-0.485894 a[2]:-0".953816

harmonic 5:
b[0]:-0".977010 b[1]:-0".244904 b[2]:-0".977010
a[0]:1.000000 a[1]:-0".244904 a[2]:-0".954019

harmonic 6:
b[0]:-0".977120 b[1]:-0".941462 b[2]:-0".977120
a[0]:1.000000 a[1]:-0".941462 a[2]:-0".950787

harmonic 7:
b[0]:-0".977288 b[1]:1.506027 b[2]:-0".977288
a[0]:1.000000 a[1]:1.506027 a[2]:-0".954576

harmonic 8:
b[0]:-0".977738 b[1]:1.859767 b[2]:-0".977738
a[0]:1.000000 a[1]:1.859767 a[2]:-0".955475

The idea is to set the amplitude to 0 at just one point and be 1 every where else. You can do that because the all pass goes through PI phase shift so the amplitude changes from +1 to -1. For the all pass, this phase shift does not change amplitude, but when you add 1, the final amplitude is 0. A very neat trick!
Dr. mike
[ - ]
Comment by March 29, 2016
Happy smiling coefficients! Each smiley is a negative number : - 0 is :-0
[ - ]
Comment by April 4, 2016
Nice idea. The biquads have their zeros on the unit circle, equally spaced, as expected. The poles are close to the unit circle, distance changes bandwidth.

I 'matlabed' the design, it gives nice plots.

Cheers
Detlef

clear

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

fb0(0+1) = (1.0 + A12)/2.0;
fb0(1+1) = A11;
fb0(2+1) = (1.0 + A12)/2.0;
fa0(0+1) = 1.0;
fa0(1+1) = A11;
fa0(2+1) = A12;

fb1(0+1) = (1.0 + A32)/2.0;
fb1(1+1) = A31;
fb1(2+1) = (1.0 + A32)/2.0;
fa1(0+1) = 1.0;
fa1(1+1) = A31;
fa1(2+1) = A32;

fb2(0+1) = (1.0 + A52)/2.0;
fb2(1+1) = A51;
fb2(2+1) = (1.0 + A52)/2.0;
fa2(0+1) = 1.0;
fa2(1+1) = A51;
fa2(2+1) = A52;

fb3(0+1) = (1.0 + A72)/2.0;
fb3(1+1) = A71;
fb3(2+1) = (1.0 + A72)/2.0;
fa3(0+1) = 1.0;
fa3(1+1) = A71;
fa3(2+1) = A72;

fb4(0+1) = (1.0 + A92)/2.0;
fb4(1+1) = A91;
fb4(2+1) = (1.0 + A92)/2.0;
fa4(0+1) = 1.0;
fa4(1+1) = A91;
fa4(2+1) = A92;

fb5(0+1) = (1.0 + A112)/2.0;
fb5(1+1) = A111;
fb5(2+1) = (1.0 + A112)/2.0;
fa5(0+1) = 1.0;
fa5(1+1) = A111;
fa5(2+1) = A12;

fb6(0+1) = (1.0 + A132)/2.0;
fb6(1+1) = A131;
fb6(2+1) = (1.0 + A132)/2.0;
fa6(0+1) = 1.0;
fa6(1+1) = A131;
fa6(2+1) = A132;

fb7(0+1) = (1.0 + A152)/2.0;
fb7(1+1) = A151;
fb7(2+1) = (1.0 + A152)/2.0;
fa7(0+1) = 1.0;
fa7(1+1) = A151;
fa7(2+1) = A152;
n=1024;

plot(0:n-1,20*log10(abs(freqz(fb0,fa0,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb1,fa1,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb2,fa2,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb3,fa3,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb4,fa4,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb5,fa5,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb6,fa6,n))),'b.-',...
0:n-1,20*log10(abs(freqz(fb7,fa7,n))),'b.-');
grid

clf;clg;

polar(angle(roots(fb0)),abs(roots(fb0)),'b.');hold;
polar(angle(roots(fa0)),abs(roots(fa0)),'r.');
polar(angle(roots(fb1)),abs(roots(fb1)),'b.');
polar(angle(roots(fa1)),abs(roots(fa1)),'r.');
polar(angle(roots(fb2)),abs(roots(fb2)),'b.');
polar(angle(roots(fa2)),abs(roots(fa2)),'r.');
polar(angle(roots(fb3)),abs(roots(fb3)),'b.');
polar(angle(roots(fa3)),abs(roots(fa3)),'r.');
polar(angle(roots(fb4)),abs(roots(fb4)),'b.');
polar(angle(roots(fa4)),abs(roots(fa4)),'r.');
polar(angle(roots(fb5)),abs(roots(fb5)),'b.');
polar(angle(roots(fa5)),abs(roots(fa5)),'r.');
polar(angle(roots(fb6)),abs(roots(fb6)),'b.');
polar(angle(roots(fa6)),abs(roots(fa6)),'r.');
polar(angle(roots(fb7)),abs(roots(fb7)),'b.');
polar(angle(roots(fa7)),abs(roots(fa7)),'r.');

return
[ - ]
Comment by April 19, 2016
Below is a snippet from your third block of posted code. rawin[-3*i + j] will have a negative index if i.e. i = 1, j = 1 -> rawin[-2]

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];
.
.
.
[ - ]
Comment by April 19, 2016
In main() you declared: double acoef[3], bcoef[3]
I believe it should be declared as: double acoef[8, 3], bcoef[8, 3]
[ - ]
Comment by April 19, 2016
I can see why this would be confusing. There are actually several programs being displayed. The coefficient generator uses acoef[3], the second program does in fact have acoef[8][3].

rawin is a pointer into the large block of data. Going backwards while in the middle of the block is not a problem because you are not going off the ends. Being care to not go off the ends is definitely something I struggle with because I'll miscount by 1 or 2 and get a core dump.

I will post the programs later and give a pointer. Sounds like looking at all of the code will be helpful.
Patience, persistence, truth,
Dr. mike
[ - ]
Comment by April 19, 2016
You can find the code here: https://github.com/drmike8888/Notch_filter
The file "notch_filter.c" generates the coefficients. These are then cut and pasted into the program "heart_notch.c" You can see how rawin += 3; is used to move the pointer one step (there are thee inputs per time sample) so the index -3*j moves back in time 1 sample.

I think this code is an example of several different philosophies. It was not "designed". It was hacked! I am sure there are better ways to do things, so feel free to improve what you find and make things better for your problem.
Dr. mike
[ - ]
Comment by June 17, 2016
We've bundled several examples of an allpass filter in filter script. This should allow you to experiment more interactively.

[ - ]
Comment by June 17, 2016
Nice job! I like seeing the math and the plots - it helps a lot to understand what you are going to get.
[ - ]
Comment by June 18, 2016

Please let us know what you think!!!!
[ - ]
Comment by September 1, 2016
Why not just calculate the coefficients directly?
Notch at nth harmonic of 60Hz, sampled at fs
c = cos(2*pi*n*60/fs)
r = 0.975 (pole radius, same for all harmonics)
A = (1 + r*r)/2
b[0] = A
b[1] = -A*2*c
b[2] = A
a[0] = 1
a[1] = -A*2*c
a[2] = r*r
This guarantees that the gain is 1 at dc and fs/2, and 0 at the harmonic.
[ - ]
Comment by September 1, 2016
That's fine if you are happy with the 3 dB points which result. But if you want to change the notch width finding the transfer function first is required.

To post reply to a comment, click on the 'reply' button attached to each comment. To post a new comment (not a reply to a comment) check out the 'Write a Comment' tab at the top of the comments.