# 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

- Comments
- Write a Comment Select to add a comment

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-]

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

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

title('Magnitude of Biquads')

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

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];

.

.

.

I believe it should be declared as: double acoef[8, 3], bcoef[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

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

http://www.advsolned.com/downloads/ASN15-DOC002.pdf

Please let us know what you think!!!!

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.

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.

Registering will allow you to participate to the forums on ALL the related sites and give you access to all pdf downloads.