matlab凯塞窗低通fir滤波器,基于Matlab的FIR滤波器设计与实现
一、摘要
前面一篇文章介绍了通过FDATool工具箱实现滤波器的设计,见“
二、实验平台
Matlab7.1
三、实验原理
以低通滤波器为例,其常用的设计指标有:
通带边缘频率fp(数字频率为Ωp)
阻带边缘频率fst (数字频率为Ωst)
通带内最大纹波衰减δp=-20log10(1-αp),单位为 dB
阻带最小衰减αs=-20log10(αs),单位为 dB
阻带起伏αs
通带峰值起伏αp
其中,以1、2、3、4条最为常用。5、6条在程序中估算滤波器阶数等参数时会用到。
数字频率 = 模拟频率/采样频率
四、实例分析
例1 用凯塞窗设计一FIR低通滤波器,通带边界频率Ωp=0.3pi,阻带边界频率 Ωs=0.5pi,阻带衰减δs不小于50dB。
方法一:手动计算滤波器阶数N和β值,之后在通过程序设计出滤波器。
第一步:通过过渡带宽度和阻带衰减,计算滤波器的阶数B和β值。
第二步:通过程序设计滤波器。
程序如下:
b = fir1(29,0.4,kaiser(30,4.55));
[h1,w1]=freqz(b,1);
plot(w1/pi,20*log10(abs(h1)));
axis([0,1,-80,10]);
grid;
xlabel('归一化频率/p') ;
ylabel('幅度/dB') ;
波形如下:
方法二:
采用[n,Wn,beta,ftype] = kaiserord(f,a,dev)函数来估计滤波器阶数等,得到凯塞窗滤波器。
这里的函数kaiserord(f,a,dev)或者kaiserord(f,a,dev,fs):
f为对应的频率,fs为采样频率;当f用数字频率表示时,fs则不需要写。
a=[1 0]为由f指定的各个频带上的幅值向量,一般只有0和1表示;a和f长度关系为(2*a的长度)- 2=(f的长度)
devs=[0.05 10^(-2.5)]用于指定各个频带输出滤波器的频率响应与其期望幅值之间的最大输出误差或偏差,长度与a相等,计算公式:
阻带衰减误差=αs,通带衰减误差=αp,可有滤波器指标中的3、4条得到。
fs缺省为2Hz。
程序如下:
fcuts = [0.3 0.5]; %归一化频率omega/pi,这里指通带截止频率、阻带起始频率
mags = [1 0];
devs = [0.05 10^(-2.5)];
[n,Wn,beta,ftype] = kaiserord(fcuts,mags,devs); %计算出凯塞窗N,beta的值
hh = fir1(n,Wn,ftype,kaiser(n+1,beta),'noscale');
freqz(hh);
波形如下:
实际中,一般调用MATLAB信号处理工具箱函数remezord来计算等波纹滤波器阶数N和加权函数W(ω),调用函数remez可进行等波纹滤波器的设计,直接求出滤波器系数。函数remezord中的数组fedge为通带和阻带边界频率,数组mval是两个边界处的幅值,而数组dev是通带和阻带的波动,fs是采样频率单位为Hz。
例2 利用雷米兹交替算法设计等波纹滤波器,设计一个线性相位低通FIR数字滤波器,其指标为:通带边界频率fc=800Hz,阻带边界fr=1000Hz,通带波动 阻带最小衰减At=40dB,采样频率fs=4000Hz。
解:在MATLAB中可以用remezord 和remez两个函数设计
程序如下:
fedge=[800 1000];
mval=[1 0];
dev=[0.0559 0.01];
fs=4000;
[N,fpts,mag,wt]=remezord(fedge,mval,dev,fs);
b=remez(N,fpts,mag,wt);
[h,w]=freqz(b,1,256);
plot(w*2000/pi,20*log10(abs(h)));
grid;
xlabel('频率/Hz') ;
ylabel('幅度/dB');
波形如下:
例3 利用MATLAB编程设计一个数字带通滤波器,指标要求如下:通带边缘频率:Ωp1=0.45pi,Ωp2=0.65pi,通带峰值起伏:δ1<=1[dB]。阻带边缘频率:Ωs1=0.3pi,Ωs2=0.8pi,最小阻带衰减:δ2>=40[dB] 。
方法一:窗函数法
程序如下:
[n,wn,bta,ftype]=kaiserord([0.3 0.45 0.65 0.8],[0 1 0],[0.01 0.1087 0.01]);%用kaiserord函数估计出滤波器阶数n和beta参数
h1=fir1(n,wn,ftype,kaiser(n+1,bta),'noscale');
[hh1,w1]=freqz(h1,1,256);
figure(1)
subplot(2,1,1)
plot(w1/pi,20*log10(abs(hh1)))
grid
xlabel('归一化频率w');ylabel('幅度/db');
subplot(2,1,2)
plot(w1/pi,angle(hh1))
grid
xlabel('归一化频率w');ylabel('相位/rad');
波形如下:
滤波器系数为:
h1 =
Columns 1 through 8
0.0041 0.0055 -0.0091 -0.0018 -0.0056 -0.0000 0.0391 -0.0152
Columns 9 through 16
-0.0381 0.0077 -0.0293 0.0940 0.0907 -0.2630 -0.0517 0.3500
Columns 17 through 24
-0.0517 -0.2630 0.0907 0.0940 -0.0293 0.0077 -0.0381 -0.0152
Columns 25 through 31
0.0391 -0.0000 -0.0056 -0.0018 -0.0091 0.0055 0.0041
如果直接用freqz(h1,1,256),得幅频特性和相频特性曲线:
方法二:等波纹法设计
程序如下:
[n,fpts,mag,wt]=remezord([0.3 0.45 0.65 0.8],[0 1 0],[0.01 0.1087 0.01]);%用remezord函数估算出remez函数要用到的阶n、归一化频带边缘矢量fpts、频带内幅值响应矢量mag及加权矢量w,使remez函数设计出的滤波器满足f、a及dev指定的性能要求。
h2=remez(n,fpts,mag,wt);%设计出等波纹滤波器
[hh2,w2]=freqz(h2,1,256);
figure(2)
subplot(2,1,1)
plot(w2/pi,20*log10(abs(hh2)))
grid
xlabel('归一化频率w');ylabel('幅度/db');
subplot(2,1,2)
plot(w2/pi,angle(hh2))
grid
xlabel('归一化频率w');ylabel('相位/rad');
h2
波形如下:
滤波器系数如下:
h2 =
Columns 1 through 9
-0.0013 0.0092 -0.0255 -0.0642 0.1177 0.0922 -0.2466 -0.0466 0.3116
Columns 10 through 17
-0.0466 -0.2466 0.0922 0.1177 -0.0642 -0.0255 0.0092 -0.0013
如果直接用freqz(h2,1,256);得幅频特性和相频特性曲线:
方法三:采用FDATool工具
这种方法需要事先计算出滤波器的阶数,bate值,然后设置相应参数,最后生成滤波器。
设置界面如下图所示:
将上述圈圈的区域设置好之后,生成滤波器,最后通过analysis菜单可以观察生成的滤波器的各种特性曲线和滤波器系数。这里的滤波器系数跟方法一的一样。
波形如下:
五、结果分析
5.1 滤波器设计总结
FIR滤波器实现一般采用窗函数法和等纹波设计法。窗函数法还包含两个分支,一种是用公式先手动算出N值和其他对应得窗函数参数值,再代入窗函数和fir1实现,一种是用函数*rord估算出N和相应参数再用fir1实现。不过要注意*rord会低估或高估阶次n,可能会使滤波器达不到指定的性能,这时应稍微增加或降低阶次。如果截止频率在0或Nyquist频率附近,或者设定的dev值较大,则得不到正确结果。
滤波器实现形式及特点:由于一般的滤波器在利用窗函数是其通带波纹和阻带波纹不同(一般为第一个阻带波纹最大)因此,在满足第一个阻带衰减旁瓣时,比其频率高的旁瓣,它们的衰减都大大超出要求。而根据阻带衰减与项数的近似关系N = P(δ2)*fs/TW,可得当阻带衰减越大,所需项数越多。
5.2 窗函数法和等波纹设计的不同之处
窗函数设计是通过最小平方积分办法来设计的,即该滤波器的误差为:
即要求最小方法来设计滤波器,这样的滤波器更忠实于理想滤波器(即滤波系数更接近于理想滤波器)。
证明如下:
因此,幅度频谱差值越小,实际滤波器就越接近理想滤波器。
而等波纹滤波器是通过最大加权误差最小化来实现,其误差为:
要求该误差最小来实现滤波器,得出来的滤波系数较窗函数设计相差较远。
以下通过对例3中的h1及h2作比较。
%sigsum是用来对数组各元素进行求和
function y=sigsum(n1,n2,n,x);
y=0;
for i=n1+1-min(n):n2+1-min(n)
y=y+x(i);
end
n=0.001:30.001;
h=2*cos(0.55*pi*(n-15)).*sin(0.175*pi*(n-15))./(pi*(n-15));
delta1=h-h1;
n=0.001:16.001;
h=2*cos(0.55*pi*(n-15)).*sin(0.175*pi*(n-15))./(pi*(n-15));
delta2=h-h2;
y1=sigsum(0,30,[0:30],(abs(delta1).^2))/31;
y2=sigsum(0,16,[0:16],(abs(delta2).^2))/17;
结果如下:
y1 =
1.9099e-004
y2 =
0.0278
由此得到用窗函数实现的滤波系数比用等波纹滤波器系数的每一项更接近于理想滤波器(y1为用窗函数实现的与理想滤波器的差值,y2为用等波纹滤波器实现的与理想滤波器的差值);
对比二者的幅度频谱可知,等波纹滤波器阻带边缘比用窗函数实现的更平滑(理想滤波器为垂直下降的)。
从设计的角度考虑,由于窗函数设计法都是通过已有的窗函数对理想滤波器的改造,因此,可以用手算的办法方便的设计滤波器。
而等波纹滤波器,其实现是通过大量的迭代运算来实现,这样的方法一般只能通过软件来设计。
项数的问题由于等波纹滤波器能较平均的分布误差,因此对于相同的阻带衰减,其所需的滤波系数比窗函数的要少。
5.3 几点说明
1.相频特性曲线形状不同说明
上面第一个图是用角度为单位画出来的,下面的图是用rad单位画出来的。从图形可以观察到在0.3到0.8数字频率间两个图都是严格的线性相位,至于下面的图为什么在这个区间会有跳变是因为rad的区间只有-pi——pi,当相位由-pi继续增加时只能跳到pi而不能大于pi,而角度表示则可以连续增大。
2.调用firl或者reme函数时,用scale(缺省方式)对滤波器进行归一化,即滤波器通带中心频率处的响应幅值为0db。用noscale不对滤波器归一化。
matlab凯塞窗低通fir滤波器,基于Matlab的FIR滤波器设计与实现相关推荐
- 基于matlab的音频信号低通滤波实验
基于matlab的音频信号低通滤波实验 一.实验任务 录制语音信号,查看其时域频谱图. 加入高频噪声信号,并查看其时域频谱图. 设计低通滤波器,将带有高频噪声的语音信号经过低通滤波器,将噪声滤除,查看 ...
- matlab 对图像进行低通滤波,matlab图像处理实现低通滤波
matlab图像处理实现低通滤波 matlab图像处理实现低通滤波 clc; clear all; img=imread('lena.jpg'); img_noise = imnoise(img, ' ...
- matlab凯塞窗设计,转载:用MATLAB设计FIR滤波器的方法解析
介绍了利用MATLAB信号处理工具箱进行FIR滤波器设计的三种方法:程序设计法.FDATool设计法和SPTool设计法,给出了详细的设计步骤,并将设计的滤波器应用到一个混和正弦波信号,以验证滤波器的 ...
- MATLAB图像处理--高斯低通滤波、高斯高通滤波(代码及示例)
本文目录 常用的滤波器 高斯低通滤波器 matlab代码 高斯高通滤波器 Butterworth 低通滤波器 常用的滤波器 常用的滤波器有以下几种: 频域滤波器 频域平滑滤波器 理想低通滤波器 巴特沃 ...
- Matlab频域高/低通滤波
建议参考书籍:数字图像处理_第三版 冈萨雷斯 写在前面: 对于给定的低通滤波器的函数表达式,可以得到高通滤波器的函数表达式: 理想高/低通滤波器 理想高通 一个二维理想高通滤波器(IHPF)定义为: ...
- matlab测量直流母线上的电压,基于模糊控制的有源滤波器直流母线电压控制
摘要针对传统直流侧电压控制方法存在超调量和静差较大的问题,设计了混合有源电力滤波器控制系统,建立了混合有源电力滤波器的系统模型,并针对直流侧电压控制设计了一套双模糊控制器,与传统直流侧电压控制方法相比 ...
- 设计低通巴特沃斯数字滤波器(matlab)去除ECG基线漂移实例
matlab代码************************** M=load ('datafile'); ecg=M(1:15000,1); 读取数据 fs=512; 采样率 fmaxd_1=5 ...
- matlab 有源电力滤波器,基于MATLAB的有源电力滤波器研究
第 3 卷第 4 期 南 阳 理 工 学 院 学 报 Vol. 3 No. 4 2 0 1 1 年 7 月 JOURNAL OF NANYANG INSTITUTE OF TECHNOLOGY Jul ...
- MATLAB高通滤波与低通滤波GUI实现代码
1.需搭建的GUI界面 两个axes,一个用于导入原始图像,一个用于输出处理后的头像,并且记好他们的Tag,在代码里会用到 2.导入图片功能代码实现 给一个String为导入图片的button添加ca ...
最新文章
- Winform程序窗体间的跳转
- Package xxx is not available, but is referred to by another package
- EPUB CFI 规范
- ajax点击更改div,jquery ajax双击div可直接修改div中的内容
- ML之XGBoost:XGBoost参数调优的优秀外文翻译—《XGBoost中的参数调优完整指南(带python中的代码)》(四)
- #翻译# 关于 Java 性能方面的 9 个谬论
- android log丢失(一)使用logd丢失log原理
- java的volatile是什么意思
- masm汇编语言堆栈段定义了却提示无堆栈段
- 安卓系统曝root权限安全漏洞 已有软件利用
- PAT 乙级 1037. 在霍格沃茨找零钱(20)Java版
- Akka和Actor一起工作的消息《ten》译
- php基础-基本语法
- 3001.Linux维护和管理培训.实操考核A
- 以锅炉安全为例,台账管理在安全建设中的作用
- 小米盒子 android 投影仪,小米盒子连接投影仪方法
- Win10显示无网络(无Internet连接)但能上网的解决办法
- Python爬虫使用lxml模块爬取豆瓣读书排行榜并分析
- ggplot去掉背景网格线和更改为白色背景
- [解决方案]罗技POWERPLAY鼠标垫无法连接G703/G903/G403等鼠标
热门文章
- 转载好文 php读取大文件
- anemometer mysql 500_Anemometer MySQL 慢查询日志监控平台
- python unicode转字符串_Python: 在Unicode和普通字符串之间转换
- 硬件笔记(2)---- 贴片电容材质NPO与X7R、X5R、Y5V、Z5U的区别
- 矩阵分解 SVD分解
- python中quad_python scipy integrate.quad用法及代码示例
- 中兴核心网服务器笔记本电脑,基于统一硬件平台的中兴通讯核心网技术介绍
- 【2017.10.08 智能驾驶/汽车电子】汽车电控常用英文缩写及功能之一
- opencv计算指定区域黑白像素占比
- Oracle用户被锁定问题