肌电信号处理程序+原始数据

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设定初始参数%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% fs=2048; %采样频率 N=512; %采样点数 n=0:N-1; t=n/fs; %时间 f1=100; f2=300; %低频 f3=800; %高频 s=cos(2*pi*f1*t)+cos(2*pi*f2*t)+cos(2*pi*f3*t);%假定信号 subplot(121);plot(t,s); title('输入信号');xlabel('t/s');ylabel('幅度');%未滤波时 时域波形 sfft=fft(s); %傅里叶变换 subplot(122); plot((1:length(sfft)/2)*fs/length(sfft),2*abs(sfft(1:length(sfft)/2))/length(sfft)); title('信号频谱');xlabel('频率/Hz');ylabel('幅度'); %未滤波时 频域波形 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%设计巴特沃斯低通滤波器%%%%%%%%%%%%%%%%%%%%%%%% Wp=900/fs;Ws=1000/fs; [n,Wn]=buttord(Wp,Ws,1,50); %阻带衰减大于50db,通带纹波小于1db %估算得到Butterworth低通滤波器的最小阶数N和3dB截止频率Wn [a,b]=butter(n,Wn); [h,f]=freqz(a,b,'whole',fs); %求数字低通滤波器的频率响应 f=(0:length(f)-1)*fs/length(f); %进行对应的频率转换 figure(2); plot(f(1:length(f)/2),abs(h(1:length(f)/2))); %绘制巴特沃斯幅频响应图 title('巴特沃斯低通滤波器');xlabel('频率/Hz');ylabel('幅度'); grid; sF=filter(a,b,s); %叠加函数s经过低通滤波器以后的新函数 figure(3); subplot(121); plot(t,sF); %绘制叠加函数s经过低通后时域图形 title('输出信号');xlabel('t/s');ylabel('幅度'); SF=fft(sF); subplot(122); plot((1:length(SF)/2)*fs/length(SF),2*abs(SF(1:length(SF)/2))/length(SF)); title('低通滤波后频谱');xlabel('频率/Hz');ylabel('幅度'); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%时域特征值计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% iemg1=sum(abs(sF))/length(sF); rms1=sqrt(sum(sF.^2)/length(sF)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%频域特征值计算%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% L1=length(sF); cx1=xcorr(sF,'unbiased'); cxk1=fft(cx1,L1); px1=abs(cxk1); %求功率谱密度 pxx1=10*log10(px1); figure(4) f1=(0:L1-1)*fs/L1; plot(f1(1:L1/2),pxx1(1:L1/2)) figure(4) xlabel('频率/Hz');ylabel('功率谱/dB'); title('平均功率谱图'); grid on %做功率谱图 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; %求中值频率 MF=MF1 MPF=MPF1 iemg=iemg1 rms=rms1