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

时间：2021-06-20 17:59:09 阅读：0 评论：0 收藏：0 [点我收藏+]

一、简介

1 滤波器
滤波器按照频带划分可以分为：低通滤波器(LPF)、高通滤波器(HPF)、带通滤波器(BPF)、带阻滤波器(BEF)。其中射频天线领域主要采用带通滤波器(BPF)。

2 带通滤波器(BPF)
技术图片
①：RF BPF：从天线中取出期望频带的有效信号，滤除不需要频带的电磁波信号和噪声。
②：RF BPF：用下级混频器(降频混频器)进行频率转换时，防止在IF频带中引入干扰信号，用BPF将不期望频率成分滤除。
③：IF BPF：从混频器进行频率转换所产生的许多频率成分中，用BPF只取出期望的IF信号，也具有限制接收信号频带的作用。
④：IF BPF：用混频器(升频转换器)进行频率转换时，用BPF防止在RF频带中引入干扰信号，也具有限制发射信号频带的作用。
⑤：RF BPF：从混频器进行频率转换所产生的许多频率成分中，用BPF只取出期望的RF信号。
⑥：RF BPF：用BPF除去功率放大器放大时失真所产生的寄生部分，以及大功率信号的输入在天线开关时所产生的寄生部分，以免寄生成分通过天线发射出去。

3 带通滤波器(BPF)选择：
(1)对于低频电路，选择LC滤波器。LC滤波器通常用于基频领域。
(2)对于100MHz~1000MH频率范围，选择SAW滤波器。SAW滤波器通常用于IF BPF领域。
(3)对于2GHz频率电路，选择介质滤波器。介质滤波器通常用于RF BPF领域。

4 带通滤波器(BPF)详解：
核心思想：BPF基本上是谐振电路。谐振电路分为串联谐振电路和并联谐振电路，对于串联谐振电路，若L和C无损耗，则在谐振频率时阻抗变为零；对于并联谐振电路，谐振频率时阻抗变为无穷大。对于串联谐振电路可知，L越大，衰减特性越陡；对于并联谐振电路可知，L越小C越大，衰减特性越陡。
谐振电路的谐振频率：
技术图片

串联谐振与并联谐振的组合电路：

5 介质滤波器

六微带线制作滤波器

二、源代码

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);