MATLAB: Applying a Notch-filter to filter the MEP-data

50hzbandstopfilternotch

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:

Best Answer

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.
Related Question