filtering problem, need help

2020 2 1
1 个回答
回答已采纳
更新时间：2020 2 8
6 次查看（30 天）

0 个投票

Dear
I have a EMG signal of 30000x4. with sampling frequency of 10KHz.
I filter noise by below code
NotchFilter = bandstop(Data01Raw,[49.9 50.1],Fs);
sEmgFilter = bandpass(NotchFilter(:,1:2),[20 500],Fs);
iEmgFilter = bandpass(NotchFilter(:,3:4),[600 2000],Fs);
Data02Filtered = [sEmgFilter iEmgFilter];
Questions
in above code bandpass or bandstop occur by which filer like is it butterwork or cheby or ellip?
I need to fiter the signal using butterworth 2nd order. How can i do this?
Thank you

2 个评论
显示 无 隐藏 无

Walter Roberson
Walter Roberson 2020-2-1
编辑：Walter Roberson 2020-2-1
https://www.mathworks.com/help/signal/ref/bandpass.html
bandpass uses a fir filter if it can achieve the desgign goals with one, and otherwise uses a iir filter.
For Butterworth filters see https://www.mathworks.com/help/signal/ref/butter.html

请先登录，再进行评论。

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

 采纳的回答

Star Strider
Star Strider 2020-2-1

1 个投票

Request two outputs from the bandpass and bandstop functions. The second output is the digital filter object. Displaying the digital filter object in your Command Window will tell you everything you need to know about it.
For example —
Fs = 10E+3;
[NotchFilter, df] = bandstop(Data01Raw,[49.9 50.1],Fs);
then:
df
displays:
df =
digitalFilter with properties:
Coefficients: [8×6 double]
Specifications:
FrequencyResponse: 'bandstop'
ImpulseResponse: 'iir'
SampleRate: 10.0000e+003
StopbandAttenuation: 60.0000e+000
PassbandRipple2: 100.0000e-003
StopbandFrequency2: 50.0843e+000
PassbandRipple1: 100.0000e-003
PassbandFrequency2: 50.1000e+000
PassbandFrequency1: 49.9000e+000
StopbandFrequency1: 49.9157e+000
DesignMethod: 'ellip'
Use fvtool to visualize filter
Use designfilt to edit filter
Use filter to filter data

12 个评论
显示 10更早的评论 隐藏 10更早的评论

