# How to design an interpolated notch filter in Python to filter powerline harmonics

I have a a time domain data set that records the magnetic field vs time, which must be processed to reveal an embedded signal. This data also contains power line harmonics (i.e. multiples of 60 Hz) that must be removed from the signal. I am processing the data with Python and am using the numpy, scipy.signal and scipy.fftpack modules to filter out information.

Instead of using a traditional notch/comb filter I am writing a function that will transform the data to the frequency domain with an array of values corresponding to amplitude and another array corresponding to frequency. The arrays are also reformatted to exist properly in the frequency spectrum. The code then determines where the frequency peaks are and the array indice that corresponds to each peak. Finally, the code analyzes each peak to see if it exists as a multiple of 60 Hz (i.e. a power line harmonic). If any of the data points corresponds to the first 50 multiples of 60 Hz, the code finds the amplitude and frequency array values for indice-1 and indice+1 and then does a linear interpolation between them to find a new amplitude. The new value of amplitude replaces the old value, thus eradicating the frequency peaks corresponding to each power line harmonic, without the depressions caused by normal notch filters. The Python function is shown below.

When I plot the frequency response I get a very nice looking plot that has clearly erased the harmonic peaks. However, when I used an ifft function to transform the data back to the time domain, I see that the power line harmonics all still exist, but the amplitude has been reduced. I suspect that python is only working on the even (real) series and is neglecting the odd (imaginary) data, but I am not sure what is causing this. Again, when I plot the frequency domain info, there are no harmonic peaks (60, 180, 300, 440, 600 Hz, etc...), but when I plot the time domain it is clearly there. Any help would be appreciated. Also, I am not including the peakdet and the linear_interpolate functions because I have tested those and am confident the problem does not lie with them.

def Interpolate_US_Power_Harmonics(XData,YData):
from scipy.fftpack import fft
import numpy as np

# Determines time between sampling assuming uniform intervals
Sample_Freq = XData[1] - XData[0]

# Resize x and y axis for frequency domain
Y_Data = fft(YData,n=YData.size)
array_length = len(Y_Data)
Y1_Data = 20*np.log10(Y_Data/1) # transform to decibels
Y1_Data = 2.0/array_length*np.abs(Y1_Data[0:array_length/2])
XTicks = np.linspace(0.0,1.0/(2*Sample_Freq),array_length/2)

# - Determine the frequency peaks (maxima) and troughs (minima) and their array
#   positions.  Odd array positions correspond to peaks and even points
#   correspond to minimas
maxima, minima, position = peakdet(Y1_Data, 200*Sample_Freq, XTicks)

# Determine array position of frequency peaks
new_position = []
for i in range(0,len(position)):
if i% 2 != 0:
new_position.append(position[i])

# - Determine if a peak is a powerline harmonic (multiple of 60 Hz) and if so
#   replace magnitude array point with one interpolated from the previous
#   and next array points.
for x in range(0,len(maxima)):
frequency = 60.0
for y in range(100):
if maxima[x][0] > frequency - 1 and maxima[x][0] < frequency + 1:
X1 = XTicks[new_position[x-1]-1]; X2 = XTicks[new_position[x-1]+1]
Y1 = Y_Data[new_position[x-1]-1]; Y2 = Y_Data[new_position[x-1]+1]
X3 = maxima[x][0]
Y_Data[new_position[x-1]] = Linear_Interpolate(X1,X2,Y1,Y2,X3)
print frequency, XTicks[new_position[x-1]]
frequency = frequency + 60

# Return frequency domain data to main program
return Y_Data

# Main program
import scipy.fftpack.ifft
# meta code to show time and amplitude array as well as plotting
Time      = insert data here
Amplitude = insert data here

# used pylab to plot original time and frequency domain signals

# Produce new frequency series filtered
new_data = Interpolate_US_Power_Harmonics(Time,Amplitude)

# - used pylab to plot filtered data.  The frequency domain showed no harmonics.
#   but the time series plots still contained 60, 180, etc.. harmonics with
#   a reduced amplitude

• Any reason why you aren't using the typical comb filter that you mentioned instead of this ad hoc approach? – Jason R Mar 4 '16 at 1:13
• Two reasons. One: The lesser of the two reasons is in the general frequency response. A comb filter will to some extent (in many cases a large extent) cause a depression in the frequency spectrum around the notched frequency, which can in some cases have an adverse affect on the signal. Two: the most important reason is I have been directed by my supervisor to use this filter. – Jon Mar 4 '16 at 6:14
• That's an important reason. I'm not sure it's a good idea, though. If you want to truly suppress the harmonics, you'll need to place zeros on the harmonic frequencies. What kind of output do you expect from the algorithm you described? It sounds about right to me; this approach will only attenuate the harmonics somewhat. Spectral modification like this also can often wreak havoc with the time-domain structure of the signal. – Jason R Mar 4 '16 at 12:57
• In the freq. domain I expect the harmonic peaks to be erased leading to a smooth transition from Y[n-1] to Y[n+1], with a value of Y[n] that is an interpolate of the two values with there corresponding X[n] values. In the time domain I expect to have a flat line with a spike due to the signal we are trying to detect and of course with a noise floor. Im convinced the theory is sound, it is just that the algorithm is not correctly handling the odd and even series correctly and I dont know why. – Jon Mar 4 '16 at 17:51
• The theory is described in Mewett, D. T., Nazeran, H., and Reynolds, K. J., "Removing Power Line Noise from Recorded EMG," Proceedings of the 23rd Annual IEEE/EMBS Conference, Oct. 25-28, which can be found on line. – Jon Mar 4 '16 at 17:55