用matlab做巴特沃斯低通滤波器
趁着暑假,做一个心电图的matlab实验,遇到了滤波器问题,网上代码比较杂乱,做了一个汇总整理。
主要做了一个简单的低通滤波器并以三角函数为例子进行低通滤波。
基本数据
fk=100 %采样频率N=1024 %采样个数n=-N/2:N/2-1f=n*fk/Nt=n/fky=sin(2*pi*10*t)
当然 f也可以表示成 f=linspace(-fk/2,fk/2,N);
第一步:做出原信号的频谱函数
注意:用fft函数作频谱分析,得到的是0~fk内的频谱 而用fftshift函数得到-fk/2~fk/2内的频谱
fft_y=fft(y);fftshift_y=fftshift(fft_y);f=linspace(-50,50,1024);plot(f,abs(fftshift_y));
第二步:做出巴特沃斯低通滤波器
用buttord函数求出阶数以及Wn
butter函数是求Butterworth数字滤波器的系数,在求出系数后对信号进行滤波时用filter函数
[g,Wn]=buttord(Wp,ws,Rp,Rs);
Wp=fp/(0.5fk)=20/50=0.4(通带截至频率)
Ws=fs/(0.5*fk)=30/50=0.6(阻带截止频率)
Rp=1; 通带最大衰减
Rs=30; 阻带最小衰减
代码如下:
[g,Wn]=buttord(0.4,0.6,1,30);
用butter 求出差分方程的系数b a
不懂差分方程及filter用法的朋友可以点以下链接,参考这位大佬的理解
差分方程介绍
[b,a]=butter(g,Wn);
至此 准备工作都做完了 ,来看看这滤波器的相频特性啥样子
滤波器的相频特性
[q,w]=freqz(b,a,256);plot(w*fs/(2*pi),abs(q))
所以说 咱写出来的滤波器长这个样子
20Hz为通带截至频率
30Hz为阻带截至频率
第三步:用滤波器过滤信号 并得出频谱图
用filter过滤
用filter函数就能把原信号y过滤成k函数 ,需要注意的是,它们都是时域函数,需要进行fft变换才能看到频谱图
k=filter(b,a,y);fft_k=fft(k);fftshift_k=fftshift(fft_k);plot(f,abs(fftshift_k));
过滤后的结果如下:
没错,没变!我们低通滤波器允许20Hz以下的信号不变,而我们原信号的频率就是10Hz!
对高频的信号的低通滤波
这里就不赘述了,对y=sin(2pi40*t)的信号,依然用上面的低通过滤的结果如下图:
不要看它长的怪,注意一下幅值,约等于0,即把40Hz的信号过滤掉了!
行文至此,点个赞呗亲!