# 60 Hz EMF filter implemented on Arduino

11 Dec

Now the digital 60 Hz filter is up and running on the Arduino. The goal was to take in all frequencies of EMF, reject anything that was not 60 Hz, and show the amplitude of the 60 Hz component. Here’s a graph of what the Arduino is doing–producing the green filtered signal from the blue raw signal. It also does a moving average over three cycles so there is a bit of lag, about 0.05 seconds. Thanks to Chris Isert for getting the Arduino interrupts set up so we could control the sampling rate (500 Hz).

Raw and filtered 60 Hz signal sampled at 500 Hz and sent over serial link

It’s very sensitive to whether the antenna is being touched–higher values are when it is touched, then it drops off when not touched as shown about halfway through the graph. But if the user is also touching ground while touching the antenna, the values drop to nearly zero.
There should really be a lowpass RC filter on the EMF detector so that the filter does not get fooled by high frequency harmonics above half the sampling frequency. This may be why the green signal is a little bit higher than the tips of the blue peaks.

//Digital 60Hz filter by CKH 12/11–report amplitude of 60Hz electromagnetic field over serial interface
int inPin = 5; //analog 5
volatile int isNewVal=0; //flag for new read, variables changed in an interrupt handler should be declared volatile
volatile int val; //where to store info from analog 5
//FilterCoefficients for 32 Tap 60HzFilter, I multiplied all of these by 10000 to have ints, this was necessary for calc to be fast enough for 500 Hz
int b[]={ -61, 105, -152, -270, -237, -7, 354, 637, 602, 171, -471,-944,-914,-330,505,1104,1104,505,-330,-914,-944,-471,171,602,637,354,-7,-237,-270,-152,105,-61};
long latestval[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0}; //Accumulate latest readings here, need longint or it wraps around
long myval;//this will hold the latest value of the filtered signal
float mrec[]={0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};//holds latest 3 cycles of filtered 60Hz signal at 500Hz (25 samples)
float myamp=0; //latest amplitude of the 60Hz filtered, rectified and moving-averaged signal
int myintamp; //this will be an integer version of myamp for saving.
unsigned int mytime=0;//used for checking how long it takes to do things in the loop

void setup()//set up interrupts for 500Hz done by Chris Isert, it looks like OCR1A is important in setting the sampling freq

//Note that if you change the sampling frequency from 500 Hz, you also need to change the filter coefficients in b.

{

Serial.begin(115200);

// initialize timer1

noInterrupts(); // disable all interrupts

TCCR1A = 0;

TCCR1B = 0;

TCNT1 = 0;

OCR1A = 32000; // compare match register 16MHz/1/500Hz

TCCR1B |= (1 << WGM12); // CTC mode

TCCR1B |= (1 << CS10); // divide-by-1 prescaler

TIMSK1 |= (1 << OCIE1A); // enable timer compare interrupt

interrupts(); // enable all interrupts

}

ISR(TIMER1_COMPA_vect) // timer compare interrupt service routine

{

isNewVal=1;//CKH: this is a flag telling the loop that there’s a new value

}

void loop()

{
if (isNewVal){
for (int i=0;i<31;i++){
latestval[i]=latestval[i+1]; //slide the latest 32 raw values back by one in the array
}
//now do the convolution with the filter coefficients in b, to get the latest filtered value
myval=0;
for (int i=0; i<32; i++){
myval=myval+latestval[i]*b[i];
}
for(int i=0;i<24;i++){
mrec[i]=mrec[i+1];//shift the latest 24 filtered and rectified values back by one in preparation for adding the newest one on the end
}
//then calculate the most recent filtered and rectified value
mrec[24]=abs(0.126*myval)/10000;//rectified (absolute value) signal, times pi/25, it gives the amplitude of the positive-only ac signal scaled for moving avg of latest 25 pts
//divided by 10000 to get back the original filter coefficients which were not integers
myamp=myamp-mrec[0]+mrec[24];//latest moving avg over 3 cyc of 60 hz at 500 hz sampling rate
myintamp=(int)myamp;//this conversion to int makes it much faster to send the value over the serial connection than a float, and code size is smaller
Serial.print(val);
Serial.print(“,”);
Serial.println(myintamp);
isNewVal=0;//clear new value flag until next interrupt
//mytime=micros(); //check whether we are staying at 500Hz, cant spend too much time number crunching between readings
//Serial.println(mytime,DEC);

}
}

