【肌电信号】基于带通滤波肌电信号处理含Matlab源码

1 简介

在光学信息处理中，滤波是一种我们实验中经常用到的方式。通常情况下， 可以利用空间滤波器来完成对特定波长的筛选工作。

2 部分代码

clear all;
file_name = 'Data1.txt';    % 文件名
fid =fopen(file_name,'r');  %打开文件
h=fgets(fid);               %去掉第一行；
C=textscan(fid,'%*f%f%f%f%f');%去掉第一列；并将其他行列写入C
fclose(fid); %关闭文件
Data=[C{1}';C{2}';C{3}';C{4}'];%将C中的数据取出放入数组Data1
fid =fopen(['new_' file_name],'w'); %打开新文件夹
fprintf(fid,'%-.10e    %-.10e   %-.10e   %-.10e\r\n  ',Data); %将取出的数据放入新文件夹，纯数据
fclose(fid);%关闭新文件夹
%load 'new_Data1.txt';
t=['new_' file_name];
[newData]=textread(t,'','delimiter',' ');
x1=newData(:,1);
x2=newData(:,2);
x3=newData(:,3);
x4=newData(:,4);
figure(1);
subplot(4,1,1)
plot(x1);
axis([0 512 0.5 1.7]);grid on;
% figure;
subplot(4,1,2)
plot(x2);
axis([0 512 0.5 1.7]);grid on;
% figure;
subplot(4,1,3)
plot(x3);
axis([0 512 0.5 1.7]);grid on;
% figure;
subplot(4,1,4)
plot(x4);
axis([0 512 0.5 1.7]);grid on;
N=512;
M=1; %调节因子
fs=1000;
%数据点数
%时间序列
n=0:N-1;
f=n*fs/N;    %频率序列
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%信号
figure(2);
y1=fft(x1,N);    %对信号进行快速Fourier变换
mag1=abs(y1);     %求得Fourier变换后的振幅
subplot(4,1,1) ,plot(f(1:N/2),mag1(1:N/2));
axis([0 500 0 20]);grid on;
y2=fft(x2,N);    %对信号进行快速Fourier变换
mag2=abs(y2);    %求得Fourier变换后的振幅
subplot(4,1,2), plot(f(1:N/2),mag2(1:N/2));
axis([0 500 0 20]);grid on;
y3=fft(x3,N);    %对信号进行快速Fourier变换
mag3=abs(y3);     %求得Fourier变换后的振幅
subplot(4,1,3) ,plot(f(1:N/2),mag3(1:N/2));
axis([0 500 0 20]);grid on;
y4=fft(x4,N);    %对信号进行快速Fourier变换
mag4=abs(y4);     %求得Fourier变换后的振幅
subplot(4,1,4), plot(f(1:N/2),mag4(1:N/2));
axis([0 500 0 20]);grid on;
%=====%part3%===========%带通滤波器%========%
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(3)
freqz(b1,a1)
title('n=4 Butterworth Bandpass Filter');
%=====%part3%===========%带通滤波器%========%
%=====%part4%===========%50hz陷波器%=======%
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(4)
freqz(b2,a2)
title('n=2 Butterworth Bandstop Filter');
%=====%part4%===========%50hz陷波器%=======%
%=====%part5%=============%数字滤波%========%
x1=filter(b1,a1,x1);
x1=filter(b2,a2,x1);
x2=filter(b1,a1,x2);
x2=filter(b2,a2,x2);
x3=filter(b1,a1,x3);
x3=filter(b2,a2,x3);
x4=filter(b1,a1,x4);
x4=filter(b2,a2,x4);
figure(5);
subplot(4,1,1)
plot(x1);grid on;
axis([0 512 -0.5 0.5])
subplot(4,1,2)
plot(x2);grid on;
axis([0 512 -0.5 0.5])
subplot(4,1,3)
plot(x3);grid on;
axis([0 512 -0.5 0.5])
subplot(4,1,4)
plot(x4);grid on;
axis([0 512 -0.5 0.5])
%=====%part5%=============%数字滤波%========%
figure(6);
y1=fft(x1,N);    %对信号进行快速Fourier变换
mag1=abs(y1);     %求得Fourier变换后的振幅
subplot(4,1,1) ,plot(f(1:N/2),mag1(1:N/2));
axis([0 500 0 11]);
grid on;
y2=fft(x2,N);    %对信号进行快速Fourier变换
mag2=abs(y2);    %求得Fourier变换后的振幅
subplot(4,1,2), plot(f(1:N/2),mag2(1:N/2));
axis([0 500 0 10]);
grid on;
y3=fft(x3,N);    %对信号进行快速Fourier变换
mag3=abs(y3);     %求得Fourier变换后的振幅
subplot(4,1,3) ,plot(f(1:N/2),mag3(1:N/2));
axis([0 500 0 10]);
grid on;
y4=fft(x4,N);    %对信号进行快速Fourier变换
mag4=abs(y4);     %求得Fourier变换后的振幅
subplot(4,1,4), plot(f(1:N/2),mag4(1:N/2));
axis([0 500 0 10]);
grid on;
%=====%part6%=============%验证效果%========%
iemg1=sum(abs(x1))/length(x1);
iemg2=sum(abs(x2))/length(x2);
iemg3=sum(abs(x3))/length(x3);
iemg4=sum(abs(x4))/length(x4);
   %求积分肌电值
rms1=sqrt(sum(x1.^2)/length(x1));
rms2=sqrt(sum(x2.^2)/length(x2));
rms3=sqrt(sum(x3.^2)/length(x3));
rms4=sqrt(sum(x4.^2)/length(x4));
  %求均方根值
L1=length(x1);
cx1=xcorr(x1,'unbiased');
cxk1=fft(cx1,L1);
px1=abs(cxk1);%求功率谱密度
pxx1=10*log10(px1);
f1=(0:L1-1)*fs/L1;
df1=fs/L1;
p1=(sum(px1(1:L1/2-1))+sum(px1(1:L1/2)))/2.*df1;
pf1=(sum(px1(1:L1/2-1).*[1:L1/2-1]'.*df1)+sum(px1(1:L1/2).*[1:L1/2]'.*df1))/2*df1;
MPF1=pf1/p1 ;%求平均功率频率
N1=1;pp1=0;
while abs(pp1-p1/2)>(px1(N1)+px1(N1+1))/2*df1
    pp1=pp1+(px1(N1)+px1(N1+1))/2*df1;
    N1=N1+1;
end
n_1=(N1+N1+1)/2;
MF1=df1*n_1; %求中值频率
L2=length(x2);
cx2=xcorr(x2,'unbiased');
cxk2=fft(cx2,L2);
px2=abs(cxk2);%求功率谱密度
pxx2=10*log10(px2);
f2=(0:L2-1)*fs/L2;
df2=fs/L2;
p2=(sum(px2(1:L2/2-1))+sum(px2(1:L2/2)))/2.*df2;
pf2=(sum(px2(1:L2/2-1).*[1:L2/2-1]'.*df2)+sum(px2(1:L2/2).*[1:L2/2]'.*df2))/2*df2;
MPF2=pf2/p2 ;%求平均功率频率
N2=1;pp2=0;

3 仿真结果

4 参考文献

[1]李向明. 基于Matlab的表面肌电信号处理软件设计与开发[D]. 上海体育学院.

博主简介：擅长智能优化算法、神经网络预测、信号处理、元胞自动机、图像处理、路径规划、无人机等多种领域的Matlab仿真，相关matlab代码问题可私信交流。

部分理论引用网络文献，若有侵权联系博主删除。