音乐信号滤波去噪——使用flattopwin窗设计的FIR滤波器
时间：2021.03.08 创作：欧阳与 学生姓名： 指导老师：黄红兵
摘 要 本次课程设计是使用Flattopwin窗设计FIR滤波器对音乐信号进行滤波去噪。通过MATLAB软件，运用窗函数法来设计滤波器。从网上下载一段满足要求的音乐，为它加入噪声信号，观察加噪前后的频谱，采用窗函数设计法，给定相应的技术指标，设计一个满足要求的滤波器，对音乐信号进行滤波去噪处理。比较原始音乐信号与滤波后的时域波形图，频谱图，回放滤波后的音乐信号，可听见滤波后的音乐信号与原始音乐信号无大致差别，成功的实现了滤波达到了设计要求。
关键词 MATLAB；滤波去噪；FIR滤波器；Flattopwin窗；
1 引 言
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
本次课程设计是通过窗函数法设计一个Flattopin的FIR滤波器对有噪声的音乐信号进行滤波去噪。在课程设计中，了解FIR滤波器的性能与原理，也了解他的设计方法和步骤。掌握了用MATLAB语言设计滤波器，通过观察音乐信号滤波前后的时域信号以及频谱更加具体的了解了滤波器的作用。 1.1 课程设计目的
通过利用MATLAB 软件来利用Flattopin设计FIR滤波器对音乐信号滤波去噪。使得我们更加熟悉MATLAB的语言环境，更加熟悉MATLAB语言的编程规则。并且在课程设计中通过观察滤波器的幅度，相位图对Flattowin有了更加深刻地了解。也在窗函数的设计过程中，对滤波器的性能，功能以及设计方法有着更具体的了解和体验。通过本次课程设计，增强了我们独立解决问题的能力，提高了自己的动手能力。 1.2 课程设计要求
从网上下载一段.wav格式的音乐，绘制观察时域波形及频谱图。对音乐信号加入噪声干扰，根据Flattopwin的性能指标合理设计FIR滤波器，再用滤波
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
器对加入干扰的音乐信号进行滤波去噪，比较滤波前后的频谱图并进行分析。再回放语音信号对比原语音信号，查看滤波器是否对语音信号进行了滤波去噪。 1.3课程设计平台
本次课程设计通过MATLAB实现，MATLA是美国MathWorks公司出品的商业数学软件，用于算法开发、数据可视化、数据分析以及数值计算的高级技术计算语言和交互式环境，主要包括MATLAB和Simulink两大部分。
MATLAB是matrix&laboratory两个词的组合，意为矩阵工厂（矩阵实验室）。是由美国mathworks公司发布的主要面对科学计算、可视化以及交互式程序设计的高科技计算环境。它将数值分析、矩阵计算、科学数据可视化以及非线性动态系统的建模和仿真等诸多强大功能集成在一个易于使用的视窗环境中，为科学研究、工程设计以及必须进行有效数值计算的众多科学领域提供了一种全面的解决方案，并在很大程度上摆脱了传统非交互式程序设计语言（如C、Fortran）的编辑模式，代表了当今国际科学计算软件的先进水平[1]。
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
MATLAB软件包括五大通用功能：数值计算功能（Nemeric）；符号运算功能（Symbolic）；数据可视化功能（Graphic）；数据图形文字统一处理功能（Notebook）和建模仿真可视化功能（Simulink）。其中，符号运算功能的实现是通过请求MAPLE内核计算并将结果返回到MATLAB命令窗口。该软件有三大特点：一是功能强大；二是界面友善、语言自然；三是开放性强。MATLAB在线性代数、矩阵分析、数值及优化、数理统计和随机信号分析、电路与系统、系统动力学、信号和图像处理、控制理论分析和系统设计、过程控制、建模和仿真、通信系统、以及财政金融等众多领域的理论研究和工程设计中得到了广泛应用。
MATLAB在信号与系统中的应用主要包括符号运算和数值计算仿真分析。由于信号与系统课程的许多内容都是基于公式演算，而MATLAB借助符号数学工具箱提供的符号运算功能能基本满足信号与系统课程的需求。例如，解微分方程、傅里叶正反变换、拉普拉斯正反变换、z正反变换等。MATLAB在信号与系统中的另一主要应用是数值计算与仿真分析，主要包括函数波形绘函数运算、冲激响应与阶跃响应仿真分
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
析、信号的时域分析、信号的频谱分析、系统的S域分析、零极点图绘制等内容。数值计算仿真分析可以帮助学生更深入理解信号与系统的理论知识，并为将来使用MATLAB进行信号处理领域的各种分析和实际应用打下基础[2]。
2 基本理论
2.1 FIR滤波器
FIR滤波器：有限长单位冲激响应滤波器，是数字信号系统中最基本的元件，它可以在保证任意幅频特性的同时具有严格的线性相频特性，同时其单位抽样响应是有限长的，因而滤波器是稳定的系统。因此，FIR滤波器在通信、图像处理、模式识别等领域都有着广泛的应用。
FIR滤波器的基本结构：横截型，级联型，频率抽样型，快速卷积型，在硬件上一般通过集成电路，DSP芯片，可编程逻辑器件，FPGA/CPLD来实现[3]。
FIR数字滤波器的特点：
优点：(1) 很容易获得严格的线性相位，避免所处理的信号产生相位失真；
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
(2) 极点全部在原点，无稳定性问题； (3) 任何一个非因果的有限长序列，总可以通
过一定的延时，转化为因果序列，所以因果性总是可以满足；
(4) 无反馈运算，运算误差小。
缺点：(1) 因为无极点，要获得好的过渡带特性，
要以较高的结束为代价；
(2) 无法利用模拟滤波器的设计结果，一般无
解析设计公式，要借助辅助程序设计完成。
一般来说滤波器的设计在于寻找一个频率响应去逼近，逼近的方法有三种窗函数设计法（时域逼近）；频率采样设计法(频域逼近)；最优化设计法（等波纹逼近）。本次课程设计采用窗函数设计法[4]。
2.2 窗口设计法
窗函数设计法是从单位脉冲响应序列着手，使得逼近理想单位脉冲响应序列
(2.1)
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
窗函数设计FIR滤波器的步骤：
(1) 根据技术要求确定滤波器的频率响应确定其对应的
单位脉冲响应。
(2.2)
(2) 根据对过渡带及阻带衰减指标的要求，选择窗函数
形式，并估计窗口长度N。
(3) 计算滤波器的单位取样响应。
(2.2)
(4) 计算滤波器的频率响应。
(2.3)
改变窗函数的形状，可改善滤波器的特性，窗函数有许多种，但要满足以下两点：
(1) 窗谱主瓣宽度要窄，已获得较陡的过渡带； (2) 相对于主瓣宽度，旁瓣要尽可能小，是能量尽量集中在主瓣中，这样就可以减小肩峰和余振，以提高阻带衰减和通带平稳性。
但实际上这两点不能兼得，一般总是通过增加主瓣
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
宽度来换取对旁瓣的抑制[5]。 常见的窗函数性能表如表2.1所示：
表2.1常见窗函数性能表
滤波器 名称 过渡带宽 矩形 巴特利特 汉宁 汉明 布莱克曼 1.8π/M 6.1π/M 6.2π/M 6.6π/M 11π/M 最小阻带衰减 滤波器 名称 过渡带宽 PARZENWIN FLATTOPWIN GAUSSWIN BARTHANNWIN BLACKMANHARRIS 6.6π/M 19.6π/M 5.8π/M 3.6π/M 16.1π/M 56dB 108dB 60dB 40dB 109dB 最小阻带衰减 21dB 25dB 44dB 51dB 74dB BOHMANWIN NUTTALLWIN 5.8π/M 15.4π/M 51.5dB 108dB CHEBWIN TUKEYWIN 15.2π/M 2.4π/M 113dB 22dB 2.3 FLATTOPWIN窗
w=Flattopwin (L) 返回L-点Flattopwin窗口中列向量。Flattopwin窗的滤波器的过渡带宽为19.6π/M，最小阻带衰减108db。
定义式：
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
2n4n6n8n)a2cos()a3cos()a4cos()NNNN(n)a0a1cos(
(2.4)
时间波形和幅度谱：
图2.1 Flattopwin窗的时间波形 图2.2 Flattopwin窗的幅度谱
2.4滤波器的结构 (1) 横截型
差分方程： (2.5)
图2.3 FIR滤波器横截型结构
(2)级联型
将差分方程分解为实系数二阶因式乘积形式： (2.6)
图2.4 FIR滤波器级联型结构
当N为偶数时，其中有一个=0（N-1个零点）。 级联型每个基本节控制一对零点，便于控制滤波器的传输零点；系数比直接型多，所需乘法运算多。
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
(3)线性相位型
FIR 的重要特点是可设计成具有严格线性相位的滤波器，此时满足偶对称或奇对称。
偶对称时：
(2.7)
N为奇对称时：
(2.8)
图2.5 N为偶数时的线性相位结构FIR滤波器 图2.6 奇为偶数时的线性相位结构FIR滤波器
(5) 频率采样型
H(z)由两部分组成： 第一部分（FIR部分）：
(2.9)
第二部分（IIR部分）：
(2.10)
图2.7 FIR滤波器的频率采样型结构
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
这一结构的最大的特点是它的系数H(k)直接就是滤波器在处的响应，因此控制滤波器响应很直接[6]。 3 设计步骤 3.1 设计流程图
根据设计的要求，从网上下载一段音乐信号，在通过系统自带的软件对音乐信号进行格式转换以及剪裁，然后对音乐信号加入噪声干扰，再利用Flattopwin窗设计合理的FIR滤波器。最后用滤波器对干扰后的音乐信号进行滤波去噪。具体设计流程图如下图3.1所示：
图3.1设计
流程图
3.2 录制音乐信号
从电脑上下载一段格式为.wav的纯音乐，并命名为yinyue.wav，通过录音机将原始音乐信号的格式更改为PCM编码，采样率为8KHz，8位，单声道，7KB/s并将其裁剪为10s的文件，调用wavread函数读取信号的参数，绘制出音乐信号的时域波形。
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
图3.2用录音机将音乐进行裁剪 图3.3用录音机将音乐格式进行转化
图3.4 原始音乐信号时域波形
图
绘制音乐信号的时域波形后，对音乐信号进行快速傅里叶变换，得到信号的频谱特性，绘制出音乐信号的频谱图，向音乐信号加入噪声，并对加入噪声的信号进行快速傅里叶变换，绘制出加入噪声后的时域波形，以及频谱图，将加入噪声后的信号时域波形图以及频谱图与未加入之前的波形图和频谱图进行对比。
>>x=x'; y=x+0.1*sin(fn*2*pi*t); %在原信号中加入噪声
>>X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱
>>X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分 >>deltaf=fs/N; % 计算频谱的谱线间隔 >>f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
图3.5 干扰前后信号时域图与频谱
图的比较 3.3 滤波器设计
根据加入噪声干扰后的频谱图，给出滤波器的性能指标，在通过窗函数法设计得到自己的滤波器，得到数字滤波器的参数a，b。其中a为系统函数的分母系数，b为系统函数的分子系数。在通过调用函数即可求得滤波器的频率响应。在计算滤波器的实际AS以及RP查看滤波器是否满足设计指标。
>>fpd=1825;fsd=2075;fsu=2325;fpu=2575;Rp=1;As=20; % 带阻滤波器设计指标
>>RP=-min(db([1:wpd/dw+1,ceil(wpu/dw)+1:end]))
%计算滤波器的
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
指标检查是否达标
>>AS=-max(db(wsd/dw+1:wsu/dw+1))
通过计算可以得出RP=0.7939，AS=20.3545满足指标要求，则滤波器符合要求。
图3.6设计的滤波器图
形
3.4 信号滤波处理
滤波器设计完成后，在MATLAB平台上用函数filter对加有噪声的音乐信号进行滤波，并求得滤波后的频谱。绘制滤波前后音乐信号的时域波形图和频谱图对比图，以及原信号的时域波形图及频谱图。
>>y_fil=filter(h_bs,1,y);%用设计好的滤波器对y进
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
行滤波
图3.7 滤波前后音乐信号的时域波形图
和频谱图对比图 3.5 结果分析
观察滤波前后的频谱图，可明显看看到干扰噪声的频谱被滤除，观察滤波后的时域波形图与原信号的波形图基本一致，在通过sound()函数对经过Flattopwin窗设计的FIR滤波器之后的音乐信号进行回访可以听出滤波之后的信号跟原始信号一样清晰，完全滤除掉了噪声的干扰。
>> sound (y_fil,fs,bits); %播放滤波后的音乐信号 >>sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放
所得结果证明了用Flattopwin窗设计的FIR滤波器对音乐信号进行滤波去噪是成功的。 3.6 滤波器结构设计
在绘制滤波器结构时由于我选择的级联型，因此要计算相关参数，可以通过相关函数来进行计算
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
>>[C,B,A]=DIR2CAS(h_bs,1); %计算级联型滤波器的参数
然后在Viso中绘制出FIR滤波器的级联型结构即可。
图3.6 FIR滤波器的级联型结构 4出现的问题及解决方法
在本次课程设计中我遇到的问题如下：
1、绘制加入噪声干扰的频谱图是不知如何调控纵坐
标，导致频谱图十分的不理想
2、在达到设计指标且完美的滤除噪声干扰信号时，滤
波器的阶数总是过大，而将阶数降低时，滤波器总是不能比较完美的滤除噪声
3、绘制滤波器结构时选择的是级联型，但是却找不到，
计算级联型参数的函数
针对以上问题，相应的解决方案如下：
1、通过向同学询问了解到有函数可以用来调控横纵
坐标
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
2、通过不断设置滤波器参数值最终设计出理想的滤
波器，可在阶数较低的情况下完美
的滤出噪声，也达到了设计指标。
3、在查阅历届的
PPT时，发现了用于计算直接型转
换级联型滤波器参数值的函数。
5 结束语
本次课程设计，我的课题是利用Flattopwin窗函数设计FIR滤波器对音乐信号滤波去噪。在经历了一个寒假的放松，脑袋里那点不多的滤波器的只是，早已魂飞天外。但是在实习的过程中，通过温习课本以及通过频谱图，对滤波器的性能以及功能有了更加深刻的了解。
每次课程设计都是一次历练，是对知识的查漏补缺，是一次给自己深刻反思的机会。每一次实习，都可以意识到自己所欠缺的，意识到在实践面前理论是如何的单薄。知识的创造永远是为了人类的发展而服务的，只存在于理论的知识是无法创造它所应有的价值，而进步却又来源于实践和人民的需求，我希望我们用时间去学习的东西能发挥它应有的价值，甚至于
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
超越。每一次实践都可以学到许多课本上所没有的只是，也了解到以前所没有注意到的方面，而在时间时，不断锻炼的对立思考的能力更是能让我们终生受益。
相对与知识，意识总是显得尤为重要。而在实践中，重温的不仅仅是知识更是，我们在于以后是的意识，去独立的思考问题，去查阅书籍，去探寻奥妙。而这些是我们平时所不能学到的。
我总是希望自己能从理论走向时间，从书本走向课本，希望自己能用自己所学的为这个发展的社会的建设献上自己的一点微薄之力。
参考文献
[1]张志涌．精通MATLAB [M]．6.5版．北京：北京航空航天大学出版社，2003．
[2]约翰·Ｇ·普罗克斯．数字信号处理[M]．西安：西安交通大学出版社，2009．
[3]张小虹．信号系统与数字信号处理[M]．第1版．西安：西安电子科技出版社，2002．
[4]谢德芳．数字信号处理[M]．北京：科学出版社，2005．
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
[5]郝建军．数字通信[M]．第2版．北京：北京邮电大学出版社，2010．
[6]张威．MATLAB 基础与编程入门[M]．西安：西安电子科技大学出版社，2010．
附录一： 音乐信号滤波去噪——使用FLARTOPWIN滤波器
%程序功能：在Matlab中，用窗口设计法设计FIR滤波器
%程序作者：蒋敏
%最后修改日期：2016-3-1
程序一------------------------(在音乐信号中加入噪声干扰)------------------------------------------
[x,fs,bits]=wavread('e:\\yinyue.wav'); % 输入参数为文件的全路径和文件名，输出的第一个参数是每个样本的值，fs是生成该波形文件时的采样率，bits是波形文件每样本的编码位数。
sound(x,fs,bits); % 按指定的采样率和每样本编码位数回放
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
N=length(x); % 计算信号x的长度 fn=2200; % 单频噪声频率
t=0:1/fs:(N-1)/fs; % 计算时间范围，样本数除以采样频率
x=x'; y=x+0.1*sin(fn*2*pi*t); %在原信号中加入噪声 sound(y,fs,bits); % 应该可以明显听出有尖锐的单频啸叫声
程序二-----------------------(绘制加噪前后音乐信号的时间波形和频谱)------------------------- X=abs(fft(x)); Y=abs(fft(y)); % 对原始信号和加噪信号进行fft变换,取幅度谱
X=X(1:N/2); Y=Y(1:N/2); % 截取前半部分 deltaf=fs/N; % 计算频谱的谱线间隔 f=0:deltaf:fs/2-deltaf; % 计算频谱频率范围 figure(); %绘制加噪前后音乐信号的时间波形和频谱
subplot(221);plot(t,x);title('原始音乐信号');xlabel('时间
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
（单位 ： s ）');ylabel('幅度'); grid on;
subplot(223);plot(t,y);title('加入干扰信号的音乐信号');xlabel('时间（单位 ： s ）');ylabel('幅度');grid on; subplot(222);plot(f,X);title('音乐信号幅度谱');xlabel('频率（单位 ： Hz ）');ylabel('幅度谱');grid on;
subplot(224);plot(f,Y);axis([0,4000,0,400]); title('加入干扰信号的音乐信号幅度谱');xlabel('频率（单位 ： Hz ）');ylabel('幅度谱');grid on;
程序三--------------------------------(设计滤波器)----------------------------------------------------
fpd=1825;fsd=2075;fsu=2325;fpu=2575;Rp=1;As=20; % 带阻滤波器设计指标
fcd=(fpd+fsd)/2;fcu=(fpu+fsu)/2;df=min((fsd-fpd),(fpu-fs
欧阳与创编 2021.03.08
u)); %
欧阳与创编 2021.03.08
wcd=fcd/fs*2*pi;wcu=fcu/fs*2*pi;dw=df/fs*2*pi; 欧阳与创编 2021.03.08
% 将
Hz
计算上下边带中心频率，和频率间隔
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
为单位的模拟频率换算为rad为单位的数字频
欧阳与创编 2021.03.08
率
wsd=fsd/fs*2*pi; wsu=fsu/fs*2*pi; wpd=fpd/fs*2*pi; wpu=fpu/fs*2*pi;
M=ceil(8.5*pi/dw)+1 ; % 计算窗设计该滤波器时需要的阶数
n=0:M-1; % 定义时间范围 w_par=(flattopwin (M)); % 产生M阶的Flattopwin窗
hd_bs=IDEAL_LP(wcd,M)+IDEAL_LP(pi,M)-IDEAL_L
欧阳与创编 2021.03.08
P(wcu,M
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
); % 调用自编函数计算理想带阻滤波器的
欧阳与创编 2021.03.08
脉冲响应
h_bs=w_par'.*hd_bs; % 用窗口法计算实际滤波器脉冲响应
[db,mag,pha,grd,w]=FREQZ_M(h_bs,1); % 调用自编函数计算滤波器的频率特性 dw=2*pi/1000;
RP=-min(db([1:wpd/dw+1,ceil(wpu/dw)+1:end])) %计算滤波器的指标检查是否达标 AS=-max(db(wsd/dw+1:wsu/dw+1)) figure() ;
subplot(2,2,1);plot(w/pi,db);title('滤波器幅度响应图');xlabel('w/pi');ylabel('db'); axis([0,1,-70,10]);
line([0,1],[-As,-As],'color','r','linestyle','--','LineWidth',2); line([0,1],[-Rp,-Rp],'color','r','linestyle','--','LineWidth',2);
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
line([wsu/pi,wsu/pi],[-70,10],'color','r','linestyle','--','LineWidth',2);
line([wsd/pi,wsd/pi],[-70,10],'color','r','linestyle','--','LineWidth',2);
subplot(2,2,2);plot(w/pi,mag);title('滤波器幅度响应图');xlabel('w/pi');ylabel('幅度mag'); grid on;
subplot(2,2,3);plot(w/pi,pha);title('滤波器相位响应图');xlabel('w/pi');ylabel('相位pha'); grid on;
subplot(2,2,4);stem(n,h_bs);title('滤波器脉冲响应图');xlabel('n');ylabel('h(n)');grid on; 程序四-------------------------------------(滤波)------------------------------------------------------
y_fil=filter(h_bs,1,y); % 用设计好的滤波器对y进行滤波
Y_fil=abs(fft(y_fil));Y_fil=Y_fil(1:N/2); % 计算频谱取前一半 figure();
subplot(321);plot(t,x);title('原始音乐信号');xlabel('时间
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
（单位 ： s ）');ylabel('幅度'); grid on;
subplot(323);plot(t,y);title('加入干扰信号的音乐信号');xlabel('时间（单位 ： s ）');ylabel('幅度');grid on; subplot(322);plot(f,X);title('音乐信号幅度谱');xlabel('频率（单位 ： Hz ）');ylabel('幅度谱');grid on;
subplot(324);plot(f,Y);axis([0,4000,0,400]); title('加入干扰信号的音乐信号幅度谱');xlabel('频率（单位 ： Hz ）');ylabel('幅度谱');grid on;
subplot(325);plot(t,y_fil); title('滤波后的音乐信号');xlabel('时间（单位 ： s ）');ylabel('幅度');grid on; subplot(326);plot(f,Y_fil);title('滤波后的音乐信号幅度谱');xlabel('频率（单位 ： Hz ）');ylabel('幅度谱'); grid on;
sound (y_fil,fs,bits);
[C,B,A]=DIR2CAS(h_bs,1); %计算级联型滤波器参数 附录二：
函数FREQZ_M.M定义：
欧阳与创编
2021.03.08
欧阳与创编 2021.03.08
function [db,mag,pha,grd,w] = freqz_m(b,a); % freqz 子程序的改进版本 % ------------------------------------ % [db,mag,pha,grd,w] = freqz_m(b,a);
% db = [0 到pi弧度]区间内的相对振幅(db) % mag = [0 到pi弧度]区间内的绝对振幅 % pha = [0 到pi弧度]区间内的相位响应 % grd = [0 到pi弧度]区间内的群迟延
% w = [0 到pi弧度]区间内的501个频率样本向量 % b = Ha(z)的分子多项式系数(对FIR b=h) % a = Ha(z)的分母多项式系数(对 FIR: a=[1]) %
[H,w] = freqz(b,a,1000,'whole');
H = (H(1:1:501))'; w = (w(1:1:501))'; mag = abs(H);
db = 20*log10((mag+eps)/max(mag));
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
pha = angle(H);
% pha = unwrap(angle(H)); grd = grpdelay(b,a,w); % grd = diff(pha); % grd = [grd(1) grd];
% grd = [0 grd(1:1:500); grd; grd(2:1:501) 0]; % grd = median(grd)*500/pi; 附录三：
函数IDEAL_LP.M定义： function hd = ideal_lp(wc,M); % 理想低通滤波器计算 % -------------------------------- % [hd] = ideal_lp(wc,M)
% hd = 0 to M-1之间的理想脉冲响应 % wc = 截止频率(弧度) % M = 理想滤波器的长度
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
%
alpha = (M-1)/2; n = [0:1:(M-1)]; m = n - alpha + eps; hd = sin(wc*m) ./ (pi*m); 附录3
function [b0,B,A] = dir2cas(b,a);
% 直接型到级联型的型式转换(复数对型) % --------------------------------------------------------- % [b0,B,A] = dir2cas(b,a)
% b = 直接型的分子多项式系数 % a = 直接型的分母多项式系数 % b0 = 增益系数
% B = 包含各bk的K乘3维实系数矩阵 % A = 包含各ak的K乘3维实系数矩阵 % compute gain coefficient b0
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
b0 = b(1); b = b/b0; a0 = a(1); a = a/a0; b0 = b0/a0;
M = length(b); N = length(a); if N > M
b = [b zeros(1,N-M)]; elseif M > N
a = [a zeros(1,M-N)]; N = M; else NM = 0; end
K = floor(N/2); B = zeros(K,3); A = zeros(K,3); if K*2 == N; b = [b 0]; a = [a 0]; end
欧阳与创编 2021.03.08
欧阳与创编 2021.03.08
broots = cplxpair(roots(b)); aroots = cplxpair(roots(a)); for i=1:2:2*K
Brow = broots(i:1:i+1,:); Brow = real(poly(Brow)); B(fix((i+1)/2),:) = Brow; Arow = aroots(i:1:i+1,:); Arow = real(poly(Arow)); A(fix((i+1)/2),:) = Arow; end
时间：2021.03.08 创作：欧阳与 欧阳与创编 2021.03.08
因篇幅问题不能全部显示，请点此查看更多更全内容