———————————

What the 32 tap filter is doing, is throwing out other components besides 60Hz. It would be ideal to test this with a real electrical signal, but meanwhile MATLAB can offer some information. In the following plot I created a 5-volt, 200 Hz sinewave and added it with a 2 volt, 60Hz in the top graph, then ran the filter in MATLAB to produce the bottom graph. The 200 Hz is gone, but the Arduino still has to figure out that the amplitude of the remaining 60Hz is 2V and send that information over the serial connection.

Top: original signal, bottom: signal filtered with digital 60 Hz bandpass filter

Below is some MATLAB code for setting up this filter (some setup is from a previous post) and creating test data to generate plots like the one above:

———————————

Fs=500;%Sampling Freq in Hz

Nyq=Fs/2;%Maximum freq i can detect

CenterFreqNorm=60/Nyq;%Normalized 60 Hz

f=[ 0 CenterFreqNorm-0.1 CenterFreqNorm-0.02 CenterFreqNorm+0.02 CenterFreqNorm+0.1 1];%set up frequency regions

a=[0 0 1 1 0 0];%set up desired amplitudes

B=firpm(31,f,a); %calc the filter coefficients for a 32 tap fir filter

[h,w]=freqz(B,1,250); %then get the freq response to see what it does

%then plot it

%plot(w/pi*250,20*log10(abs(h)))

%Now lets make up some sensor data with 60 Hz, 200 Hz, and 10 Hz and random

%noise, see what happens when the B filter is convolved with it

MyTime=0:(2^14-1); %make it a length a power of 2 in case that is useful for anything

MyElapsedTime=(2^14-1)/Fs; %about 32 seconds at 500 Hz

My60Hz=2*sin(60*MyTime*2*pi/Fs);% a 2Vpp 60 Hz signal

My200Hz=0.5*sin(200*MyTime*2*pi/Fs); %and some other signals to filter out

My10Hz=2*sin(10*MyTime*2*pi/Fs);

noise=0.5*randn(1,length(MyTime)); %bwahahaha mess it up

MySignal=My60Hz+My200Hz+My10Hz+noise;

%MySignal=My60Hz+10*My200Hz+My10Hz;

%MySignal=My60Hz;

MyC=conv(MySignal,B);%do the filter in the time domain

%check it out

subplot(2,1,1)

plot(MySignal(1:500))

subplot(2,1,2)

plot(MyC(1:500))

%You still have to rectify and find the amplitude

%What if the orig signal were >0 like we saw with emf detector

%MyPosSignal=MySignal;

%MyPosSignal(MyPosSignal<0)=0;

%figure ()

%plot(MyPosSignal(1:500))

%MyPosC=conv(MyPosSignal,B);%filter the lopped off signal

%subplot(2,1,1)

%plot(MyPosSignal(1:500))

%subplot(2,1,2)

%plot(MyPosC(1:500))

%Result, filtering the lopped signal gives a 60Hz sine with half the p-p

%amplitude. The filtered sine is centered at 0.

%what is a quick way to calc the amplitude of a sine wave–rectify it

%multiply neg values by -1, keep a running average

MyRecC=MyC;%create a copy to rectify

MyRecC(MyRecC<0)=-1*MyRecC(MyRecC<0);%flip it positive if it’s neg

MyAvg=filter(ones(1,100),100,MyRecC);%this will avg over about 12 cycles of the

%60 Hz sine wave with an evenly weighted 100-point moving average

%the avg value of a half sine is 2*amplitude/pi

%so to get the original 60 Hz amplitude back, you mult avg by pi/2 or 1.57, maybe not super important

MyCalc60HzAmp=MyAvg*1.57;

hold on

plot (MyCalc60HzAmp(1:500),’r-‘)%works, goes right atop the extracted 60Hz signal peaks