Applying a Notch-filter to filter my MEP-data

5 次查看（过去 30 天）
Tom Siebring
Tom Siebring 2020-5-15
评论： Star Strider 2020-5-26
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:

请先登录，再回答此问题。

采纳的回答

Star Strider
Star Strider 2020-5-15
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.
  4 个评论
显示 2更早的评论
Tom Siebring
Tom Siebring 2020-5-26
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
Star Strider 2020-5-26
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!)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Digital Filter Design 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by