Remove 50 Hz noise

I am recording 6 channels from a force transducer, but there is a 50 Hz noise I want to remove. I was considering to use a simple solution as iirnotch, but I don't have this function on my Matlab R2015a. How can I create a filter (filtfilt and fir1?) to eliminate this 50 Hz without a signal processing toolbox?

Thank you

[Rik]

I added the 'no toolbox' tag, because this is almost trivial otherwise. My Google search (for 'bandstop filter') didn't turn up something either. There are several contributors more skilled in signal processing than me, so they might be able to help here.

[Paola]

Hi, Thanks for adding the tag. I hope some expert will help.

[Walter Roberson]

https://www.music.mcgill.ca/~gary/618/week1/signals.html

[Walter Roberson]

How wide can the stop band be? What is your sampling frequency?

[Paola]

20000 Hz. Precisely my fs is:

fs = Fdata.streams.Wav1.fs that is defined my acquisition system.

I have already tried usign a butterworth filter, but of course, it doesn't help so much.

 ForceData = double(ForceData);
       [B,A]=butter(2,50/fs/2); 
    for ch=1:6
        ForceData(ch,:) = filtfilt(B, A, ForceData(ch,:));
    end

[Walter Roberson]

That is not the correct way to create a bandstop butterworth filter.

https://www.mathworks.com/help/signal/ref/butter.html#bucse3u-ftype

[Paola]

What am I doing wrong?

According to the "cut off frequency" paragraph: "If Wn is scalar, then butter designs a lowpass or highpass filter with cutoff frequency Wn." So apparently for a bandstop I should use a two element vector. In this case what value should I put, to remove only the 50 Hz? And how do I choose the order?

ForceData = double(ForceData); 
[B,A]=butter(2,[50]/(fs/2),'stop'); 
for ch=1:6 
ForceData(ch,:) = filtfilt(B, A,ForceData(ch,:));  
end

Thank you

[Walter Roberson]

You cannot remove only 50 Hz exactly with a digital filter of that type. Instead you effectively construct a lowpass filter that starts to drop down near 50 Hz, cutting off high frequencies, "plus" (in some sense) a highpass filter that starts coming up near 50 Hz, cutting off low frequencies. The downwards dip heading into 50 Hz and the upwards dip heading after 50 Hz together diminish 50 Hz, along with diminishing frequencies approaching 50 Hz on other side.

The higher the order that you use, the sharper the cutoff can be, but also the longer the filter representation comes out when expressed as a convolution, so the processing rate is higher for lower orders.

[Paola]

Thank you for your clear explanation.

ex:
ForceData = double(ForceData); 
[B,A]=butter(2,[49 51]/(fs/2),'stop'); 
for ch=1:6 
ForceData(ch,:) = filtfilt(B, A,ForceData(ch,:));  
end

Is it correct in this way or should I express the cutoff frequencies in radians? I will try different orders to check what could be the best for me.

Thanks again

[Walter Roberson]

I think that looks okay, but I hardly do any filter work so I am not sure.

[Paola]

Hi, I tried this code but it doesn't change much. The ideal would be a notch but I don't have DSP toolbox and I am not able to calculate it by myself. Any suggestion to talk to a filter expert?

[Walter Roberson]

Tests I did with butter suggest that it is not very good at narrow band suppression :( With low orders it still left a fair bit of signal in the target range, and with orders from about 5 up, it created quite large transients.

[Paola]

Definitely! I am trying to write down a notch but I find it quite difficult since I am not an expert...

[Star Strider]

I thought you said that you do not have the Signal Processing Toolbox! I agree with Rik Wisselink that with it this is essentially trivial. However, you appear to be using its functions, so I will add this:

Fs = 20000;
Fn = Fs/2;
Wp = [47 53]/Fn;
Ws = [49 51]/Fn;
Rp =  1;
Rs = 50;
[n,Wn] = buttord(Wp,Ws,Rp,Rs);
[z,p,k] = butter(n,Wn,'stop');
[sos,g] = zp2sos(z,p,k);

figure
freqz(sos, 2^14, Fs)
set(subplot(2,1,1), 'XLim',[0 100])
set(subplot(2,1,2), 'XLim',[0 100])

The freqz call lets you see what the filter is doing in the frequency domain. There is also no need for the loop, because the filter functions operate column-wise. You simply need to transpose your matrix.

This should work:

ForceDataFilt = filtfilt(sos, g, ForceData');                       % Transpose & Filter
ForceDataFilt = ForceDataFilt';

I transposed ‘ForceDataFilt’ so it will match the rest of your code. Experiment with the filter passband and stopband frequencies to get the result you want. Remember that for a notch or bandstop filter, the stopband frequencies (here ‘Ws’) are within the passband limits, ‘Wp’. The passband and stopband attenuation values are ‘Rp’ and ‘Rs’ respectively.

[Walter Roberson]

Looks like you should be able to invoke designfilt in R2015a.

[Paola]

Actually, I don't have the DSP toolbox and I have an old version of the Signal Processing toolbox.

Answer 1

[Kouichi C. Nakamura]

This looks like an answer.

Remove the 60 Hz Hum from a Signal

https://uk.mathworks.com/help/signal/ug/remove-the-60-hz-hum-from-a-signal.html