How do I eliminate multiple 60Hz harmonics within 1 set of data in an ECG

I have an ECG signal and I'm trying to filter out the 60 Hz harmonics at 60, 120, 180 etc. I am currently using a digital biquad filter. I am currently running into the problem that they are not within one data set so when I plot them with the ECG singal it's not plotting one graph with both 60 AND 120 minimized, it's plotting 1 graph with 60 Hz minimized and 1 graph with 120 Hz minimized on top of each other. I've tried a for loop and just rewriting the code using hold on to continue to plot. How do I get them so that they are all in within the same set of data and not different plots on top of each other for each? Here is some of my current code that is showing JUST the 60 Hz harmonic minimized:

    dz= .0005;   % d = 1/Q
    dp= 0.3;
    
    w= 2*pi*(60);  % notch frequency in radians/sec
    B =[1 dz*w w^2];
    A =[1 dp*w w^2];
    h = tf(B,A);
    hd =c2d(h,1/Fs,'tustin'); %Fs= 2000 Hz
    filter_x = filter(hd.Numerator{1},hd.Denominator{1},x);
    plot(filter_x)

Answer 1

[Star Strider]

The Signal Processing Toolbox is best for this. See: Remove the 60 Hz Hum from a Signal (link). You’ll need to design a separate filter for each harmonic.

If you have R2018a or later, the bandstop (link) function makes such filter design and implementation easy.

[Kathleen Nicole Rice]

How do i implement each seperate filter for the harmonics all onto 1 graph/1 set of data though? I've read a lot of articles regarding it but many only include removing 1 harmonic at a time. Would I have to remove one, save the data, and then remove another? Is there no easier way to do this?

[Star Strider]

The easiest way would be to filter the 60 Hz fundamental with one filter, then use that output to filter the 120 Hz harmonic. Do the same for higher harmonics.

You could design a FIR filter to filter the fundamental and the first two harmonics at the same time:

Fs = 1000;                                                  % Sampling Frequency (Hz)
fcuts = [56 58 62 64 116 118 122 124 176 178 182 184];      % Frequencies
mags = [1 0 1 0 1 0 1];                                     % Passbands & Stopbands
devs = [0.05 0.01 0.05 0.01 0.05 0.01 0.05];                % Tolerances
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs,Fs);          % Kaiser Window FIR Specification
n = n + rem(n,2);
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');           % Filter Realisation

figure
freqz(hh,1,2^14,Fs)
set(subplot(2,1,1), 'XLim',[0 200]);                        % Zoom Frequency Axis
set(subplot(2,1,2), 'XLim',[0 200]);                        % Zoom Frequency Axis

Experiment with the ‘fcuts’ frequency vector to get the respoinse you want. If you add frequencies to it, you will have to change the ‘mags’ and ‘devs’ vectors as well.

This should be robust to different sampling frequencies, providing that the sampling frequency of your signal is at least 400 Hz.

I always use filtfilt to do the filtering. (With some FIR filters, filter works. This is likely one such filter, since it appears to have a linear phase response.)

My apologies for the delay. I’ve not designed FIR filters in a while, and had to refresh my memory.

[Kathleen Nicole Rice]

No problem- thank you so much for your help!! :)

[Star Strider]

As always, my pleasure!