用matlab做巴特沃斯低通滤波器

用matlab做巴特沃斯低通滤波器

趁着暑假,做一个心电图的matlab实验,遇到了滤波器问题,网上代码比较杂乱,做了一个汇总整理。
主要做了一个简单的低通滤波器并以三角函数为例子进行低通滤波。

基本数据

fk=100  %采样频率
N=1024  %采样个数
n=-N/2:N/2-1
f=n*fk/N
t=n/fk
y=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的信号过滤掉了!
行文至此,点个赞呗亲!

  • 73
    点赞
  • 273
    收藏
    觉得还不错? 一键收藏
  • 7
    评论
巴特沃斯低通滤波器是一种常用的数字滤波器,可用于信号处理中的低频信号滤波。MATLAB提供了butter函数来设计巴特沃斯低通滤波器。下面是一个基于MATLAB实现的巴特沃斯低通滤波器的例子。 首先,我们需要定义一些参数,包括采样频率、通带截止频率、阻带截止频率和通带最大衰减(dB): ```matlab fs = 1000; % 采样频率 fpass = 100; % 通带截止频率 fstop = 150; % 阻带截止频率 Ap = 1; % 通带最大衰减(dB) ``` 然后,我们可以使用butter函数来设计巴特沃斯低通滤波器: ```matlab [n, Wn] = buttord(fpass/(fs/2), fstop/(fs/2), Ap, 60); [b, a] = butter(n, Wn, 'low'); ``` 其中,`buttord`函数用于计算滤波器的阶数和截止频率,`butter`函数用于计算滤波器的系数。 最后,我们可以使用`filter`函数来应用滤波器: ```matlab x = sin(2*pi*50*(0:1/fs:1)); y = filter(b, a, x); ``` 其中,`x`是一个包含50Hz正弦波的信号,`y`是应用了巴特沃斯低通滤波器后的信号。 完整的MATLAB代码如下: ```matlab % 定义参数 fs = 1000; % 采样频率 fpass = 100; % 通带截止频率 fstop = 150; % 阻带截止频率 Ap = 1; % 通带最大衰减(dB) % 设计滤波器 [n, Wn] = buttord(fpass/(fs/2), fstop/(fs/2), Ap, 60); [b, a] = butter(n, Wn, 'low'); % 应用滤波器 x = sin(2*pi*50*(0:1/fs:1)); y = filter(b, a, x); % 绘制结果 subplot(2,1,1); plot((0:length(x)-1)/fs, x); title('原始信号'); xlabel('时间(秒)'); ylabel('幅值'); subplot(2,1,2); plot((0:length(y)-1)/fs, y); title('滤波后信号'); xlabel('时间(秒)'); ylabel('幅值'); ``` 运行此代码将生成一个包含原始信号和滤波后信号的图形。

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 7
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值