filter.m

针对EMG（肌电信号）的频率范围

PI=3.14159;
%=====%part1%=========%原始信号的产生%======%
%y=sin（wx），其中：w=2∏/T=2∏f
%对于模拟频率，需要对其采样，根据采样点数N与采样频率fs决定，y（k）=sin（wx*k/fs）,等于是将原来的信号1s采样2000个，而采样频率为
%4000hz，等于说每2hz采一个信号
%============copyright by Then=============%
y=2*sin(2*PI*50)+sin(2*PI*40)+sin(2*PI*49)+sin(2*PI*60)+sin(2*PI*80)+sin(2*PI*100)+sin(2*PI*120)+sin(2*PI*150)+sin(2*PI*160)+sin(2*PI*300)+sin(2*PI*400)+sin(2*PI*500)+sin(2*PI*600);
%=====%part1%=========%原始信号的产生%======%
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%=====%part2%============%信号采样%=========%
%============copyright by Then=============%
N=1000;
fs=2000;
dt=1/fs;
for k=1:N
    y(k)=2*sin(2*PI*50*k*dt)+sin(2*PI*40*k*dt)+sin(2*PI*49*k*dt)+sin(2*PI*60*k*dt)+sin(2*PI*80*k*dt)+sin(2*PI*100*k*dt)+sin(2*PI*120*k*dt)+sin(2*PI*150*k*dt)+sin(2*PI*160*k*dt)+sin(2*PI*300*k*dt)+sin(2*PI*400*k*dt)+sin(2*PI*500*k*dt)+sin(2*PI*600*k*dt);
end
%=====%part2%============%信号采样%=========%
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%=====%part3%===========%带通滤波器%========%
%硬件：截止频率为1000hz，根据来奎斯特定理，fs>=2fc(工程上一般取4fc)
%Wp=[Wp1 Wp2]/来奎斯特频率 （这里fn=fs/2=4fc/2=2fc）
%通带衰减小于1db，阻带衰减大于50db
%============copyright by Then=============%
Wp = [50 150]/1000; Ws = [10 500]/1000;
Rp = 1; Rs = 50;
[n,Wn] = buttord(Wp,Ws,Rp,Rs)
[b1,a1] = butter(n,Wn);
figure(1)
freqz(b1,a1,1000,2000)
title('n=4 Butterworth Bandpass Filter')
%=====%part3%===========%带通滤波器%========%
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%=====%part4%===========%50hz陷波器%=======%
%============copyright by Then=============%
Wp = [10 500]/1000; Ws = [49 51]/1000;
Rp = 1; Rs = 50;
[n,Wn] = buttord(Wp,Ws,Rp,Rs)
[b2,a2] = butter(n,Wn,'stop');
figure(2)
freqz(b2,a2,1000,2000)
title('n=2 Butterworth Bandstop Filter')
%=====%part4%===========%50hz陷波器%=======%
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%=====%part5%=============%数字滤波%========%
%============copyright by Then=============%
y1=filter(b1,a1,y);
y2=filter(b2,a2,y1);
%=====%part5%=============%数字滤波%========%
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%==========================================&
%=====%part6%=============%验证效果%========%
%============copyright by Then=============%
y=fft(y,N);                                %原始信号周期谱
pyy=y.*conj(y);
figure(3)
plot(1:2:2*(N/2),pyy(1:N/2));
xlabel('Hz')
y1=fft(y1,N);                              %滤波后的周期谱
pyy1=y1.*conj(y1);
figure(4)
plot(1:2:2*(N/2),pyy1(1:N/2));
xlabel('Hz')
y2=fft(y2,N);                              %滤波后的周期谱
pyy2=y2.*conj(y2);
figure(5)
plot(1:2:2*(N/2),pyy2(1:N/2));
xlabel('Hz')
%=====%part6%=============%验证效果%========%