Applying a Notch-filter to filter my MEP-data

Hi everyone,

I've been trying to apply a 50 Hz Notchfilter to my data (Motor evoked potentials) for some time but I can't get it to work.

One example of what I've tried is this:

Fs = 5000 %Sampling frequency is 5000 Hz
t = (0:length(trialtime)-1)/Fs; % My trace is 400 ms long
plot(t,trialdata(:,1))
d = designfilt('bandstopiir', 'FilterOrder',2,'HalfPowerFrequency1',49,'HalfPowerFrequency2',51,'DesignMethod','butter','SampleRate',Fs);
fvtool(d,'Fs',Fs)
buttLoop = filtfilt(d,trialdata(:,1));
plot(t,trialdata(:,1),t,buttLoop)

However, doing this I get the following error:

Warning: SOS structures not supported for FieldTrip drop-in filtfilt

Usage: y=filtfilt(b,a,x)

Whatever other way of filtering I try, it always comes back to using the filtfilt function and I get the error above.

It must have something to do with the size of b but I don't understand how to correct this or where it goes wrong.

function y = filtfilt(b, a, x)

if (nargin ~= 3)
    error('Usage: y=filtfilt(b,a,x)');
end

if (min(size(b)) > 1)
    warning('SOS structures not supported for FieldTrip drop-in filtfilt');
    error('Usage: y=filtfilt(b,a,x)');
end

My data looks as following:

trialdata = array of my datapoints (voltage change)

trialtime = time array over which the voltage change is measured

plot(trialtime,trialdata) looks like this:

Answer 1

[Star Strider]

Using second-order-section representation is certainly the best option. However if FieldTrip (that I have no experience with) forces you to use transfer-function representation, I doubt a second-order Butterworth design is going to be adequate.

Here is an elliptic design that provides steep rolloffs to the rejection region and appears to be stable as a transfer function. It designs a third-order filter:

Fs = 5000;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [46  54]/Fn;                                           % Normalised Passband (Passband = 46 Hz To 54 Hz)
Ws = [49  51]/Fn;                                           % Normalised Stopband (Passband = 49 Hz To 51 Hz)
Rp =  1;                                                    % Passband Ripple/Attenuation
Rs = 50;                                                    % Stopband Ripple/Attenuation
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                          % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'stop');                      % Elliptic Filter
[b,a] = zp2tf(z,p,k);                                       % Convert To Transfer Function

figure
freqz(b,a,2^16,Fs)
set(subplot(2,1,1), 'XLim',[40 60])                         % ‘Zoom’ Frequency Axis
set(subplot(2,1,2), 'XLim',[40 60])                         % ‘Zoom’ Frequency Axis

This is as steep as I could get the filter rolloffs using transferfunction representation, and produce a stable filter.

[Tom Siebring]

Thank you very much, this way it works on my data.

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Tom Siebring]

I have an additional question.

Besides a notch-filter, I also need to apply a bandpass filter [20 2000] on my data.

I've tried the function y = bandpass(x,fpass,fs) but it gave me the same error as what I got with the notch-filter. So I figured (and since I saw your reply to someone elses question about a bandpassfilter https://nl.mathworks.com/matlabcentral/answers/361348-how-i-applly-a-bandpass-filter-in-a-signal) I can apply a filter in the similar way as what you told me for the notch filter.

so I did the following:

Fs = 5000;                                                  % Sampling Frequency (Hz)
Fn = Fs/2;                                                  % Nyquist Frequency (Hz)
Wp = [20  2000]/Fn;                                           % Normalised Passband (Passband = 46 Hz To 54 Hz)
Ws = [19  2001]/Fn;                                           % Normalised Stopband (Passband = 49 Hz To 51 Hz)
Rp =  1;                                                    % Passband Ripple/Attenuation
Rs = 50;                                                    % Stopband Ripple/Attenuation
[n,Wp] = ellipord(Wp, Ws, Rp, Rs);                          % Calculate Elliptic Filter Optimum Order
[z,p,k] = ellip(n, Rp, Rs, Wp,'bandpass');                      % Elliptic Filter
[b,a] = zp2tf(z,p,k);                                       % Convert To Transfer Function
figure
freqz(b,a,2^16,Fs)

filtered_trials = filtfilt(b,a,trial_notch);

However, this did not work on my data. I'm not sure if all the values I filled in are correct.

Plotting the filter looked like this:

applying this filter to my data made the traces look way worse than before:

Do you know how to get this right also?

Thank you!

[Star Strider]

As always, my pleasure!

I’ll do my best!

When I tried it, I was able to get a stable transfer function filter with passbands of [20 2000] and stopbands of [10 2050]. No other changes were necessary. (Good choice in changing my previous code to use an elliptic filter instead!)