Ali Asghar
Ali Asghar 2020-2-1
dear
wrote code for following work on data 30000x4
1- apply notch filter of 50Hz
2- apply band pass of 2nd order, zero-phase shift of 20-500hz for column 1-2.
3- apply band pass of 2nd order, 60-2000hz for column 3-4.
[b,a] = butter(2,[49.9 50.1]/(10000/2),'stop');
filter1 = filter(b,a,dath001);
[b,a] = butter(2,[20 500]/(10000/2),'bandpass');
filter2 = filter(b,a,filter1(:,1:2));
[b,a] = butter(2,[60 2000]/(10000/2),'bandpass');
filter3 = filter(b,a,filter1(:,3:4));
filter_signal = [filter2 filter3];
I got the values without error but am i did it correctly?
Thank you
Star Strider
Star Strider 2020-2-1
I do not understand the overlapping bandpass filters. (The notch filter for the 50 Hz mains interference is obviously required.)
I would instead use elliptic filters, because they are shorter and therfore more computationally efficient, and among other things, they will produce the result you want.
Example —
[n,Wp] = ellipord([20 500]/Fn, [5 515]/Fn, 1, 50);
[z,p,k] = ellip(n, 1, 50, Wp,'bandpass');
[sos2,g2] = zp2sos(z,p,k);
figure('Name','Filter #2')
freqz(sos2, 2^24, Fs)
set(subplot(2,1,1), 'XLim',[0 1000]) % Optional
set(subplot(2,1,2), 'XLim',[0 1000]) % Optional
Then use the filtfilt function to do the actual filtering:
filter1 = filtfilt(sos2,g2,filter1(:,1:2));
and so for the others. The filtfilt function produces a phase-neutral result. so the signal has no phase distortion.
My apologies for the delay. I am having problems with MATLAB being very slow this morning for some reason.
Ali Asghar
Ali Asghar 2020-2-1
Dear
Thank you for your advice....
Actually, in column 1-2 i have sEMG data and in column 3-4 i have iEMG data.
As per my literature review butterworth give good result.
So kindly check that the above coding is correct or not.....
Star Strider
Star Strider 2020-2-1
It is, however it could be improved. Specifically, the filter order needs to be increased, since order 2 filters will not do what you want them to do. Use the buttord function to determine the optimal filter order, and specify a stopband as well as a passband in the filter design, such as I did with the second filter here. The second-order-section representation will be much more stable than the transfer-function representation. Also, it is best to use the freqz function to be certain the filters are doing what you want them to do.
Consider this version:
Fs = 10000;
Fn = Fs/2;
[z,p,k] = butter(2,[49.9 50.1]/(10000/2),'stop');
[sos1,g1] = zp2sos(z,p,k);
filter1 = filtfilt(sos1,g1,dath001);
[n,Wn] = buttord([20 500]/Fn, [5 515]/Fn, 1, 50);
[z,p,k] = butter(n, Wn,'bandpass');
[sos2,g2] = zp2sos(z,p,k);
filter2 = filtfilt(sos2,g2,filter1(:,1:2));
[z,p,k] = butter(2,[60 2000]/(10000/2),'bandpass');
[sos3,g3] = zp2sos(z,p,k);
filter3 = filtfilt(sos3,g3,filter1(:,3:4));
filter_signal = [filter2 filter3];
figure('Name','Filter #1')
freqz(sos1, 2^24, Fs)
set(subplot(2,1,1), 'XLim',[0 100]) % Optional
set(subplot(2,1,2), 'XLim',[0 100]) % Optional
figure('Name','Filter #2')
freqz(sos2, 2^24, Fs)
set(subplot(2,1,1), 'XLim',[0 700]) % Optional
set(subplot(2,1,2), 'XLim',[0 700]) % Optional
figure('Name','Filter #3')
freqz(sos3, 2^24, Fs)
set(subplot(2,1,1), 'XLim',[0 2500]) % Optional
set(subplot(2,1,2), 'XLim',[0 2500]) % Optional
Star Strider
Star Strider 2020-2-6
As always, my pleasure!
Ali Asghar
Ali Asghar 2020-2-7
Hey star
signal is 38000x1
I need to find the rms values of signal range 1-2000.
then rms of 2001-4000
till 38000 values of signal
means i get 19 values of rms.
How can i do this?
Star Strider
Star Strider 2020-2-7
Try this:
sigseg = buffer(signal, 2000);
RMS_sigseg = sqrt(mean(sigseg.^2));
See the documentation for buffer for informaiton about it.
Ali Asghar
Ali Asghar 2020-2-8
编辑：Ali Asghar 2020-2-8
Dear
It works but when i use it in function then it give only one value
function MAV=jMAV(X)
y=buffer(X, 2000);
MAV=mean(abs(y));
end
MAV=jMAV(x);
give one value, instead of giving multiple values.
Star Strider
Star Strider 2020-2-8
It gives one value for each column of the matrix. That is what it is defined to do.
Ali Asghar
Ali Asghar 2020-2-8
Dear, Can I find mutiple values of mav with this function?
Star Strider
Star Strider 2020-2-8
I am not certain what you are asking. You can have buffer break the signal up into shorter segments if you want.
With respect to calculating the RMS value, if you have R2016a or later, you can use the movmean function instead of mean.
This computes the RMS value over a sliding window of 200 samples:
RMS_sigseg = sqrt(movmean(sigseg.^2, 200));
It results in a matrix the same size as the original matrix.
Note that the code you posted gives the mean of the absolute value. This is not the same as RMS value.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Digital Filtering 的更多信息

产品

标签

Community Treasure Hunt

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

Start Hunting!

Translated by