学期末,要求一篇关于滤波器的课程设计,不得不把已经忘记的差不多的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('时间');