双线性变换法的心跳信号滤波器（带通滤波器）

学期末，要求一篇关于滤波器的课程设计，不得不把已经忘记的差不多的MATLAB拿过来 “预习预习”。😂以此篇来记录一下完成这篇小论文的过程。（等以后再回头看看以前的自己🤞）。

滤波器

滤波器就是过滤掉我们不想要的，留下我们想要的信号。滤波器可以分为时域滤波和频域滤波。此设计中我们采用频域滤波。

设计详解

傅里叶变换

通过傅里叶变换，从时域转到频域。
MATLAB内置函数 fft() 快速傅里叶变换。傅里叶变换后的最高频率为采样频率的一半。（可参考抽样定理）

带通滤波器

允许特定频段的波通过同时屏蔽其他频段的。
通带为可以通过的频域范围；阻带是屏蔽的信号的频率范围。该设计中设置通带的频率为[1825 1950]；阻带的频率为[1675 2100]；这个频率并不是真正的截止频率，而是一个范围值。需要通过 buttord() 函数求得滤波器阶数n 和截止频率wn；

双线性变换

从s域转换到z域。
bilinear() 函数
[bz az]=bilinear(b a Fs)
b,a分别是s域的系统函数的分子，分母的系数序列；
bz,az分别是z域的系统函数的分子，分母的系数序列；
Fs是采样频率

模拟滤波器

现将设置的通带的截止频率[1825 1950]，阻带的截止频率[1675 2100]转换成角频率，用2pi截止频率/Fs 进行转换；再将角频率归一化处理，在转换成模拟滤波器的指标，即2/T*tan(截止频率/2)，得到的结果用于下面的函数。

[n wn]=buttord(wp,ws,rp,rs,‘s’);
wp是通带的截止频率；ws是阻带的截止频率；rp是通带的最大衰减；rs是阻带的最小衰减。
[b a]=butter(n,wn,‘s’)
这样就求得了模拟滤波器的系统函数。

数字滤波器

通过bilinear()函数将模拟滤波器转换成数字滤波器

滤波

filtfilt() 函数
outdata=filtfilt(bz az Fs)
得到的结果outdata就是滤波后的结果。

运行结果

源代码

这部分代码直接运行是不可以的，因为还缺少一个数据包heartsignal_part.mat，换成自己的数据就可以运行了。

clc;
clear;
%[data,fs]=audioread('13262.wav');   %读取音频信号
load heartsignal_part.mat;    %导入数据
data=data_pata;
len=length(data);
FS=8000; %设置抽样频率
 
%对信号进行傅里叶变换
fftdata=fft(data);
fftmax=abs(fftdata(1:len/2));%计算幅值
w=(1:len/2)*FS/len;


%绘制原始信号的信息图像
figure(1);
subplot(2,1,1);
plot(data);
title('原始信号');
ylabel('幅值');
xlabel('时间');
subplot(2,1,2);
plot(w,fftmax);
title('抽样信号的频率特性');
ylabel('幅值');
xlabel('频率');

%设置带通滤波器
rp=3;    %通带最大衰减
rs=40;   %阻带最小衰减

%频率转化为角频率
tonglow1=2*pi*1825/FS;   %通带最小     
tonghigh1=2*pi*1950/FS;    %通带最大
zulow1=2*pi*1675/FS;     %阻带最小
zuhigh1=2*pi*2100/FS;      %阻带最大

%归一化处理
%预畸变处理，将角频率参数转为模拟滤波器参数
tonglow=2*FS*tan(tonglow1/2);
tonghigh=2*FS*tan(tonghigh1/2);
zulow=2*FS*tan(zulow1/2);
zuhigh=2*FS*tan(zuhigh1/2);

tong=[tonglow tonghigh];
zu=[zulow zuhigh];

[n,Wn]=buttord(tong,zu,rp,rs,'s');
[b,a]=butter(n,Wn,'s');
[bz,az]=bilinear(b,a,FS);  %双线性变换


%滤波器的幅频，相频特性
[H1,h1]=freqz(b,a);
[H,h]=freqz(bz,az);
figure(2);
subplot(2,1,1);
plot(h,abs(H));
title('数字滤波器幅频特性曲线');
ylabel('幅值');
xlabel('频率');
subplot(2,1,2);
plot(h1,abs(H1));
title('模拟滤波器幅频特性曲线');
ylabel('幅值');
xlabel('频率');

figure(9);
plot(h,angle(H),'color','r');
hold on;
plot(h1,angle(H1),'color','g');
title('相频特性曲线');
ylabel('相位角');
xlabel('频率');
legend('数字滤波器','模拟滤波器');

%滤波
outdata=filtfilt(bz,az,data);

%对滤波后信号进行傅里叶变换
fftoutdata=fft(outdata);
fftoutmax=abs(fftoutdata(1:len/2));%计算幅值


figure(3);
subplot(2,1,1);
plot(w,fftmax);
title('滤波前信号频谱');
ylabel('幅值');
xlabel('频率');
subplot(2,1,2);
plot(w,fftoutmax);
title('滤波后信号频谱');
ylabel('幅值');
xlabel('频率');



figure(4);
subplot(2,1,1);
plot(data);
title('滤波前信号');
ylabel('幅值');
xlabel('时间');
subplot(2,1,2);
plot(outdata);
title('滤波后信号');
ylabel('幅值');
xlabel('时间');