如何利用matlab处理音频信号_百度文库
您的浏览器Javascript被禁用，需开启后体验完整功能，
享专业文档下载特权
&赠共享文档下载特权
&10W篇文档免费专享
&每天抽奖多种福利
两大类热门资源免费畅读
续费一年阅读会员，立省24元！
如何利用matlab处理音频信号
阅读已结束，下载本文需要
想免费下载本文？
定制HR最喜欢的简历
下载文档到电脑，同时保存到云知识，更方便管理
加入VIP
还剩16页未读，
定制HR最喜欢的简历
你可能喜欢matlab中如何生成一个时间数组，比如我想生成9：15：00到11：29：59的序列！_百度知道
matlab中如何生成一个时间数组，比如我想生成9：15：00到11：29：59的序列！
我有更好的答案
是不是要这样的？09:15:0009:15:0109:15:0209:15:0309:15:0409:15:05&（中间省略很多）11:29:5711:29:5811:29:59&代码：s=datenum('9:15:00');e=datenum('11:29:59');d=1/3600/24;T=datestr((s:d:e)','HH:MM:ss')
采纳率：90%
来自团队：
为您推荐：
其他类似问题
matlab的相关知识
换一换
回答问题，赢新手礼包
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。只需一步，快速开始
扫一扫，访问微社区
查看: 5263|回复: 5|关注: 0
用matlab产生的声音信号
<h1 style="color:# 麦片财富积分
新手, 积分 19, 距离下一级还需 31 积分
关注者： 1
Fs = 44100; % 采样频率
T = 4;& &&&% 时间长度
n = Fs*T;&&% 采样点数
f = 500;& &% 声音频率
y = sin(2*pi*f*T*linspace(0,1,n+1));
sound(y,Fs)复制代码例如上面的代码，我产生了500HZ频率的sin声音信号，那么，在y=sin()这里面，我不明白为什么会需要乘T，另外linspace是做什么用的？仅仅是0和1之间取那么多点吗？&&代表了什么？&&
如果我要让一个矩形波（周期性）产生声音，是不是也是只要让fs满足奈奎斯特定理后，再0和1之间取点就可以了？
论坛优秀回答者
帖子最佳答案
关注者： 273
T*linspace(...)是你的各个采样点时间
<h1 style="color:# 麦片财富积分
关注者： 1
T*linspace(...)是你的各个采样点时间
奥，如果我要取别的函数，也可以这样定义吗
论坛优秀回答者
帖子最佳答案
关注者： 273
奥，如果我要取别的函数，也可以这样定义吗
不懂你的问题，什么函数
<h1 style="color:# 麦片财富积分
关注者： 1
不懂你的问题，什么函数
矩形波~ 那种。&&另外 我生成一段0.01秒的声音，但是我要让它延长至1秒，如果不改变时间t，还能用什么函数？&&采样率是44.1KHZ
论坛优秀回答者
帖子最佳答案
关注者： 273
|此回复为最佳答案
矩形波~ 那种。&&另外 我生成一段0.01秒的声音，但是我要让它延长至1秒，如果不改变时间t，还能用什么函 ...
sin的本质就是sin(2*pi*f*t)。你只要保证你的函数满足这个条件就可以了。至于矩形波，可以看看rectpuls和pulstran
站长推荐 /3
筑起功能安全的堡垒 - 基于模型设计的软件开发
MATLAB中文论坛是全球最大的 MATLAB & Simulink 中文社区。用户免费注册会员后，即可下载代码，讨论问题，请教资深用户及结识书籍作者。立即注册加入我们吧！
MATLAB官方社交平台
MATLAB中文论坛微社区MATLAB中导入的声音如何看它的频率节数_百度知道
MATLAB中导入的声音如何看它的频率节数
我有更好的答案
1）用mp3录音，生成文件cricket.wav，把该文件放到matlab文件夹里面（就是你打开matlab后中间顶部的地址）。2）使用如下程序，做波形显示以及fft变换。[y,Fs,bits]=wavread(&#39;cricket.wav&#39;);%读出信号，采样率和采样位数。y=y(:,1);%我这里假设你的声音是双声道，我只取单声道作分析，如果你想分析另外一个声道，请改成y=y(:,2)sigLength=length(y);Y = fft(y,sigLength);Pyy = Y.* conj(Y) / sigLhalflength=floor(sigLength/2);f=Fs*(0:halflength)/sigLplot(f,Pyy(1:halflength+1));xlabel(&#39;Frequency(Hz)&#39;);t=(0:sigLength-1)/Fs;plot(t,y);xlabel(&#39;Time(s)&#39;);3）频率看频谱就有了，声音间隔看声音波形，周期看声音波形。
采纳率：75%
来自团队：
为您推荐：
其他类似问题
matlab的相关知识
&#xe675;换一换
回答问题，赢新手礼包&#xe6b9;
个人、企业类
违法有害信息,请在下方选择后提交
色情、暴力
我们会通过消息、邮箱等方式尽快将举报结果通知您。当前位置： >>
MATLAB实现
1.9 1.9.1离散信号和系统分析的 MATLAB 实现 利用 MATLAB 产生离散信号用 MATLAB 表示一离散序列 x[k]时，可用两个向量来表示。其中 一个向量表示自变量 k 的取值范围，另一个向量表示序列 x[k]的值。 例如序列 x[k]={2,1,1,-1,3,0,2}可用 MATLAB 表示为 K=-2:4;x=[2,1,1,-1,3,0,2] 可用 stem(k,f)画出序列波形。 当序列是从 k=0 开始时，可以只用一个向量 x 来表示序列。由于 计算机内寸的限制，MATLAB 无法表示一个无穷长的序列。 例1-38 利用 MATLAB 计算单位脉冲序 δ [k ? 2]在 ? 4 ≤ k ≤ 4 范围内各点的取值。 解： %progran 1_1 Ks=-4;ke=4;n=2; K=[ks:ke]; X=[(k-n)==0]; Stem(k,x):xlabel(‘k’); 程序产生的序列波形如图 1-49 所示。 例1-39 利用 MATLAB 画出信号 X[k]=10sin(0.02 πk )+n[k],0 ≤ k ≤ 100产生单位脉冲序列的波形。其中 n[k]表示为均值为 0 方差为 1 的 Gauss 分布随机信号。 解：1MALAB 提供了两个产生（伪）随机序列的函数。Rand(1,N)产生 1 行 N 列的[0，1]均匀分布随机数。Randn(1,N)产生 1 行 N 列均值为 0 方差为 1 的 Gauss 分布随机数。 %program 1_2 N=100;k=0:N; X=10*sin(0.02*pi*k)+randn(1,N+1); Plot(k,x); Xlabel(‘k’); Ylabel(‘x[k]’); 程序产生序列如图 1-50 所示。 1.9.2 离散卷积的计算 产生受噪声干扰的正弦信号离散卷积是数字信号处理中的一个基本运算，MTLAB 提供的计 算两个离散序列卷积的函数是 conv，其调用方式为 y=conv(x,h) 其中调用参数 x,h 为卷积运算所需的两个序列， 返回值 y 是卷积结果。 MATLAB 函数 conv 的返回值 y 中只有卷积的结果，没有 y 的取 值范围。由离散序列卷积的性质可知，当序列 x 和 h 的起始点都为 k=0 时， 的取值范围为 k=0 至 length(x)+length(h)-2。 y 当序列 x 或 （和） h 的起始点不是 k=0 时，由例 1-3 知，y 的非零值范围可根据例 1-4 的结论进行计算。 例 1-40 试用 MATLAB 函数 conv 计算例 1-2 中序列的卷积。 解：2program 1_3 x=[-0.5,0,0.5,1]; kx=-1:2; h=[1,1,1]; kh=-2:0; y=conv(x,h);计算离散卷积 %序列 x 的值%序列 x 的取值范围%计算卷积k=kx(1)+kh(1):kx(end)+kh(end); %计算 y 的取值范围 stem(k,y); xlabel(‘k’);ylabel(‘y’); 程序的运行结果如图 1-51 所示。 1.9.3 离散 LTI 系统响应 MATLAB 求解 许多离散 LTI 都可用如下的线性常系数的差分方程描述∑ any[k ? n] = ∑ bnx[k ? n]n =0 n=0NM(1-151)其中 x[k]、y[k]分别系统的输入和输出。在已知差分方程的 N 个 初始状态 y[k], ? N ≤ k ≤ ?1 和输入 x[k],就可由下式迭代计算出系统的 输出。 y[k]=- ∑ (an / a0) y[k ? n] + ∑ (bn / b0) x[k ? n]n =1 n=0 N M利用 MATLAB 提供的 filter 函数， 可方便地计算出上述差分方程 的零状态响应。filter 函数调用形式为 y=filter(b,a,x) 其中3a=[a0,a1,…,aN], b=[bo,b1,…,bM] 分别表示差分方程系数。X 表示输入序列，y 表示输出序列。输 出序列的长度和序列相同。 例1-40 y[k]=1 M试用 M=9 点滑动平均系统M ?1 n=0∑ x[k ? n]处理例 1-39 中的受噪声干扰的正弦信号。 解： 由式（1-151）可知，M 点滑动平均系统可看成是 N=0 的 差分方程。在调用 filter 函数时，调用参数 a=1,b 为有 M 个元素的向 量，b 中每个元素的值均为 1/M。 %program 1_4 M=9; N=100;k=0:N; s=10*sin(0.02*pi*k); x=s+randn(1,N+1); b=ones(M,1)/M; a=1; y=filter(b,a,x); plot(k,y,s,’:’); xlabel(‘k’); ylabel（ ‘幅度’ ） legend(‘y[k]’,’s[k]’); 程序的运行结果如图 1-52 所示。图中的虚线表示未受噪声干扰的原 信号 s[k],实线为 9 点滑动的响应 y[k]。比较图 1-50 的信号可以看出，4滑动平均去噪系统滤出了受干扰信号中的大部分噪声，输出信号相对输入信号有 4 个样本的延迟。 1.9.4 利用 MATLAB 计算 DTET当序列的 DTET 可写成 e ? j? 的有理多项式时，可用 MATLAB 信号处 理工具箱提供的 freqz 函数计算 DTFT 的抽样值。另外，可用 MATLAB 提供的 abs、angle、real、imag 等基本函数计算 DTFT 的 幅度、相位、实部、虚部。若 X（ e j? ）可表示为B(e j? ) bo + b1e ? j? + ... + bMe? j?M = X（ e ）= A(e j? ) a 0 + a1e? j? + ... + aNe? j?Nj?则 freqz 的调用形式为 X=freqz(b,a,w) (1-153)式（1-153）中的 b 和 a 分别是表示式（1-152）中分子多项式和分 母多项式系数的向量，即 b=[b0,b1,…,bM] a=[a0,a1,…,aN] w 为抽样的频率点，在以式（1-153）形式调用 freqz 函数时，向量 w 的长度至少为 2。返回值 X 就是 DTFT 在抽样点 w 上的值。注意 一般情况下，函数 freqz 的返回值 X 是复数。 例1-41 H（z）= 已知离散系统的 H（z）为(0.5009 ? 1.0019 z ?1 + 0.5009 z ?2 )(0.5320 + 1.0640 z ?1 + 0.5320 z ?2 ) (1 ? 0.8519 z ?1 + 0.4167 z ?2 )(1 + 0.8519 z ?1 + 0.4167 z ?2试画出该系统的幅度响应。 解：5%program 1_5计算系统幅度响应b1=[0.9 0.5009]; b2=[0.0 0.5320]; a1=[1.9 0.4167]; a2=[1.9 0.4167]; b=conv(b1,b2);%计算 H（z）的分子多项式 a=conv(a1,a2);%计算 H（z）的分母多项式 w=linspace(0,pi,512); H=freqz(b,a,w); plot(w/pi,abs(H)); ylabel(‘幅度’); xlabel(‘Normalized frequency’); 系统幅度响应如图 1-53 1.9.5 部分分式法的 MATLAB 实现MATLAB 的信号处理工具箱提供了一个计算 X（z）的部分分 式展开的函数，它的调用形式如下 [r,p,k]=resifuez(b,a) 其中调用参数 b,a 分别表示用 z ?1 表示 X（z）的分子和分母多项式。 如果 X（z）的部分分式展开为 X(z)=r1 r2 r3 r4 + + + + k1 + k 2 z ?1 ?1 ?1 ?1 1 ? p1z 1 ? p2 z 1 ? p3 z (1 ? p3 z ?1 ) 2则 residuez 的返回参数 r,p,k 分别为 r=[r 1 r 2 r 3 r 4 ]6p=[p1 p2 p3 p3] k=[k1 k2] 同一极点 p3 在向量 p 中出现了两次，表示 p3 是一个二阶的重极点。 Residuez 也用于由 r,p,k 计算 z ?1 表示的 X（z）的分子和分母多项式， 其调用形式为 [b,a]=residuez(r,p,k) 例 1-43 X（z）= （ ） 试用 MATLAB 计算1.5 + 0.98 z ?1 ? 2.608 z ?2 z + 1.2 z ?3 ? 0.144 z ?4 z & 0 .6 1 ? 1.4 z ?1 + 0.6 z ? 2 ? 0.072 z ?3的部分分式展开。 解：计算部分分式展开的 MATLAB 程序如下 %program 1_6 部分分式展开 b=[1.5,0.98,-2.608,1.2,-0.144]; a=[1,-1.4,0.6,-0.072]; [r,p,k]=residuez(b,a); disp(‘留书’)；disp(r’) disp(‘极点’)；disp(p’) disp(‘常数’)；disp(k’); 程序运行结果为 留数 0.0i 0.0i 0.3000 极点 0.0i 0.0i 0.20007常数 0 2 部分分式展开的结果为 X(z)=0 .7 0 .5 0 .3 + + + z ?1 ?1 ?1 2 ?1 1 ? 0 .6 z (1 ? 0.6 z ) 1 ? 0 .2 zz 反变换的结果为 x[k]=(0.7x0.6 k +0.5x(k+1)0.6 k +0.3x0.2 k )u[k]+2 δ [k-1] 1.9.6 用 MATLAB 计算系统的零极点 MATLAB 信号处理工具箱提供的 tf2zp 和 zp2tf 函数可进行系统 函数的不同表示形式的转换。z 正幂有理多项式表示的系统函数为H ( z) = b(1) z M + b(2) z M ?1 + ... + b( M ) a (1) z N + a (2) z N ?1 + ... + a ( N )(1-154)用零点、极点和常数表示的一阶因子形式的系统函数为H ( z) = k ( z ? z (1))( z ? z (2))...( z ? z ( M )) ( z ? p (1))( z ? p (2))...( z ? p ( N ))（1-155）MATLAB 函数[z,p,k]=tf2zp(b,a)将有理多项式表示的系统函数转换为 一阶因子形式的系统函数。 [b,a]=zp2tf(z,p,k)将一阶因子形式的系统函 数转换为有理多项式表示的系统函数。 例 1-44 试将下面的系统函数表示为一阶因子形式。H ( z) = 1 + 0.04 z ?2 1 ? 0.8 z ?1 + 0.16 z ?2 ? 0.128 z ?3（1-156）解：用 z 正幂有理多项式表示的系统函数为H ( z) = z 3 + 0.04 z z 3 ? 0.8 z 2 + 0.16 z ? 0.128%program 1_7 确定一阶因子形式的系统函数8b=[1 0 0.04 0]; a=[1 -0.8 0.16 -0.128]; [z,p,k]=tf2zp(b,a); disp(‘零点’);disp(z’); disp(‘极点‘)；disp(p’); disp(‘常数’)；disp(k’); 程序运行结果为 零点 0 极点 0.0+0.4000i 常数 1 系统函数的一阶因子形式为H ( z) = z ( z ? 0.2 j )( z + 0.2 j ) ( z ? 0.8)( z ? 0.4 j )( z + 0.4 j )0+0.2000i0-0.2000i0.0i（1-157）为了在 H（z）的表达式中不出现复数，也可将式（1-157）等价的写 成二阶因子的形式1 + 0.04 z ?2 H ( z) (1 ? 0.8 z ?1 )(1 + 0.16 z ?2 )（1-158）当 H（z）是表示为式（1-154）的形式时，可用[z,p,k]=tf2zp(b,a) 求出系统的零极点，从而可根据极点的位置判断系统的稳定性，也可 用函数 zplane(b,a)画出 z 平面的零极点分布来判断系统的稳定性。9例１-45已知一因果的ＬＴＩ系统的函数为z 2 + 2z + 1 z 3 ? 0.5 z 2 ? 0.005 z + 0.3H ( z) =试用函数 zplane 画出系统的零极点分布，并判断系统的稳定性。 解： %program 1_8 系统零极点分布 b=[1 2 1]; a=[1 -0.5 -0.005 0.3]; zplane(b,a); 图１－５４画出了系统零极点分布。图中符号 Ο 表示零点。符号 Ο 旁 的数字表示零点的阶数， 符号 x 表示极点， 图中的虚线画的是单位圆。 由图可知，该系统的极点全在单位圆内，故系统是稳定的。 １.９.７ ９７ 例 1-46 简单数字滤波器的设计 已知信号 x[k]中含有角频率分别为 ?1 = 0.015π 和 ? 2 = 0.4π的离散正弦信号。试设计一高通滤波器，滤出信号 x[k]中低频分量， 保留信号 x[k]高频分量。 解：为了简单起见，假设数字滤波器是一个具有如下形式的长度 为３的ＦＩＲ系统 h[0]=h[2]= α ,h[1]= β 系统的查分方程y[k ] = αx[k ] + βx[k ? 1] + αx[k ? 2]系统的频率响应为10H (e j? ) = h[0] + h[1]e ? j? + h[2]e ? j 2?= α (1 + e ? j 2? ) + β e ? j? = (2α cos ? + β )e ? j?由上式可知系统的群延迟为τ g (? ) = 1为了滤除低频分量，系统需满足?2 cos(?1 ) ?2 cos(? ) 2 ?1? 1? ??α ? ?0? ? β ? ＝ ?1 ? ? ? ? ?可用ＭＡＴＬＡＢ命令Ａ／Ｂ求解上面的的方程组得：α = ?0.7248β = 1.4479下面的 MATLAB 程序实现了滤波器的设计及信号的滤波。 ％progam 1_9 简单数字滤波器的设计 ％确定滤波器系数 w1=0.015*W2=0.4*N=50; A=[2*cos(w1) 1;2*cos(w1) 1];B=[0;1] x=A/B; a=1;b=[x(1) x(2) x(1)]; %产生两个余弦信号 k=0:N; x1=cos(w1*k); x2=cos(w2*k); ％信号滤波 y=filter(b,a,x1+x2);11plot(k,y,`r`,k,x2,`b--`,x1+x2,`k:`); ylabel(`幅度`);xlabel(`k`) legend(`y[k]`,`x2[k]`,`x1[k]+x2[k]`); 图１－５５画出了程序运行结果。由图可以看出在 k=0，１时， 输出响应有波动。这是因为系统响应中存在瞬态响应。当 k ≥ 2 时，系 统进入稳态响应。滤波后的信号相对与原信号有一个样本的延迟，着 是因为在 ? = ? 2 时，τ p (? 2 ) = τ g (? 2 ) = 1 １.９.８ 语音信号的滤波 ９８ 例 1-47 已知一段语音信号中混入了一频率 f n =1200Hz 正弦型干扰信号。语音信号的抽样频率 f sam =22050Hz。试用二阶带阻滤波器 滤除语音信号中的噪音信号。 解：干扰信号的数字频率 ? 0 为? 0 = 2πf n / f sam = 0.3419由式（1-110）可得β = cos ? 0 = 0.9421取 ?? 3dB = 0.06 ，由式（1-111）可得α 2 ? 2α /(?? 3dB ) + 1 = 0解上述方程得 α 的值为 1.0619 和 0.9417。当 α 的值为 1.0619 时，系 统有一个极点在单位圆外，当 α 的值为 0.941７时，系统的二个极点 在单位圆内。为了保证系统的稳定，取 α ＝0.9417。由 α 和 β 的值可 得系统函数为H BS ( z ) = 0.9709 ? 1.8293 z ?1 + 0.9709 z ?2 1 ? 1.8293 z ?1 + 09417αz ?212%program 1_10 语音信号滤波 Fn=1200;w 3dB =0.06; %读入语音信号 ［x,Fs］=wavread(`sample`); %播放语音信号 sound(x,Fs); %设计滤波器 w0=2*pi*Fn/Fs; beta=cos(w0); alpha=min(roots([1-2/cos(w 3dB ) 1])); a=[1,-beta*(1+alpha),alpha]; b=[1 -2 *beta 1]*(1+alpha)/2 ; %信号滤波 y=filter(b,a,x); %播放滤波后语音信号 sound(y,Fs); 图１－５６画出了滤波前后语音信号的频谱。 滤波前的语音信号 的频谱在１２００Ｈz 有一高峰，这是干扰所造成的。由滤波后语音 信号的频谱可看出，滤波器成功地制造了语音信号中的干扰。132.6利用 MATLAB 实现信号 DFT 的计算2.6.1 利用 MATLAB 计算信号 DFT 在 MATLAB 信号处理工具箱中， 函数 dftmtx(N)可用来产生 N × N 的 DFT 矩正 D N 。N × N 的 IDFT 矩正 D N ?1 可用函数 conj(dfmtx(N))/N 来确定。 此外， MATLAB 提供了 4 个内部函数用于计算 DFT 和 IDFT， 它们分别是： fft(x), fft(x,N), ifft(X), ifft(X,N)fft(x) 计算 M 点的 DFT。M 是序列 x 的长度，即 M=length(x)。 fft(x,N) 计算 N 点的 DFT。若 M&N，则将原序列截短为 N 点序 列，再计算其 N 点 DFT；若 M&N，则将原序列补零至 N 点，然后计 算其 N 点 DFT。 ifft(X) 计算 M 点的 IDFT。M 是序列 X 的长度。 ifft(X,N) 计算 N 点 IDFT。若 M&N，则将原序列截短为 N 点序 列，再计算其 N 点 IDFT；若 M&N，则将原序列补零至 N 点，然后 计算其 N 点 IDFT。 为了提高 fft 与 ifft 的运算效率， 应尽量使序列长度 M=2 r ，或对 序列补零使 N=2 r 。实现例 2-7 的程序如下： %program 2_1 计算有限长余弦信号频谱 N=30；%数据长度 L=512；%DFT 的点数 f1=100; f2=120; fs=600; T=1/14ws=2*pi* t=(0:N-1)*T; f=cos(2*pi*f1*t)+cos(2*pi*f2*t); F=fftshift(fft(f,L)); w=(-ws/2+(0:L-1)*ws/L)/(2*pi); plot(w,abs(F)); ylabel(`幅度谱`) 实现例 2-8 的程序为： %program 2_2 利用 Hamming 计算有限长余弦信号频谱 N=50；%数据长度 L=512；%DFT 的点数 f1=100; f2=120; fs=600; %N=25； T=1/ ws=2*pi* t=(0:N-1)*T; f=cos(2*pi*f1*t)+cos(2*pi*f2*t); wh=(hamming(N))`; f=f.* F=fftshift(fft(f,L)); w=(-ws/2+(0:L-1)*ws/L)/(2*pi); plot(w,abs(F));15ylabel(`幅度谱`) 例 2-11 已知一长度为 16 点的有限序列x[k ] = cos( 2π rk ), N = 16, r = 4 N试用 MATLAB 计算序列 x[k]的 16 点和 512 点 DFT。 解： %progam 2_3 Numerical Computation of DTFT Using FFT k=0:15; f=cos(2*pi*k*4./16); F_16=fft(f); F_512=fft(f,512); L=0:511; plot(L/512,abs(F_512)); plot(k/16,abs(F_16),`0`); set(gca,`xtick`,[0,0.25,0.5,0.75,1]); set(gca,`ytick`,[0,2,4,6,8]); xlabel(`Normalized frequency`); ylabel(`Magnirude`); hold off 从图 2-19 可见，对长度为 16 的序列 x[k]进行 512 点 DFT，可以 获得频谱函数 X（ e j? ）更多的细节，因为序列 N 点的离散傅立叶变16换 X N [m] 就是序列 X（ e j? ）一个周期的 N 个等间隔抽样，尽管从信 息表达的角度，由序列 x[k]16 点 DFT X[m]足以恢复原序列 x[k]。 2.6.2 利用 MATLAB 实现由 DFT 计算线性卷积 并与 例 2-12 试利用 MATLAB 实现由 DFT 计算有限序列线性卷积， 线性卷积直接计算的结果进行比较。 在 序列线性卷积的直接结果是通过函数 conv(x,h) 解： MALAB 中， 来实现的。 %program 2_4 linear Convolution through DFT x=[1 2 0 1]; h=[2 2 1 1]; %Determine the length of the result of convolution L=length(x)+length(h)-1; %Compute the DFT by zero-padding XE=fft(x,L); HE=fft(h,L); %Determine the IDFT of the product y1=ifft(XE.*HE); %plot the sequence generated by circular convolution %and the error from direct linear convolution k=0:L-1; subplot(2,1,1); setm(k,real(y1));axis([0 6 0 7]);17title(`Result of Linear Convolution`); xlabel(`Time index k`);ylabel(`Amplitude`); y2=conv(x,h); error=y1-y2; subplot(2,1,2); stem(k,abs(error)); xlabel(`Time index k`);ylabel(`Amplitude`); title(`Error Magnitude`); 由图 2-20（b）可见，由 DFT 计算的线性卷积的误差非常小，这 些误差主要是计算机在计算 DFT 时，由计算机有限字长而产生的舍 入误差与计算误差。 此外，MATLAB 信号处理工具箱中提供了 fftfilt 函数，该函数是 基于重叠相加法的原理，可以实现长序列与短序列之间的线性卷积， 其调用形式为 y=fftfilt(h,x,n) 其中 h:短序列； x：长序列； n：DFT 的点数。一般选择 n 为 2 正整次幂，以提高 DFT 运算的 效率。183.6 线性调频 z 变换算法 例 3-3 已知序列 x[k]= (0.8) k ,( 0 ≤ k ≤ 12 ),试计算序列 x[k]在单位 圆上的 CZT，其中 θ 0 = 解： 由于是计算序列 x[k]在单位圆上的 CZT，故 A0=1，W0=1。计算 序列 czt 的 MATLAB 程序如下： %program 3_1 % Compute CZT of sequence x[k] N=13; n=0:N-1; xn=(0.8).^n; sita=pi/4; fai=pi/6; A=exp(j*sita); W=exp(-j*fai); M=8; y=czt(xn,M,W,A); subplot(2,1,1) stem(n,xn); xlabel(`k`);title(`sequence x[k`]); m=0”M-1; subplot(2,1,2);19π4,? 0 =π6,M = 8。%length of x[k]%phase of start point %phase interval %complex starting point %complex ratio between points %length of czt %call czt functionstem(m,abs(y)); xlabel(`m`);title(`CZT of x[k]`); 运行结果如图 3-22 所示。 4.5 用 MATLAB 实现滤波器设计4.5.1 用 MATLAB 实现模拟低通滤波器的设计 MATLAB 信号处理工具箱提供了常用的设计模拟低通滤波器的 MATLAB 函数。在实现应用中，可方便地调用这些函数完成模拟滤 波器的设计。关于这些函数现分别介绍如下： Butterworth 滤波器 [N，wc]=buttord(wp,ws,Ap,As,`s`) [num,den]=butter(N,wc,`s`) 根据 BW 型滤波器的设计指标， 利用 MATLAB 函数 buttord 获得 BW 型滤波器参数 N 和 wc。函数 buttord 的输入参数 wp 和 ws(rad/s) 分别表示滤波器的通带和阻带截频，Ap 和 As(dB)表示滤波器的通带 和阻带衰减。`s`表示所设计的是模拟滤波器。函数 buttord 返回参数 N 为 BW 滤波器的阶数，wc(rad/s)等于 BW 滤波器 3 dB 截频 ω c 。由 于 wc 由阻带方程确定，故由参数 N、wc 得出的滤波器在阻带刚好满 足设计指标，在通带将超过设计指标。 当 BW 型滤波器的参数 N、 确定后， wc 可用 MATLAB 函数 butter 获得 BW 型滤波器系统函数 H LP (s) 的分子多项式（num）和分母多项 式（den） 。 Chebyshev I 型滤波器20[N，wc]=cheblord(wp,ws,Ap,As,`s`) [num,den]=cheby1(N,Ap,wc,`s`) MATLAB 函数 cheblord 返回参数 N 表示 CB I 型滤波器的阶数， wc(rad/s)等于 wp。参数 N、wc 的取值可使 cheby1 设计出的 CB I 型 、 滤波器在通带刚好满足设计指标。 cheby1 函数利用参数 N、wc 和 Ap 确定 CB I 型滤波器系统函数H LP (s ) 的分子多项式（num）和分母多项式（den） 。Chebyshev II 型滤波器 [N，wc]=cheb2ord(wp,ws,Ap,As,`s`) [num,den]=cheby2(N,As,wc,`s`) MATLAB 函数 cheb2ord 返回参数 N 表示 CBII 型滤波器的阶数， wc(rad/s)的取值可使 cheby2 设计出的滤波器在通带刚好满足设计指 标。 cheby2 函数利用参数 N、wc 和 As 确定 CBII 滤波器系统函数H LP (s ) 的分子多项式（num）和分母多项式（den） 。椭圆滤波器 [N，wc]=ellipord(wp,ws,Ap,As,`s`) [num,den]=ellip(N,Ap,wc,`s`) MATLAB 函数 ellipord 返回参数 N 表示椭圆滤波器的阶数， wc=wp。 MATLAB 函数 ellip 确定阶数为 N，通带衰减为 Ap(dB),阻带衰 减为 As(dB),通带截频为 wc 的椭圆滤波器系统函数 H LP (s) 的分子多项21式和分母多项式。在滤波器的阶数确定后，采用了 4.1.3 节中方案 A 重新计算滤波器的参数 k1。故设计出的椭圆滤波器阻带衰减超过设 计指标，而其他指标刚好满足设计要求。 例 4-15 设计一个满足下列指标的模拟 BW 低频滤波器。f p = 1kHz, f s = 3kHz , A p ≤ 1dB, As ≥ 50dB解： %program 4_1 设计模拟 Butterworth 低频滤波器 %滤波器指标 wp=2*pi*1000;ws=2*pi*3000;Ap=1;As=50; %设计滤波器 [N，wc]=buttord(wp,ws,Ap,As,`s`) [num,den]=butter(N,wc,`s`) fprintf(`Order of the fiter=%.Of\n`,N) disp(`Numerator polynomial`); fprintf(`%.4e\n`,num); disp(`Denominator polynomial`); fprintf(`%.4e\n`,den); %计算 Ap,As 及增益响应 omega1=linspace(0,wp,500); omega2=linspace(wp,ws,200); omega3=linspace(ws,5*1000*pi*2,500); h1=20*log10(abs(freqs(num,den,omega1)));22h2=20*log10(abs(freqs(num,den,omega2))); h3=20*log10(abs(freqs(num,den,omega3))); fprintf(`Ap=%.4f\n`,max(-h1)); fprintf(`As=%.4f\n`,min(-h3)); plot([omega1 omega2 omega3]/(2*pi),[h1 h2 h3]);xlabel(`Frequency in Hz`);ylabel(`Gain in dB`); 程序运行后的结果如下： order of the filter=6 Numerator polynomial 0....... Denominator polynomial 1.000e+000 2...... Ap=0.7488 As=50.0000 图 4-25 画出了该滤波器的增益响应。 例 4-16 利用模拟椭圆低通滤波器重新设计例 4-15。 解：只需将例 4-15 程序的第 6 行和第 7 行改为 [N，wc]=ellipord(wp,ws,Ap,As,`s`) [num,den]=ellip(N,Ap,wc,`s`) 即可。程序的运行结果如下： Order of the filter=423Numerator polynomial 3..... Denominator polynomial 1.. Ap=1.0000 As=51.9e+003 5..图 4-26 画出了该滤波器的增益响应。比较例 4-15 和例 4-16 可知，实 现同样设计指标的滤波器，椭圆滤波器所需阶数较低。 4.5.2 用 MATLAB 实现模拟域频率变换MATLAB 信号处理工具箱提供了实现四种模拟域频率变换的函 数，它们分为： 低通到低通的变换 [numt,dent]=1p21p(num,den,w0) 低通到高通的变换 [numt,dent]=1p2hp(num,den,w0) 低通到带通的变换 [numt,dent]=1p2bp(num,den,w0,B) 低通到带阻的变换 [numt,dent]=1p2bs(num,den,w0,B) 其中 num、den 分别表示变换前模拟滤波器系统函数的分子多项式和 分母多项式。numt、dent 分别表示变换后模拟滤波器系统函数的分子 多项式和分母多项式。w0 和 B 为变换中的参数。24例 4-47试设计一个满足下列指标的模拟 BW 型带阻滤波器A p ≤ 1dB, As ≥ 10dB, ω p1 = 6rad / s, ω p 2 = 13rad / s, ω s1 = 9rad / s, ω s 2 = 11rad / s解： （1）由式（4-77）确定低通到带阻变换中的参数 B=2， （2）由式（4-79）? p1 = 0.1905,? p 2 = ?0.3714 ω 0 = 9.9499由式（4-80）得原型低通滤波器的通带截频 wp_bar 为 wp_bar=0.3714 （3）由 [N，wc]=buttord(wp_bar,1,Ap,As,`s`); [num,den]=butter(N,wc,`s`); 得原型低通滤波器为H L ( s) = 0.3333 s + 0.8165s + 0.33332（4）由 [numt,dent]=1p2bp(num,den,w0,B)； 得带阻低通滤波器为H BS ( s ) =s 4 + 1.98 × 10 2 s 2 + 9.801 × 10 3 s 4 + 4.899 s 3 + 2.1 × 10 2 s 2 + 4.85 × 10 2 s + 9.801 × 10 3图 4-27 画出了该滤波器的增益响应曲线。该滤波器的通带和阻带衰 假分别为 Ap1=0.05dB, Ap2=0.68dB, As1=As2=10dB。 4.5.3 脉冲响应不变的 MATLAB 实现25MATLAB 信号处理提供的 impinvar(num,den,Fs)函数， 可实现用 脉冲响应不变法将模拟滤波器转化为数字滤波器，起调用形式为： [numd,dend]=impinvar(num,den,Fs) 式中 num 和 den 分别表示模拟滤波器系数函数 H（s）的分子多项式 和分母多项式，Fs 是脉冲响应不变法中的抽样频率，单位是 Hz。输 出变量 numd 和 dend 分别表示数字滤波器的系统函数 H(z)分子多项 式和分母多项式。 若某个因果模拟滤波器的系统函数为H (s) = 2 s + 14 s + 2s + 52取 T=1，则由 [numd,dend]=impinvar([2 14],[1 2 5],1) 得 numd=2.0000 dend=1.3 0.3所以由脉冲响应不变法获得的数字滤波器为H ( z) = 2 + 2.3133 z ?1 1 + 0.3062 z ?1 + 0.1353例 4-18 利用 Butterworth 低通滤波器及脉冲响应不变法设计满 足下列指标的数字滤波器? p = 0.1πrad , ? s = 0.4πrad , A p ≤ 1dB, As ≥ 25dB解：设计满足上述指标数字滤波器的 MATLAB 程序如下： %program 4_2 脉冲响应不变法设计 BW 型 LP DF %DF BW LP 指标26Wp=0.1*Ws=0.4*Ap=1;As=25; Fs=1;%抽样频率（Hz） %确定模拟 BW 指标 Wp=Wp*Fs;Ws=Ws*Fs; %确定 AF 阶数 N=buttord(Wp,Ws,Ap,As,`s`); %由通带指标确定 3-dB 截频 Wc=Wp/(10^(0.1*Ap)-1)^(1/2/N); %确定 BW AF [numa,dena]=butter(N,Wc,`s`); %确定 DF[numd,dend]=impinvar(numa,dena,Fs); w=linspace(0,pi,512); h=freqz(numd,dend,w); %幅度归一化 DF 的幅度响应 norm=max(abs(h)); numd=numd/ plot(w/pi,20*log10(abs(h)/norm)); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); disp(`Numerator polynomial`); fprintf(`%.4e\n`,numd);27disp(`Denominator polynomial`); fprintf(`%.4e\n`,dend); %计算 Ap 和 As w=[Wp Ws];h=freqz(numd,dend,w); fprintf(`A p=%.4f\n`,-20*log10(abs(h(1)))); fprintf(`As=%.4f\n`,-20*log10(abs(h(2)))); 运行结果为 Numerator polynomial -5.... Denominator polynomial 1. -2.. -4. Ap=0.9985 As=30.3240 Ap=0.9985 As=30.324028图 4-28 画出了该滤波器的增益响应 4.5.4 双线性变换法的 MATLAB 实现MATLAB 信号处理工具箱提供的（num,den,Fs）函数可以用来实 现双线性变换，其调用形式为： [numd,dend]=bibinear(num,den,Fs) 式中 num 和 den 分别表示模拟滤波器系统函数 H（s）的分子多项式 和分母多项式，Fs=1/T. 输出变量 numd 和 dend 分别数字滤波器系统函数 H（s）的分子 多项式和分母多项式。 例 4-19 三阶归一化 Butterworth 滤波器的系统函数为H ( s) = 1 ( s + 1)( s 2 + s + 1)解：取 T=2，则由 den=conv([1 1],[1 1 1]); [numd,dend]=bilinear([1],den,0.5) 得经双线性变换后的数字滤波器的分子多项式和分母多项式为 numd=0.0 0.7 dend=1.0 0.0 注意 0.1667 是 1/6 的近似值，0.3333 是 1/3 的近似值，所以双线 性变换后的数字滤波器的系统函数 H（z）为H ( z) = (1 / 6)(1 + z ?1 ) 3 1 + (1 / 3) z ?2例 4-20 利用 CB I 型模拟滤波器和双线性变换法，设计一满足29下列指标的数字低通滤波器。? p = 0.1πrad , ? s = 0.4πrad , A p ≤ 1dB, As ≥ 10dB解：取 T=2，由 ω =? 2 tan( ) 得模拟滤波器的频率指标 T 2ω p = 0.1584, ω s = 0.7265由 [N,wc]=cheb1ord(0.5,1,10,`s`) [num,den]=cheby1(N,1,wc,`s`) 得模拟滤波器的分子多项式和分母多项式为 num= 0 0 0.7den=1.9 即H a (s) =0.0246 s + 0.1739 s + 0.02772再由 [numd,dend]=bilinear(num,den,0.5) 得双线性变换后的数字滤波器的分子多项式和分母多项式 numd=0.0 0.0205 dend=1.5 0.7106 即H ( z) 0.0205(1 + z ?1 ) 2 1 ? 1.6185 z ?1 + 0.7106 z ?2上式表示的数字滤波器的增益响应如图 4-29 所示。该滤波器的通带 和阻带衰减分别为 AP ≈ 1dB, As ≈ 27dB 。 4.5.5 用 MATLAB 实现数字滤波器设计30MATLAB 也提供了直接设计 IIR 数字滤波器的函数。设计的基 本思想是用式（4-103）将数字滤波器的频率指标转化为模拟滤波器 的指标，然后设计模拟滤波器，最后用双线性变换把模拟滤波器转换 为数字滤波器。由于在设计中用的模拟滤波器的类型不同，所以给出 了以下四组不同的函数 BW 型数字低通滤波器 [N,Wc]=buttord(Wp,Ws,Ap,As) [num,den]=butter(N,Wc) 在函数 buttord 函数中，调用参数 Wp、Ws 是数字低通滤波器的 归一化的通带和阻带截频。例如要求数字滤波器的? p = 0.1π , ? s = 0.4π ，则 Wp=0.1,Ws=0.4。Ap 和 As 为滤波器的通带和阻带的衰减（dB） 。返回的参数为Ｎ和 Wc，其中 N 为滤波器阶数。 函数 butter 中，调用参数 N 和 Wc 由函数 buttord 确定。返回参 数 num 和 den 分别是数字滤波器的分子多项式和分母多项式的系数。 CB I 型数字低通滤波器 [N,Wc]=cheb1ord(Wp,Ws,Ap,As) [num,den]=cheby1(N,Ap,Wc) CB II 型数字低通滤波器 [N,Wc]=cheb1ord(Wp,Ws,Ap,As) [num,den]=cheby1(N,Ap,Wc) 椭圆型数字低通滤波器 [N,Wc]=ellipord(Wp,Ws,Ap,As)31[num,den]=ellip(N,Ap,As,Wc) 利用上述四组 MATLAB 函数，也可设计数字高通、带通和带阻 滤波器。调用方法可参见 MATLAB 手册。 例 4-21 器。? p = 0.2πrad , ? s = 0.4πrad , A p ≤ 0.5dB, As ≥ 20dB设计满足下列指标的 BW 型和椭圆型数字低通滤波解： %program 4_3 设计数字低通滤波器 %DF 指标 Wp=0.2;Ws=0.4; Ap=0.5;As=20; %设计 BW 型 DF LP [N,Wc]=buttord(Wp,Ws,Ap,As); [b1,a1]=butter(N,Wc); %设计 C 型 DF LP [N,Wc]=ellipord(Wp,Ws,Ap,As); [b2,a2]=ellip(N,Ap,As,Wc); disp(`BW 型分子多项式`)； fprintf(`%.4e\n`,b1); disp(`BW 型分母多项式`)； fprintf(`%.4e\n`,a1); disp(`C 型分子多项式`)；32fprintf(`%.4e\n`,b2); disp(`C 型分母多项式`)； fprintf(`%.4e\n`,a2); w=linspace(0,0.6*pi,200); h1=freqz(b1,a1,w); h2=freqz(b2,a2,w); plot(w/pi,20*log10(abs(h1)),w/pi,20*log10(abs(h2))); legend(`BW 型`，`C 型`)； xlabel(`Normalized frequency`); ylabel(`Gain,dB`); axis([0 0.6 -50 1]); 程序运行的结果为 BW 型分子多项式 4....2.. BW 型分母多项式 1.. -2. -5.. -1.C 型分子多项式 9. -1. -1.. C 型分母多项式 1. -1.. -4.33图 4-30 画出了所设计的滤波器的增益响应。由运行结果知对同样 的设计指标，用 BW 型模拟滤波器设计出的数字滤波器是五阶的，而用 C 型模拟滤波器设计出的数字滤波器只有三阶。 5.5 利用 MATLAB 实现 FIR 滤波器设计5.5.1 窗函数法的 MATLAB 实现 窗函数的计算 窗函数的计算 MATLAB 提供了许多常用的窗函数， 其中部分窗函数的调用形 式为 w=hanning(N) w=hamming(N) w=blackman(N) w=Kaiser(N,beta) 其中 N 是窗函数的长度,beta 是控制 kaiser 窗形状的参数。返回的变 量 w 是一个长度为 N 的列向量，给出窗函数在 N 点的取值。 窗函数法设计 FIR 滤波器 窗函数法设计 FIR 滤波器一般分为 3 个步骤：第 1 步估计 FIR 滤波器阶数 M（或长度 N） 。如果用 Kaiser 窗时可用式（5-41）估计 FIR 滤波器阶数；第 2 步确定所用的窗函数并计算出窗函数的值；第 3 步计算理想滤波器的单位脉冲响应并用窗函数将其截断即得所设计 的 FIR 滤波器的 h[k]。 例 5-5 用 I 型线性相位滤波器设计一满足下列指标的 FIR 高通滤 波器。 ? p = 0.8πrad , ? s = 0.7πrad , Ap ≤ 0.3dB, As ≥ 40dB ，34解：表 5-2 可知用 Hann 窗可使滤波器满足指标。由过渡带宽宽 度得滤波器长度 N 需满足N≥ 6.2π = 62 ?p ? ?s由于要求用 I 型滤波器，取 N=63。取理想高通滤波器的截频为?c = (?p + ?s ) / 2 = 0.75πrad由式（5-28）得hd [k ] = δ [k ? 31] ? 0.75S a [0.75π (k ? 31)]用长度 N=63 的 Hamming 窗截断 hd [k ] 即得所要求的 FIR 高通滤波器 的 h[k]为h[k ] = hd [k ]w[k ]完成上述计算过程的 MATLAB 程序如下： %program 5_1 窗函数法设计 FIR 高通 Wp=0.8*Ws=0.7*Ap=0.3;As=40; %确定滤波器阶数 N=ceil(6.2*pi/(Wp-Ws)); N=mod(N+1,2)+N; M=N-1; fprintf(`滤波器阶数=%.of\n`,M); w=hanning(N)`; %理想低通截频 Wc=(Wp+Ws)/2; k=0:M;35hd=-(Wc/pi)*sinc(Wc*(k-0.5*M)/pi),hd(0.5*M+1)=hd(0.5*M+1)+ 1; h=hd.*w; omega=linspace(0,pi,512); mag=freqz(h,[1],omega); magdb=20*lg(abs(mag)); plot(omega/pi,magdb); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 图 5-24 画出了用 hann 窗设计的 FIR 高通滤波器的增益响应。 Hann 窗 设 计 的 FIR 高 通 滤 波 器 的 通 带 和 阻 带 衰 减 分 别 为 Ap=0.034dB、As=44dB。作为比较，图 5-25 画出了用 Kaiser 窗设计 的 FIR 高通滤波器的增益响应。用 Kaiser 窗设计出的滤波器长度 N=47，Kaiser 窗的参数 β =3.3953。滤波器通带和阻带衰减分别为 Ap=0.057dB,As=42dB。比较两种设计结果可知，所设计出的滤波器 都满足设计指标。用 Kaiser 窗得出的滤波器阶数较低，用 Hann 设计 出的滤波器阻带的波纹衰减较快。 用 Kaiser 窗设计 FIR 滤波器阶数的 MATLAB 函数基本调用形式 为： [M,Wc,beta,ftype]=Kaiserord(f,a,dev) h=fir1(M,Wc,‘ftype’,window)36函数 Kaiserord 完成估计滤波器阶数 M 及 Kaiser 窗的参数 beta。函数 fir1 根据所用窗函数及理想滤波器的截频计算出 FIR 滤波器的 h。 函数 Kaiserord 调用参数 f 表示需设计的 FIR 滤波器的频带。如 FIR 滤波器的 B 个频带分别为0 ≤ ? ≤ πf 1 , πf 2 ≤ ? ≤ πf 3πf 4 ≤ ? ≤ πf 5 ,..., πf 2( B ?1)? 2 ≤ ? ≤ πf 2( B ?1) ?1 , πf 2 B ?2 ≤ ? ≤ π则 f 是一个有 2B-2 个元素的向量（B ≥ 2） ，其值为f = [ f1 , f 2 , f 3 , f 4 , f 5 ,... f 2 B ? 2 ]函数 kaiserord 调用参数 a 是一个 B 个元素的向量， 分别表示 FIR 滤波器在 B 个频带中的幅度值。 一般对通带取值为 1， 阻带取值为 0。 函数 kaiserord 调用参数 dev 是一个 B 个元素的向量，分别表示 FIR 滤波器在 B 个频带中的波动值。由前面的分析可知，Kaiser 窗设 计的 FIR 滤波器在各带的波动是相等的，所以 kaiserord 是利用最小 的波动值来估计滤波器的阶数。 函数 kaiserord 返回参数 M 及 beta,分别表示 FIR 滤波器阶数 M 及 kaiser 窗的参数 beta。返回参数 Wc 和 ftype 是函数 firl 的调用参 数。Wc 是 B-1 个元素的向量，如 Wc=[W1,W2,…,W B ?1 ] 则表示理想 FIR 滤波器的 B 个频带分别为0 ≤ ? ≤ πw1 , πw1 ≤ ? ≤ πw2 , πw2 ≤ ? ≤ πw3 ,..., πwB ?1 ≤ ? ≤ π返回参数 ftype 是一个字符串。 ftype 为空时， 表示滤波器为低通， ftype= ‘high’滤波器为高通，ftype=‘stop’滤波器为带阻，ftype=‘DC-0’ 多带滤波器第一个频带为阻带，ftype=‘DC-1’多带滤波器第一个频37带为通带。例如 f=[0.2,0.6], a=[1,0] 时，函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc=0.4, ftype=`` 当 f=[0.2,0.6], a=[0,1] 时，函数 kaiserord 返回值 Wc 和 ftype=`high` 当 f=[0.2,0.4,0.6,0.8], a=[1,0,1] 时，函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc=0.4, 当 f=[0.2,0.4,0.6,0.8], a=[0,1,0] 时，函数 kaiserord 返回值 Wc 和 ftype 分别为 Wc=[0.3,0.7], ftype=`DC-0` 函数 firl 的调用参数 M 表示滤波器的阶数，Wc 表示理想 FIR 滤 波器的 B 个频带，ftype 表示滤波器的类型，缺省值为空白。window 是一个长度为 M+1 的向量。如果调用是没给窗函数，程序 firl 自动 使用 Hamming 窗。 例 5-6 试用 Kaiser 窗设计满足下列指标的 FIR 高通滤波器? p = 0.4πrad , ? p = 0.6πrad , δ s = 0.01ftype=`stop`解：38%program 5_2 Kaiser 窗设计 FIR 高通滤波器 Rs=0.01; f=[0.4,0.6]; a=[0,1]; dev=Rs*ones(1,length(a)); [M,Wc,beta,ftype]=kaiserord(f,a,dev); %使滤波器为 I 型 M=mod(M,2)+M; h=firl(M,Wc,ftype,Kaiser(M+1,beta)); omega=linspace(0,pi,512); mag=freqz(h,[1],omega); plot(omega/pi,20*log10(abs(mag))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 图 5-26 画出了所设计的 FIR 高通滤波器的增益响应。 例 5-7 试用 Kaiser 窗设计满足下列指标的具有两个通带 FIR 滤 波器? s1 = 0.1πrad , ? p1 = 0.2πrad , ? p 2 = 0.4πrad , ? s 2 = 0.5πrad , ? s 3 = 0.6πrad , ? p 3 = 0.7πrad , ? p 4 = 0.8πrad , ? s 4 = 0.9πrad , δ s = 0.01解：只需将 program 5_2 第 4 和第 5 改为 f=[0.1 0.2 0.4 0.5 0.6 0.7 0.8 0.9]; a=[0,1,0,1,0];39并去掉其第 9 行即可。图 5-27（a）画出了取 Rs=0.01 时所设计的 FIR 滤波器的幅度响应。 由图可见所设计出的滤波器仅有一个阻带满足设 计指标， 其他两个阻带不满足设计指标。 5-27 b） 图 （ ） 画出了取 Rs=0.008 时 FIR 滤波器的幅度响应，这时滤波器的三个阻带均满足设计指标。 MATLAB 信号处理工具箱提供的另一个用窗函数法设计 FIR 滤 波器的函数为 h=fir2(M,f,a,npt,window) fir2 实现的基本思想是利用线性内插得出理想滤波器的频率响应H (e j? ) 。 然后对 H (e j? ) 均匀抽样得出 H (e j? ) 在 npt 个点 （缺点值为 512）上的抽样值。 H (e j? ) 的 npt 点的 IDFT 得 h d [k],最后用窗函数将 h d [k] 对 截断为长度为 M+1 的序列。 函数 fir2 的调用参数 N 表示的阶数，npt 表示 IDFT 的点数，缺 省值为 512。window 为所用窗函数序列，缺省时使用 Hamming 窗。 如果所设计的滤波器的 B 频率点上的幅度值为H (e j 2πfm ) = a m ,0 = f 1 ≤ f 2 ,..., ≤ f B = 1则 f 和 a 的值分别为 f=[f1,f2,..., f B ], a=[a1,a2,…, a B ] 函数 fir2 的返回参数 h 为长度为 M+1 的 FIR 滤波器系数。 例 5-8 试设计一个线性相位 FIR 滤波器，使其幅度响应逼近?1 ? 5 ? / π ? 0.6 H (e j? ) = ? ?00.4π ≤ ? ≤ 0.8π其他解： H (e j? ) 在 0.4π ≤ ? ≤ 0.6π 和 0.6π ≤ ? ≤ 0.8π 范围内都是 ? 的线40性函数，所以可用 fir2 进行设计。 %program 5_3 %通带响应为线性函数的 FIR 滤波器设计 M=input(`M=`); f=[0 0.4 0.6 0.8 1]; a=[0 0 1 0 0]; h=fir2(M,f,a,ones(1,M+1)); omega=linspace(0,pi,512); mag=freqz(h,[1],omega); h d =plot(omega/pi,abs(mag)); xlabel(`Normalized frequency`); ylabel(`幅度`)； 图 5-28 画出了 M=40 时所设计的 FIR 滤波器的幅度响应。 5.5.2 频率取样法的 MATLAB 实现用 ATLAB 实现频率取样法过程非常简单。例如对 I 型线性相位 滤波器，首先可以根据需设计的滤波器的频率响应由式（5-49）确定 取样点上的值 H d [m] ，然后对 H d [m] 做 M+1 点 IDFT 即可得到 h[k]。 例 5-9 用 频 率 取 样 法 设 计 一 个 M=63 （ II 型 ） ，?p = 0.5πrad , ?s = 0.6πrad 线性相位低通滤波器。解：由前面分析可知，如果对截频 ?c = 0.5π 的理想低通滤波器 进行取样，则设计出的滤波器在阻带的衰减不会超过 20dB。若采用 过渡带[ ?p, ?s ]之间加入一个幅度在 0 和 1 之间的过渡点来设计滤波41器，则过渡点的幅度值 T1 的大小将影响设计出的滤波器性能。能使 滤波器阻带衰减达到最大的过渡点幅度值 T1 称为最优值，最优值可 通过优化算法或实验得出。参考文献[2]列出了用优化算法得出最优 值。由于在通带内有 17 个取样点，利用参考文献[2]给出的表格得出 近似最优的 T1=0.38。确定了过渡点的幅度值，就可根据式（5-50） 进行编程。 %program 5-4 频率取样法设计 II 型线性相位低通 FIR M=63;Wp=0.5* m=0:(M+1)/2; Wm=2*pi*m./(M+1); mtr=floor(Wp*(M+1)/(2*pi))+2; Ad=[Wm&=Wp]; Ad(mtr)=0.38; Hd=Ad.*exp(-j*0.5*M*Wm); Hd=[Hd conj(fliplr(Hd(2:(M+1)/2)))]; h=real(ifft(Hd)); w=linspace(0,pi,1000); H=freqz(h,[1],w); plot(w/pi,20*log10(abs(H))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 图 5-29 画出了所设计的 FIR 滤波器的频率特性。由图可见，增加一42个过渡点后使滤波器的阻带衰减达到了 43dB。如要进一步增加阻带 衰减可考虑增加两个过渡点。查参考文献[2]中表格可知，最优值过 渡点的幅度值分别为 T2=0.59 和 T1=0.11。 5-30 画出了用两个过渡 图 点设计滤波器的频率特性， 由图可见滤波器的阻带衰减达到了 62dB。 例 5-10 利用频率取样法设计一个满足下列指标的 I 型线性相位高 通滤波器。?s = 0.5πrad , ?p = 0.6πrad解：为使在过渡带中有一个样本点，本例取 M=32。由于参考文 献[2]中表格给出的最优值是按低通滤波器得出的，因此并不适合高 通滤波器。由实验确定本例的近似最优值 T1=0.28。根据（5-49）编 制的 MATLAB 程序如下： %program 5_5 频率取样法设计 I 型线性相位高通 FIR M=32;Wp=0.6* m=0:M/2; Wm=2*pi*m./(M+1); mtr=ceil(Wp*(M+1)/(2*pi)); Ad=[Wm&=Wp];Ad(mtr)=0.28; Hd=Ad.*exp(-j*0.5*M*Wm); Hd=[Hd conj(fliplr(Hd(2:M/2+1)))]; h=real(ifft(Hd)); w=linspace(0,pi,1000); H=freqz(h,[1],w);43plot(w/pi,20*log10(abs(H))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 图 5-31 画出了所设计的 FIR 高通滤波器的频率特性。 由图可见， 用增加一个过渡点的方法， 可使线性相位 FIR 高通滤波器的阻带衰减 达到 40dB。 例 5-11 用频率取样设计一个线性相位数字微分器。 解：由例 5-6 知，理想数字微分器的频率响应为H DIF (e j? ) = j?, ? ≤ π由于微分器的幅度函数关于 ? =0 奇对称，故只能用 III 性或 IV 型线 性相位 FIR 滤波器。考虑到 III 型滤波器在 ? = π 的幅度值为零，故选 用 IV 型滤波器。由式（5-52）和微分器定义可知，在 0 ≤ m ≤ ( M + 1) / 2 范围内 Hd[m]值为H d [ m] = j 2πm ? j M +1m e ,0 ≤ m ≤ ( M + 1) / 2 M +1πmHd[m]在 ( M + 3) / 2 ≤ m ≤ M 范围内的值可由 Hd[m]在 1 ≤ m ≤ ( M = 1) / 2 范 围内的确定。根据以上分析就可得出下面的设计程序。 %program 5_6 频率抽样法设计数字微分器 M=15； m=0:(M+1)/2; Wm=2*pi*m/(M+1); Hd=j*Wm.*exp(-j*0.5*M*Wm); Hd=[Hd conj(f liplr(Hd(2(M+1)/2)))];44hd=real(ifft(Hd)); 图 5-32 画出了 M=15 时，FIR 数字微分器幅度函数和单位脉冲响应 h[k]. 5.5.3 优化设计的 MATLAB 实现离散加权平方误差准则设计 FIR由式（5-77）可知，一旦确定了矩阵 W、Q、C 及失量 d，就可 、 、 解方程得到最小加权平方误差的解 g。由 g 即可确定滤波器的 h[k]。 矩阵 W 和 Q 是对角阵， 其对角线的元素由 MATLAB 函数 diag 得出。 C 是 k × (J+1)矩阵，如 K=6，J=3 时有?? 1 ? ?? ? ? 2? ?? 3 ? ? ? [0 ?? 4 ? ?? 5 ? ? ? ?? 6 ? ? ? ?0 ?0 ? ?0 3]= ? ?0 ?0 ? ?0 ? ?1 ?2 ?3 ?4 ?5 ?6 2? 1 2? 2 2? 3 2? 4 2? 5 2? 6 3?1 ? 3? 2 ? ? 3? 3 ? ? 3? 4 ? 3? 5 ? ? 3? 6 ? ?1 2将上式作为函数 cos 的自变量即可获得 C 矩阵。 例 5-12 离 散 加 权 平 方 误 差 准 则 设 计 M=63 （ II 型 ） ，? p = 0.5πrad , ? s = 0.6πrad 的线性相位低通滤波器。解：由式（5-70）知，II 型线性相位低通滤波器的 Q（ ? ）为Q (?) = cos(? / 2)当 g[k]确定后，可由式（5-63）得到 II 型线性相位低通滤波器 的 h[k。完成上述设计过程的 MATLAB 程序如下： %progam 5_7 离散加权平方误差准则设计 FIR M=63;L=(M-1)/2; Wp=0.5*Ws==0.6*45%通带和阻带中的样本点数 Np=100;Ns=100; %通带和阻带的加权系数 weightp=1;weights=1; WpN=linspace(0,Wp,Np); WsN=linspace(Ws,pi-0.05,Ns); Wm=[WpN WsN]`; d=[ones(1,NP),zeros(1,Ns)]`; W=diag([weightp*ones(1,Np),weghts*ones(1,Ns)]); C=cos(Wm*[0:L]); Q=diag(cos(Wm*0.5)); g=(W*Q*C)\(W*d);g=g`; h=(g(L:-1:2)+g(L+1:-1:3))/4; h=[g(L+1)/4h g(1)*0.5+g(2)/4]; h=[h f liplr(h)]; w=linspace(0,pi,1000); mag=freqz(h,[1],w); plot(w/pi,20*log10(abs(mag))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 图 5-33（a）画出了通带和阻带的加权系数为 1，通带和阻带中 的离散频率点数为 100 时所设计的滤波器幅度响应， 由图可知滤波器46阻带衰减为 54dB。为增加滤波器阻带衰减，将阻带的加权系数增加 到 100，阻带中的离散频率点数增加到 200。所得的滤波器幅度响应 如图 5-33（b）所示，滤波器阻带衰减达到了 80dB。 积分加权平方误差准则设计 FIR MATLAB 信号处理工具箱提供的积分加权误差准则设计 FIR 滤 波器函数的基本调用形式为 h=firls(M,f,a,w) 函数 firls 的调用参数 M 表示滤波器的阶数。如 FIR 滤波器的 B 个频带分别为πf 1 ≤ ? ≤ πf 2,πf 3 ≤ ? ≤ πf 4 ,..., πf 2 B ?1 ≤ ? ≤ πf 2 B则 f 是一个有 2B 个元素的向量，其值为 f=[f1,f2,..., f 2 B ] a 也是一个有 2B 个元素的向量，表示滤波器在上述各频带边界的幅 度。 程序 firls 是在假设滤波器各频带的幅度是 2 个边界的幅度的线性 内插的条件下，计算出式（5-84）和式（5-85）的积分的解析表达式 而得出算法。w 是 B 个元素的向量，表示各频带的加权值。缺省时表 示各频带的加权值相同。例 5-13 用 积 分 加 权 平 方 误 差 准 则 设 计 M=63 （ II 型 ） ，? p = 0.5πrad , ? s = 0.6πrad 的线性相位低通滤波器。解： %program 5_8 积分加权平方误差准则设计 FIR M=63；47Fp＝0.5;Fs=0.6; h=firls(M,[0 Fp Fs 1],[1 1 0 0]); w=linspace(0,pi,1000); mag=freqz(h,[1],w); plot(w/pi,20*log10(abs(mag))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 作为比较图， 5-34 分别画出了积分加权平方误差准则和离散加权 平方误差准则设计的 FIR 滤波器的增益响应。 用离散加权平方误差准 则设计时通带取了 100 个点，阻带取了 200 个点。由图可以看出当离 散的数据点取得较多时， 两种方法设计出来的滤波器的增益响应几乎 相同。 等波纹 FIR 滤波器设计 MATLAB 信号工具箱提供的等波纹 FIR 滤波器设计函数的基本 调用形式为 [M,fo,ao,w]=remezord(f,a,dev) h=remez(M,fo,ao,w) 函数 remezord 估计滤波器阶数 M。函数 remez 完成等波纹 FIR 滤波 器设计。 函数 remezord 调用参数和 kaiserord 中的调用参数相同。FIR 滤 波器有 B 个频带时,f,a,dev 分别为 2B-2、B 和 B 个元素的向量。 函数 remezord 返回参数 M 表示滤波器的阶数。fo,ao 是 2B 个元48素的向量，分别表示 B 个频带的 2B 边界频率及幅度值。w 是 B 个 元素的向量， 表示各频带的加权值。 缺省时表示各频带的加权值相同。 函数 remez 的返回值为 M 阶等波纹 FIR 滤波器的系数。 例 5-14 设计满足下列指标的等波纹线性相位 FIR 低通滤波器。? p = 0.5πrad , ? s = 0.6πrad , δ p = δ s = 0.0017解： %program 5_9 等波纹 FIR 滤波器设计 Fp=0.5;Fs=0.6;ds=0.0017;dp= f=[Fp Fs];a=[1 0];dev=[dp ds]; [M,fo,ao,w]=remezord(f,a,dev); h=remez(M,fo,ao,w); w=linspace(0,pi,1000); mag=freqz(h,[1],w); hd=plot(w/pi,20*log10(abs(mag))); xlabel(`Normalized frequency`); ylabel(`Gain,dB`); 滤波器的阶数 M=59。图 5-35 画出了该滤波器的增益响应。例 5-15 设计满足下列指标的等波纹线性相位 FIR 带通滤波器。? s1 = 0.2πrad , ? p1 = 0.3πrad , ? p 2 = 0.6πrad , ? s1 = 0.7πrad , δ p = 0.1, δ s = 0.01解：只需将 program 5_9 的第 2 行和第 3 行改为 Fs1=0.2;Fp1=0.3;Fp2=0.6;Fs2=0.7;dp=0.1;ds=0.01; f=[Fs1 Fp1 Fp2 Fs2];a=[0 1 0];dev=[ds dp ds];49即可。由程序 remezord 确定的滤波器阶数 M=25，其增益响应如图 5-36（a）所示。由图可知滤波器在阻带的衰减没有满足指标。图 5-36 （b）为将阶数增加五阶后由 remez 得出的滤波器的幅度响应。 7.4 7.4.1 直接型 数字滤波器结构的 MATLAB 实现IIR 直接型结构由两个矢量 b 和 a 描述，b 包含{bn}系数，a 包含 {an}系数。它由 filter 函数实现，调用格式为 filter(b,a,x)。 FIR 直接型结构由包含{bn}系数的矢量 b 描述，它也由 filter 函 数实现，只是把矢量 a 置为 1，既调用格式为 filter(b,1,x)。 线性相位 FIR 数字滤波器的数字滤波器的系统函数形式相同， 只 是系数具有对称关系。因此，在 MATLAB 实现上，线性相位结构可 以由直接型结构实现。 7.4.2 级联型对于 IIR 级联型实现， 系统函数 H(z)通常可表示为有限个实系数 二阶有理分式之积，即H ( z) = Π b0i + b1i z ?1 + b2i z ?2 i =1 a + a z ?1 + a z ? 2 0i 1i 2iL（7-41）MATLAB 用一个 L × 6 矩阵 sos 来表示这个二阶分式? b01 ?b sos= ? 02 ? M ? ?b0 L b11 b12 M b1L b12 b22 M b2 L a 01 a 02 M a1L a11 a12 M a2L a12 ? a 22 ? ? M ? ? a2 L ?（7-42）若系统函数 H(z)是以零极点增益表示的，即H ( z) = k ( z ? z1 )( z ? z 2 ) L ( z ? z n ) ( z ? p1 )( z ? p 2 ) L ( z ? p n )（7-43）50则可用 zp2cos 函数将其转化式（7-41）表示的二阶分式之积，调用格 式为 sos=zp2sos(z,p,k) 将式（7-41）中的分母设为 1，即得 FIR 系统函数的有限个二阶因式 乘积表示。因此 FIR 系统函数的每一个二阶因式也可用 sos 来表示， 只需设式（7-42）中的 a 0i =1,其余 a 系数为零即可。 例 7-5 实现下面系统的级联型结构。5 ?1 2 ?2 z + z 3 3 H ( z) = 1 1 1 1 + z ?1 + z ?2 ? z ?3 6 3 6 3+解： %program 7_1 Cascade Realization of an IIR format rat num=[3,5/3,2/3]; den=[1,1/6,1/3,-1/6]; [z,p,k]=tf2zp(num,den); sos=zp2sos(z,p,k); 程序运行结果为 sos= 2/3 0 0 1 -/3 09/2 5/2 1 1 1/2 1/2 由此可得实系数有理分式之积的形式512 9 5 ?1 + z + z ?2 3 2 2 H ( z) 1 ?1 1 1 1 ? z 1 + z ?1 z + z ?2 3 2 2一般，系统函数的实数极点对应系数一阶有理分式，共轭复数极点对 应系数二阶有理分式。 7.4.3 并联型对于 IIR 并联型实现， 系统函数 H(z)通常可以表示为有限个实系 数一阶或二阶有理分式之和的形式。这种形式在 MATLAB 中可以通 过 residuez 函数实现，如例 7-5 所述系统，可以下面程序将 H(z)展开 成部分分式之和。 %program 7_2 parallel Realization of an IIR format rat num=[3,5/3,2/3]; den=conv([1,-1/3],[1,1/2,1/2]); [r,p,k]=residuez(num,den) 程序运行结果为 r=1/2+i, p=-1/4-506/765i, k=[ ] 该系统有一对共轭复数极点，利用[b1,a1]=residuez(R1,P1,0)语句可获 得其所对应的实数二阶分式的分子、分母多项式系数，其中 R1 为共 轭复数留数所构成的向量，P1 为共轭复数极点所构成的向量，b1,a1 为有理分式分子和分母多项式的系数向量。运行下面语句 1/2-i, -1/4+506/765i, 2 1/352R1=[r(1),r(2)]; P1=[p(1),p(2)]; [b1,a1]=residuez(R1,P1,0); 可得实系数二阶因子式的分子、分母多项式系数分别为 b1=1 a1=1 1 1/2 0 1/2由 b1 和 a1 及 r(3),p(3)可写出并联形式的系统函数，即H ( z) = 1 + z ?1 2 + 1 1 1 1 + z ?1 + z ? 2 1 ? z ?1 2 2 37.4.4 格型 全零点 AZ 和全极点 AP 系统的格型结构由包含 K 参数的行矢量 K 描述， 具有零点和极点系统的格型结构由包含 K 参数的行矢量和包 含 c 参数的行矢量 C 描述。 如果已知系统的直接型结构，那么 AZ 和 AP 系统格型结构的 K 参数可通过 ploy2rc 函数得到，调用格式为 K=ploy2rc(a) 如果已知 K 参数，通过 rc2ploy 函数可以得到直接型的多项式系数， 调用格式为 a=rc2ploy(K) 对于零极点系统的格型结构，如果已知系统的直接型结构，可通 过 tf2latc 函数求出 K 参数和 c 参数，调用格式为 [K,C]=tf2latc(num,den)53如果已知 K 参数和 c 参数，通过 latc2tf 函数可以得到直接型的分子 分母多项式系数，调用格式为 [num,den]=latc2tf(K,C) 全零点和全极点系统格型结构的 K 参数也可以通过 tf2latc 函数得到， 调用格式分别为 K=tf2latc(1,den)和 K=tf2latc(num) 格型结构可由 latcfilt 函数实现， 对于全零点、 全极点和零极点系 统的格型结构，其调用格式分别为 y=latcfilt(K,x),y=latcfilt(K,1,x)和 y=latcfilt(K,C,x)例 7-6解：用 MATLAB 实现例 7-4 所述系统的格型结构。%program 7_3 Lattice Realization of an IIR num=[3 5/3 2/3]; %分子多项式系数向量 den-=[1 1/6 1/3 -1/6]; %分母多项式系数向量 [K,C]=tf2latc(num,den); %求格型结构 K 参数和 C 参数x=[1 0 0 0]; %输入信号 y1=filter(num,den,x) %直接型结构实现 y2=latcfilt(K,C,x) MATLAB 计算结果为 K=0.4 -0.1667 C=2.3 0.667 0..7 -0.154%格型实现y2=3.7 -0.155
(&#39;噪声图像&#39;) 腐蚀和膨胀的 matlab 实现: map1=imread(&#39;p.jpg&#39;); [row,col,dep]=size(map1); %行,列,深度值 map=zeros(row,col); pixsum=row*col...?1 2 1 ? 2.3.1 特征向量及其归一化的 MATLAB 实现 MATLAB 中求矩阵特征值和特征向量的函数是 eig, 其调用的格式为[V,D]=eig(A), 其中,V 为特征向量...背景差分法MATLAB实现 - 程序 1 背景差分法 MATLAB 实现 function temp3 d=60; b=&#39;d\6\capfile3.avi&#39;; e=&#39;.bmp&#39;; for i=0...波前法及matlab实现_计算机软件及应用_IT/计算机_专业资料。有限元大型方程求解技术!有限元二维热传导波前法 MATLAB 程序 ? ? ? ? ? ? ? 二维热传导有限元 ...KMV与MATLAB实现_数学_自然科学_专业资料。目 录 一、课题背景 ... 1 1目 录 一、课题背景 ......matlab程序_计算机软件及应用_IT/计算机_专业资料。气象统计应用作业程序,各种程序clear clc fid=fopen(&#39;CO2.txt&#39;,&#39;r&#39;); tline=fgetl(fid); co=fscanf(fid,&#39;...matlab实现波特图 - matlab画波特图,用matlab画波特图。... matlab实现波特图_信息与通信_工程科技_专业资料。matlab画波特图,用matlab画波特图。 ...数字信号处理 Matlab 实现实例第1章 离散时间信号与系统 例 1-1 用 MATLAB 计算序列{-2 0 1 C1 3}和序列{1 2 0 -1}的离散卷积。 解 MATLAB 程序如...MATLAB程序大全_计算机软件及应用_IT/计算机_专业资料。MATLAB程序大全 1.全景图到穹景图 这个程序我最初是用 FreeImage 写的,这两天改成了 matlab,再不贴上来,...matlab程序大全答案_工学_高等教育_教育专区。matlab程序大全答案精彩的人生每天都是现场直播!iBO 频率特性类题目 1 一个系统的开环传递函数为 G ( s) ? k s...
All rights reserved Powered by
www.tceic.com
copyright &copyright 。文档资料库内容来自网络，如有侵犯请联系客服。