这份是我在网备滚卖络上找的：
（1）.周期图法： 思想：周期图法是为了得到功率谱估值，先取信号序列的离散傅里叶变换，然后取其幅频特性的平方并除以序列长度 N。由于序列 x(n)的离散傅里叶变换 X(k) 具有周期性，因而这种功率谱也具有周期性，常称为周期图。周期图是信号功率谱的一个有偏估值；而且，当信号序列的长度增大到无穷时，估值的方差不趋于零。因此，随着所取的信号序列长度的不同，所得到的周期图也不同，这种现象称为随机起伏。由于随机起伏大，使用周期图不能得到比较稳定的估值。
程序：
%首先，生成输入信号的程序为：
clear fs=20000 n=0:1/ fs: 0. 1N=lengt h( n)W=2000*pi%因方波频率 F=1000HZ 所以角频率 W=2000pi X1n=square( W*n) %方波信号 X2n=randn( 1, N)%白噪声信号 xn=X1n+0. 2*X2nsubplot (3, 1, 1) plot (n, xn) xlabel( ' n' ) ylabel('输入信仿逗号') %其次，开始用周期图法进行估计； clear all fs=20000n=0:1/fs:0.1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nsubplot(2,1,1) plot(n,xn)Nfft=256N=256%傅里叶变换的采样点数256 Pxx=abs(fft(xn,Nfft).^2)/Nf=(0:length(Pxx)-1)*fs/length(Pxx)subplot(2,1,2), plot(f,10*log10(Pxx)),%转成DB 单位
（2）.自相关函数法： 思想： 随机信号 x(n)的相关函数是在时间域内描述随机过程的重要特征。自相关函数是随机信号在不同时刻的值之间的依赖性的量度，是一个很有用的统计平均量，其定义为自相关函数 (1) 式中E【·】表示数学期望，*表示共轭值,m 为时间滞后数。在随机信号处理中，自相关函数可以用来检测淹没在随机噪声干扰中的信号，随机信号的自功率谱等于它的自相关函数的傅里叶变换。因此，通过自相关估计可求得信号的功率谱。利用计算机计算自相关估值有两种方法。一种是直接方法，先计算出随机信号和它的滞后序列的乘积，再取其平均值即得相关函数的估计值。另一种是间接方法，先用快速变换算法计算随机序列的功率谱密度，再作反变换计算出相关函数。
程序： n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(2) subplot(3,1,1) plot(n,xn)%输入信号 axis([0 0.01 -5 5]) m =-100:100 [r,lag]=xcorr(xn,100,'biased')%求XN 的自相关函数R，biased 为有偏估计,lag 为R 的序列号 subplot(3,1,2) hndl=stem(m,r)%绘制离散图，分布点从-100—+100 set(hndl,'Marker','.') set(hndl,'MarkerSize',2)ylabel('自相关函数R(m)') %利用间接法计算功率谱 k=0:1000%取1000 个点 w=(pi/500)*kM=k/500X=r*(exp(-li*pi/500).^(m'*k))%对R 求傅里叶变换,li 就是j magX=abs(X)subplot(3,1,3) plot(M,10*log10(magX))xlabel('功率谱的改进直接法估计')
（3）.自协方差法： 思想：在实际中,有时用随机信号在两个不同时刻t 1 ,t 2 的取值X(t 1 )和X(t 2 )之间的二阶混合中心矩来描述随机信号 X(t)在任意两个时刻取值起伏变化的相依程度, 即自协方差。自协方差函数与自相关函数描述的特性是一致的，所以其原理与相关函数法近似。
clear allfs=20000n=0:1/fs:0.1P=2000*piy=square(P*n)xn=y+0.2*randn(size(n))%绘制信号波形 figure(3) subplot(2,1,1) plot(n,xn) xlabel('时间（s）') ylabel('幅度') title('y+randn(size(n))') ymax_xn=max(xn)+0.2ymin_xn=min(xn)-0.2axis([0 0.3 ymin_xn ymax_xn]) %使用协方差法估计序列功率谱 p=floor(length(xn)/3)+1nfft=1024[xpsd,f]=pcov(xn,p,nfft,fs,'half')%绘制功率谱估计 pmax=max(xpsd)xpsd=xpsd/pmaxxpsd=10*log10(xpsd+0.000001)subplot(2,1,2) plot(f,xpsd) title('基于协方差的功率谱估计') ylabel('功率谱估计（db）') xlabel('频率') grid onymin=min(xpsd)-2ymax=max(xpsd)+2axis([0 fs/2 ymin ymax])
（4）.最大熵法 思想：上网查的最大熵法原理是取一组时间序列，使其自相关函数与一组已知数据的自相关函数相同，同时使已知自相关函数以外的部分的随机性最强，以所取时间序列的谱作为已知数据的谱估值。它等效于根据使随机过程的熵为最大的原则,利用 N 个已知的自相关函数值来外推其他未知的自相关函数值所得到的功率谱。最大熵法功率谱估值是一种可获得高分辨率的非线性谱估值方法，特别适用于数据长度较短的情况。最大熵法谱估值对未知数据的假定 ，一个平稳的随机序列，可以用周期图法对其功率谱进行估值。这种估值方法隐含着假定未知数据是已知数据的周期性重复。现有的线性谱估计方法是假定未知数据的自相关函数值为零，这种人为假定带来的误差较大。最大熵法是利用已知的自相关函数值来外推未知的自相关函数值，去除了对未知数据的人为假定，从而使谱估计的结果更为合理。熵在信息论中是信息的度量，事件越不确定，其信息量越大，熵也越大。对于上述问题来说，对随机过程的未知的自相关函数值，除了从已知的自相关函数值得到有关它的信息以外，没有其他的先验知识。因而，在外推时，不希望加以其他任何新的限制，亦即使之“最不确定”。换言之，就是使随机过程的熵最大。
程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(5) subplot(3,1,1) plot(n,xn)asis([0 0.01 -5 5]) Nfft=256%分段长度256 [Pxx,f]=pmem(xn,14,Nfft,fs)%调用最大熵函数pmem，滤波器阶数14 subplot(2,1,2), plot(f,10*log10(Pxx)), title(' 最大熵法，滤波器14'),xlabel('频率'),ylabel('功率谱db')
（5）.最大似然法： 思想：最大似然法原理是让信号通过一个滤波器，选择滤波器的参数使所关心的频率的正弦波信号能够不失真地通过，同时，使所有其他频率的正弦波通过这个滤波器后输出的均方值最小。在这个条件下，信号经过这个滤波器后输出的均方值就作为其最大似然法功率谱估值。可以证明,如果信号x 是由一个确定性信号S 加上一个高斯白噪声n 所组成，则上述滤波器的输出是信号S 的最大似然估值。如果n 不是高斯噪声，则上述滤波器的输出是信号S 的最小方差的线性的无偏估值，而且它能得到比使用固定的窗口函数的周期图法更高的分辨率。
程序： fs=20000n=0:1/fs:1N=length(n)W=2000*pix1n=square(W*n)x2n=randn(1,N)xn=x1n+0.2*x2nfigure(4) subplot(3,1,1) plot(n,xn)axis([0 0.01 -5 5]) %估计自相关函数 m=-500:500[r,lag]=xcorr(xn,500,'biased')R=[r(501) r(502) r(503) r(504)r(500) r(501) r(502) r(503)r(499) r(500) r(501) r(502)r(498) r(499) r(500) r(501)][V,D]=eig(R)V3=[V(1,3),V(2,3),V(3,3),V(4,3)].'V3=[V(1,4),V(2,4),V(3,4),V(4,4)].'p=0:3wm=[0:0.002*pi:2*pi]B=[(exp(-li)).^(wm'*p)]%li 就是虚数单位j A=B%最小方差功率谱估计 z=A*inv(R)*A'Z=diag(z')pmv=1./Zsubplot(2,1,2) plot(wm/pi,pmv)title('基于最大似然的功率谱估计') ylabel('功率谱幅度（db）') xlabel('角度频率w/pi')
自己分下段就能运行了。
直接法：直接法又称周期图法，它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。
间接法：
间接法先由序列x(n)估计出自相关函数R(n)，然后对R(n)进行傅立叶变换，便得到x(n)的功率谱估计。
改进的直接法：
对于直接法的功率谱估计，当数据长度N太大时，谱曲线起伏加剧，若N太小，谱的分辨率又不好，因此需要改进。
Bartlett平均周期派搏运图的方法是将N点的有限长序列x(n)分段求周期图再平均。
功率谱估计在现代信号处理中是一个很重要的课题，涉及的问题很多。在这里，结合matlab，做一个粗略介绍。功率谱估计可以分为经典谱估计方法与现代谱估计方法。经典谱估计中最简单的就是周期图法，又分为直接法与间接法。直接法先取N点数据的傅里叶变换（即频谱），然后取频谱与其共轭的乘积银昌，就得到功率谱的估计；间接法先计算N点样本数据的自相关函数，然后取自相关函数的傅里叶变换，即得到功率谱的估计.都可以编程实现，很简单。在matlab中，周期图法可以用函数periodogram实现。
但是周期图法估计出的功率谱不够精细，分辨率比较低。因此需要对周期图法进行修正，可以将信号序列x(n)分为n个不相重叠的小段，分别用周期图法进行谱估计，然后将这n段数据估计的结果的平均值作为整段数据功率谱估计的结果。还可以将信号序列x(n)重叠分段，分别计算功率谱，再计算平均值作为整段数据的功率谱估计。
这2种称为分段平均周期图法，一般后者比前者效果好。加窗平均周期图法是对分段平均周期图法的改进，即在数据分段后，对每段数据加一个非矩形窗进行预处理，然后在按分段平均周期图法估计功率谱。相对于分段平均周期图法，加窗平均周期图法可以减小频率泄漏，增加频峰的宽度。welch法就是利用改进的平均周期图法估计估计随机信号的功率谱，它采用信号分段重叠尘梁，加窗，FFT等技术来计算功率谱。与周期图法比较，welch法可以改善估计谱曲线的光滑性，大大提高谱估计的分辨率。
现代谱估计主要针对经典谱估计分辨率低和方差性不好提出的，可以极大的提高估计的分辨率和平滑性。可以分为参数模型谱估计和非参数模型谱估计。参数模型谱估计有AR模型，MA模型，ARMA模型等；非参数模型谱估计有最小方差法和MUSIC法等。由于涉及的问题太多，这里不再详述，可以参考有关资料。
matlab中，现代谱估计的很多方法都可以实现。music方法用pmusic命令实现；pburg函数利用burg法实现功率谱估计；pyulear函数利用yule-walker算法实现功率谱估计等等。
另外，sptool工具箱也具有功率谱估计的功能。窗口化的 *** 作界面很方便，而且有多种方法可以选择。
1.1 实验目的1．了解数字信号处理系统的一般构成；
2．掌握奈奎斯特抽样定理。
1.2 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
1.3 实验原理
一个典型的DSP系统除了数字信号处理部分外，还包括A/D和D/A两部分。这是因为自然界的信号，如声音、图像等大多是模拟信号，因此需要将其数字化后进行数字信号处理，模拟信号的数字化即称为A/D转换。数字信号处理后的数据可能需还原为模拟信号，这就需要进行D/A转换。一个仅包括A/D和D/A两部分的简化数字信号处理系统功能如图1所示。
A/D转换包括三个紧密相关的过程，即抽样、量化和编码。A/D转换中需解决的以下几个重要问题：抽样后输出信号中还有没有原始信号伍铅烂的信息？如果有能不能把它取出来？抽样频率应该如何选择？
奈奎斯特抽样定理（即低通信号的均匀抽样定理）告诉我们，一个频带限制在0至fx以内的低通信号x(t)，如果以fs≥2fx的抽样速率进行均匀抽样，则x(t)可以由抽样后的信号xs(t)完全地确定，即xs(t)包含有x(t)的成分，可以通过适当的低通滤波器不失真地恢复出x(t)。最小抽样速率fs=2fx称为奈奎斯特速率。
低通
译码
编码
量化
抽样
输入信号 样点输出 滤波输出
A/D（模数转换） D/A（数模转换）
图1 低通采样定理演示
为方便实现，实验中更换了一种表现形式，即抽样频率固定（10KHz），通过改变输入模拟信号的频率来展示低通抽样定理。我们可以通过研究抽样频率和模拟信号最高频率分量的频率之间的关系，来验证低通抽样定理。
1.4 实验内容
1．软件仿真实验：编写并调试MATLAB程序，分析有关参数，记录有关波形。
2．硬件实验：输入不同频率的正弦信号，观察采样时钟波形、输入信号波形、样点输出波形和滤波输出波形。
1.5 MATLAB参考程序和仿真内容
%*******************************************************************%
%f—余弦信号的频率
% M—基2 FFT幂次数 N=2^M为采样点数，这样取值是为了便于作基2的FFT分析
%2. 采样频率Fs
%*******************************************************************%
function samples(f,Fs,M)
N=2^M% fft点数=取样总点数
Ts=1/Fs% 取样时间间隔
T=N*Ts% 取样总时间=取样总点数*取样时间间隔
n=0:N-1
t=n*Ts
Xn=cos(2*f*pi*t)
subplot(2,1,1)
stem(t,Xn)
axis([0 T 1.1*min(Xn) 1.1*max(Xn)])
xlabel('t -->')
ylabel('Xn')
Xk=abs(fft(Xn,N))
subplot(2,1,2)
stem(n,Xk)
axis([0 N 1.1*min(Xk) 1.1*max(Xk)])
xlabel('frequency -->')
ylabel('!Xk!')
%*******************************************************************%
假如有一个1Hz的腔漏余弦信号y=cos(2*π*t)，对其用4Hz的采样频率进行采样，共采样32点，只需执行samples(1,4,5)，即可得到仿真结果。
软件仿真实验内容如下表所示：
仿真参数
f
Fs
Wo(计算)
Xn（图形）
Xk（图形）
（1，4，5）
另外记录图形，并标图号
（1，8，5）
（2，8，6）
自 选
1.6 硬件实验步骤
本实验箱采样频率fs固定为10KHz，低通滤波器的截止频率约为4.5KHz。
1、用低频信号源产生正弦信号，正弦信号源频率f自定，并将其接至2TP2（模拟输入）端，将示波器通道一探头接至2TP6（采样时钟）端观察采样时钟波形，示波器通道二探头接至2TP2观察并记录输入信号波形。
2、将示波器通道二探头接至2TP3观激铅察并记录样点输出波形。
3、将示波器通道二探头接至2TP4观察并记录滤波输出波形。
4、根据采样定理，分f=fs /8、f=fs/4、f=fs/2等3种情况更改正弦信号频率，重复步骤2至步骤3。
5、用低频信号源产生方波信号，重复步骤1至步骤4。
1.7 思考题
1、 讨论在仿真实验中所计算的数字域频率Wo和Xk的图形中非零谱线位置之间的对应关系。
2、 讨论在仿真实验中自选参数的意义。
3、将在2TP2端加方波信号后的恢复波形，与相同频率的正弦信号的恢复波形相比，能够得出哪些结论？
2 FFT频谱分析实验
2.1 实验目的
1．通过实验加深对快速傅立叶变换（FFT）基本原理的理解。
2．了解FFT点数与频谱分辨率的关系，以及两种加长序列FFT与原序列FFT的关系。
2.2 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
2.3 实验原理
离散傅里叶变换（DFT）和卷积是信号处理中两个最基本也是最常用的运算，它们涉及到信号与系统的分析与综合这一广泛的信号处理领域。实际上卷积与DFT之间有着互通的联系：卷积可化为DFT来实现，其它的许多算法，如相关、滤波和谱估计等都可化为DFT来实现，DFT也可化为卷积来实现。
对N点序列x(n)，其DFT变换对定义为：
在DFT运算中包含大量的重复运算。FFT算法利用了蝶形因子WN的周期性和对称性，从而加快了运算的速度。FFT算法将长序列的DFT分解为短序列的DFT。N点的DFT先分解为2个N/2点的DFT，每个N/2点的DFT又分解为2个N/4点的DFT。按照此规律，最小变换的点数即所谓的“基数（radix）。”因此，基数为2的FFT算法的最小变换（或称蝶形）是2点DFT。一般地，对N点FFT，对应于N个输入样值，有N个频域样值与之对应。一般而言，FFT算法可以分为时间抽取（DIT）FFT和频率抽取（DIF）两大类。
在实际计算中，可以采用在原来序列后面补0的加长方法来提高FFT的分辨率；可以采用在原来序列后面重复的加长方法来增加FFT的幅度。
2.4 实验内容
1．软件仿真实验：分别观察并记录正弦序列、方波序列及改变FFT的点数后的频谱；分别观察并记录正弦序列、方波序列及2种加长序列等信号的频谱。
2．硬件实验：分别观察并记录正弦信号、方波信号及改变FFT的点数后的频谱。
2.5 MATLAB参考程序和仿真内容
%*******************************************************************%
function[x]=ffts(mode,M)
Nfft=2^M
x=zeros(1,Nfft)%定义一个长度为Nfft的一维全0数组
if mode= =1 for n=0:Nfft-1 x(n+1)=sin(2*pi*n/Nfft)end
end %定义一个长度为Nfft的单周期正弦序列
if mode= =2 for n=0:Nfft-1 x(n+1)=sin(4*pi*n/Nfft)end
end %定义一个长度为Nfft的双周期正弦序列
if mode= =3 for n=0:Nfft/2-1 x(n+1)=sin(4*pi*n/Nfft)end
end %定义一个长度为Nfft/2的正弦序列，后面一半为0序列。
if mode= =4 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft)end
end
if mode= =5 for n=0:Nfft-1 x(n+1)=square(2*pi*n/Nfft)end
end
if mode= =6 for n=0:Nfft/2-1 x(n+1)=square(4*pi*n/Nfft)end
end
n=0:Nfft-1
subplot(2,1,1)
stem(n,x)
axis([0 Nfft-1 1.1*min(x) 1.1*max(x)])
xlabel('Points-->')
ylabel('x(n)')
X=abs(fft(x,Nfft))
subplot(2,1,2)
stem(n,X)
axis([0 Nfft-1 1.1*min(X) 1.1*max(X)])
xlabel('frequency-->')
ylabel('!X(k)!')
%*******************************************************************%
假设需观察方波信号的频谱，对一个周期的方波信号作32点的FFT，则只需在MATLAB的命令窗口下键入：[x]=ffts(21,5) ，程序进行模拟，并且输出FFT的结果。
关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个以上的例子进行记录。例如要观察后面补0的加长方法来提高FFT的分辨率的现象，可以仿真ffts(4,5)和ffts(6,6)两个例子。
2.6 硬件实验步骤
1．将低频信号源输出加到实验箱模拟通道1输入端，将示波器探头接至模拟通道1输出端。
2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“FFT频谱分析”子菜单，出现显示FFT频谱分析功能提示信息的窗口。
3．用低频信号产生器产生一个1KHz的正弦信号。
4．选择FFT频谱分析与显示的点数为64点，开始进行FFT运算。此后，计算机将周期性地取回DSP运算后的FFT数据并绘图显示
5．改信号源频率，观察并记录频谱图的变化。
6．选择FFT的点数为128点，观察并记录频谱图的变化。
7．更改正弦信号的频率，重复步骤4 ~步骤6。
8．用低频信号产生器产生一个1KHz的方波信号，重复步骤4 ~步骤7。注意：应根据实验箱采样频率fs为10KHz和方波信号的频带宽度选择方波信号的频率。
本硬件实验要进行两种信号，每个信号两种频率，每个信号两种点数等共8次具体实验内容，性质能够体现实验要求的4个以上的例子进行记录。
2.7 思考题
1．对同一个信号，不同点数FFT观察到的频谱图有何区别？
2．序列加长后FFT与原序列FFT的关系是什么，试推导其中一种关系。
3．用傅立叶级数理论，试说明正弦信号频谱和方波信号频谱之间的关系。
3 IIR滤波器设计实验
3.1 实验目的
1．通过实验加深对IIR滤波器基本原理的理解。
2．学习编写IIR滤波器的MATLAB仿真程序。
3.2 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
3.3 实验原理
IIR滤波器有以下几个特点：
1．IIR数字滤波器的系统函数可以写成封闭函数的形式。
2．IIR数字滤波器采用递归型结构，即结构上带有反馈环路。IIR滤波器运算结构通常由延时、乘以系数和相加等基本运算组成，可以组合成直接型、正准型、级联型、并联型四种结构形式，都具有反馈回路。由于运算中的舍入处理，使误差不断累积，有时会产生微弱的寄生振荡。
3．IIR数字滤波器在设计上可以借助成熟的模拟滤波器的成果，如巴特沃斯、契比雪夫和椭圆滤波器等，有现成的设计数据或图表可查，其设计工作量比较小，对计算工具的要求不高。在设计一个IIR数字滤波器时，我们根据指标先写出模拟滤波器的公式，然后通过一定的变换，将模拟滤波器的公式转换成数字滤波器的公式。
4．IIR数字滤波器的相位特性不好控制，对相位要求较高时，需加相位校准网络。
在MATLAB下设计IIR滤波器可使用Butterworth函数设计出巴特沃斯滤波器，使用Cheby1函数设计出契比雪夫I型滤波器，使用Cheby2设计出契比雪夫II型滤波器，使用ellipord函数设计出椭圆滤波器。下面主要介绍前两个函数的使用。
与FIR滤波器的设计不同，IIR滤波器设计时的阶数不是由设计者指定，而是根据设计者输入的各个滤波器参数（截止频率、通带滤纹、阻带衰减等），由软件设计出满足这些参数的最低滤波器阶数。在MATLAB下设计不同类型IIR滤波器均有与之对应的函数用于阶数的选择。
一、巴特沃斯IIR滤波器的设计
在MATLAB下，设计巴特沃斯IIR滤波器可使用butter函数。
Butter函数可设计低通、高通、带通和带阻的数字和模拟IIR滤波器，其特性为使通带内的幅度响应最大限度地平坦，但同时损失截止频率处的下降斜度。在期望通带平滑的情况下，可使用butter函数。
butter函数的用法为：
[b,a]=butter(n,Wn,/ftype/)
其中n代表滤波器阶数，Wn代表滤波器的截止频率，这两个参数可使用buttord函数来确定。buttord函数可在给定滤波器性能的情况下，求出巴特沃斯滤波器的最小阶数n，同时给出对应的截止频率Wn。buttord函数的用法为：
[n,Wn]= buttord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。
不同类型（高通、低通、带通和带阻）滤波器对应的Wp和Ws值遵循以下规则：
1．高通滤波器：Wp和Ws为一元矢量且Wp>Ws；
2．低通滤波器：Wp和Ws为一元矢量且Wp<Ws；
3．带通滤波器：Wp和Ws为二元矢量且Wp<Ws，如Wp=[0.2,0.7],Ws=[0.1,0.8]
4．带阻滤波器：Wp和Ws为二元矢量且Wp>Ws，如Wp=[0.1,0.8],Ws=[0.2,0.7]。
二、契比雪夫I型IIR滤波器的设计
在期望通带下降斜率大的场合，应使用椭圆滤波器或契比雪夫滤波器。在MATLAB下可使用cheby1函数设计出契比雪夫I型IIR滤波器。
cheby1函数可设计低通、高通、带通和带阻契比雪夫I型滤IIR波器，其通带内为等波纹，阻带内为单调。契比雪夫I型的下降斜度比II型大，但其代价是通带内波纹较大。
cheby1函数的用法为：
[b,a]=cheby1(n,Rp,Wn,/ftype/)
在使用cheby1函数设计IIR滤波器之前，可使用cheblord函数求出滤波器阶数n和截止频率Wn。cheblord函数可在给定滤波器性能的情况下，选择契比雪夫I型滤波器的最小阶和截止频率Wn。
cheblord函数的用法为：
[n,Wn]=cheblord(Wp,Ws,Rp,Rs)
其中Wp和Ws分别是通带和阻带的拐角频率（截止频率），其取值范围为0至1之间。当其值为1时代表采样频率的一半。Rp和Rs分别是通带和阻带区的波纹系数。
3.4 实验内容
1．软件仿真实验：编写并调试MATLAB程序，选择不同形式，不同类型的4种滤波器进行仿真，记录幅频和相频特性，对比巴特沃斯滤波器和契比雪夫滤波器。
2．硬件实验：设计IIR滤波器，在计算机上观察冲激响应、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。
3.5 MATLAB参考程序和仿真内容
%*******************************************************************%
%mode: 1--巴特沃斯低通；2--巴特沃斯高通；3--巴特沃斯带通；4--巴特沃斯带阻
% 5--契比雪夫低通；6--契比雪夫高通；7--契比雪夫带通；8--契比雪夫带阻
%fp1,fp2： 通带截止频率，当高通或低通时只有fp1有效
%fs1, fs2： 阻带截止频率，当高通或低通时只有fs1有效
%rp: 通带波纹系数
%as: 阻带衰减系数
%sample: 采样率
%h: 返回设计好的滤波器系数
%*******************************************************************%
function[b,a]=iirfilt(mode,fp1,fp2,fs1,fs2,rp,as,sample)
wp1=2*fp1/samplewp2=2*fp2/sample
ws1=2*fs1/samplews2=2*fs2/sample
%得到巴特沃斯滤波器的最小阶数N和3bd频率wn
if mode<3[N,wn]=buttord(wp1,ws1,rp,as)
elseif mode<5[N,wn]=buttord([wp1 wp2],[ws1 ws2],rp,as)
%得到契比雪夫滤波器的最小阶数N和3bd频率wn
elseif mode<7[N,wn]=cheb1ord(wp1,ws1,rp,as)
else[N,wn]=cheblord([wp1 wp2],[ws1 ws2],rp,as)
end
%得到滤波器系数的分子b和分母a
if mode= =1[b,a]=butter(N,wn)end
if mode= =2[b,a]=butter(N,wn,/high/)end
if mode= =3[b,a]=butter(N,wn)end
if mode= =4[b,a]=butter(N,wn,/stop/)end
if mode= =5[b,a]=cheby1(N,rp,wn)end
if mode= =6[b,a]=cheby1(N,rp,wn,/high/)end
if mode= =7[b,a]=cheby1(N,rp,wn)end
if mode= =8[b,a]=cheby1(N,rp,wn,/stop/)end
set(gcf,/menubar/,menubar)
freq_response=freqz(b,a)
magnitude=20*log10(abs(freq_response))
m=0:511
f=m*sample/(2*511)
subplot(3,1,1)plot(f,magnitude)grid%幅频特性
axis([0 sample/2 1.1*min(magnitude) 1.1*max(magnitude)])
ylabel('Magnitude')xlabel('Frequency-->')
phase=angle(freq_response)
subplot(3,1,2)plot(f,phase)grid%相频特性
axis([0 sample/2 1.1*min(phase) 1.1*max(phase)])
ylabel('Phase')xlabel('Frequency-->')
h=impz(b,a,32)%32点的单位函数响应
t=1:32
subplot(3,1,3)stem(t,h)grid
axis([0 32 1.2*min(h) 1.1*max(h)])
ylabel('h(n)')xlabel('n-->')
%*******************************************************************%
假设需设计一个巴特沃斯低通IIR滤波器，通带截止频率为2KHz，阻带截止频率为3KHz，通带波纹系数为1，阻带衰减系数为20，采样频率为10KHz，则只需在MATLAB的命令窗口下键入：
[b,a]=iirfilt(1,2000,3000,2400,2600,1,20,10000)
程序进行模拟，并且按照如下顺序输出数字滤波器系统函数
的系数
b= b0 b1 ……bn
a= a0 a1 ……an
关于软件仿真实验内容，建议在完成大量仿真例子的基础上，选择能够体现实验要求的4个例子进行记录，系统函数只要记录系统的阶数。
3.6 硬件实验步骤
1．根据实验箱采样频率fs为10KHz的条件，用低频信号发生器产生一个频率合适的低频正弦信号，将其加到实验箱模拟通道1输入端，将示波器通道1探头接至模拟通道1输入端，通道2探头接至模拟通道2输出端。
2．在保证实验箱正确加电且串口电缆连接正常的情况下，运行数字信号处理与DSP应用实验开发软件，在“数字信号处理实验”菜单下选择“IIR滤波器”子菜单，出现提示信息。
3．输入滤波器类型、滤波器截止频率等参数后，分别点击“幅频特性”和“相频特性”按钮，在窗口右侧观察IIR滤波器的幅频特性和相频特性。此时提示信息将消失，如需查看提示信息，可点击“设计说明”按钮。
4．点击“下载实现”按钮，IIR滤波器开始工作，此时窗口右侧将显示IIR滤波器的幅频特性。
5．根据输入滤波器类型，更改低频信号源的频率，观察示波器上输入输出波形幅度的变化情况，测量IIR滤波器的幅频响应特性，看其是否与设计的幅频特性一致。
6．更改滤波器类型、滤波器截止频率等参数（共4种），重复步骤3至步骤5。所选择的例子参数最好和MATLAB仿真程序的例子一样。
7．用低频信号产生器产生一个500Hz的方波信号，分别设计3种滤波器，完成如下表要求的功能，并且记录参数和波形。
功 能
滤波器类型
参 数
输出波形
fp1
fp2
fs1
fs2
通过3次及以下次数的谐波
另外记录图形，并标图号
滤除5次及以下次数的谐波
通过3次到5次的谐波
3.7 思考题
1．在实验箱采样频率fs固定为10KHz的条件下，要观察方波信号频带宽度内的各个谐波分量，方波信号的频率最高不能超过多少，为什么？
2．硬件实验内容7中输出信号各个谐波分量，与原来方波信号同样谐波分量相比，有没有发生失真？主要发生了什么类型的失真？为什么？
4 窗函数法FIR滤波器设计实验
4.1 实验目的
1．通过实验加深对FIR滤波器基本原理的理解。
2．学习使用窗函数法设计FIR滤波器，了解窗函数的形式和长度对滤波器性能的影响。
4.2 实验仪器
1．YBLD智能综合信号源测试仪 1台
2．双踪示波器 1台
3．MCOM－TG305数字信号处理与现代通信技术实验箱 1台
4．PC机（装有MATLAB、MCOM－TG305配套实验软件） 1台
4.3 实验原理
数字滤波器的设计是数字信号处理中的一个重要内容。数字滤波器设计包括FIR（有限单位脉冲响应）滤波器与IIR（无限单位脉冲响应）滤波器两种。
与IIR滤波器相比，FIR滤波器在保证幅度特性满足技术要求的同时，很容易做到严格的线性相位特性。设FIR滤波器单位脉冲响应h(n)长度为N，其系统函数H（z）为：
H（z）是z－1的N－1次多项式，它在z平面上有N－1个零点，原点z=0是N－1阶重极点，因此H（z）是永远稳定的。稳定和线性相位特性是FIR滤波器突出的优点。
FIR滤波器的设计任务是选择有限长度的h(n)。使传输函数H( )满足技术要求。FIR滤波器的设计方法有多种，如窗函数法、频率采样法及其它各种优化设计方法，本实验介绍窗函数法的FIR滤波器设计。
窗函数法是使用矩形窗、三角窗、巴特利特窗、汉明窗、汉宁窗和布莱克曼窗等设计出标准响应的高通、低通、带通和带阻FIR滤波器。
一、firl函数的使用
在MATLAB下设计标准响应FIR滤波器可使用firl函数。firl函数以经典方法实现加窗线性相位FIR滤波器设计，它可以设计出标准的低通、带通、高通和带阻滤波器。firl函数的用法为：
b=firl(n,Wn,/ftype/,Window)
各个参数的含义如下：
b—滤波器系数。对于一个n阶的FIR滤波器，其n+1个滤波器系数可表示为：b(z)=b(1)+b(2)z－1+…+b(n+1)z－n。
n—滤波器阶数。
Wn—截止频率，0≤Wn≤1，Wn=1对应于采样频率的一半。当设计带通和带阻滤波器时，Wn=[W1 W2],W1≤ω≤W2。
ftype—当指定ftype时，可设计高通和带阻滤波器。Ftype=high时，设计高通FIR滤波器；ftype=stop时设计带阻FIR滤波器。低通和带通FIR滤波器无需输入ftype参数。
Window—窗函数。窗函数的长度应等于FIR滤波器系数个数，即阶数n+1。
二、窗函数的使用
在MATLAB下，这些窗函数分别为：
1．矩形窗：w=boxcar(n)，产生一个n点的矩形窗函数。
2．三角窗：w=triang(n)，产生一个n点的三角窗函数。
当n为奇数时，三角窗系数为w(k)=
当n为偶数时，三角窗系数为w(k)=
3．巴特利特窗：w=Bartlett(n)，产生一个n点的巴特利特窗函数。
巴特利特窗系数为w(k)=
巴特利特窗与三角窗非常相似。巴特利特窗在取样点1和n上总以零结束，而三角窗在这些点上并不为零。实际上，当n为奇数时bartlett(n)的中心n－2个点等效于triang(n－2)。
4．汉明窗：w=hamming(n)，产生一个n点的汉明窗函数。
汉明窗系数为w(k+1)=0.54－0.46cos（ ） k=0,…,n－1
5．汉宁窗：w=hanning(n)，产生一个n点的汉宁窗函数。
汉宁窗系数为w(k)=0.5[1－cos( )] k=1,…,n
6．布莱克曼窗：w=Blackman(n),产生一个n点的布莱克曼窗函数。
布莱克曼窗系数为w(k)=0.42－0.5cos(2π )+0.8cos（4π ）] k=1,…,n
与等长度的汉明窗和汉宁窗相比，布莱克曼窗的主瓣稍宽，旁瓣稍低。
7．凯泽窗：w=Kaiser(n,beta)，产生一个n点的凯泽窗数，其中beta为影响窗函数旁瓣的β参数，其最小的旁瓣抑制α与β的关系为:
0.1102(α－0.87) α>50
β= 0.5842(α－21)0.4+0.07886(α－21) 21≤α≤50
0 α<21
增加β可使主瓣变宽，旁瓣的幅度降低。
8．契比雪夫窗：w=chebwin(n,r)产生一个n点的契比雪夫窗函数。其傅里叶变换后的旁瓣波纹低于主瓣r个db数。
4.4 实验内容
1．软件仿真实验：编写并调试MATLAB程序，观察不同窗，不同类型滤波器不同点数等共4种FIR滤波器的h(n)，并记录幅频特性和相频特性。
2．硬件实验：用窗函数法设计标准响应的FIR滤波器，在计算机上观察窗函数幅频特性、幅频特性和相频特性，然后下载到实验箱。用示波器观察输入输出波形，测试滤波器的幅频响应特性。
4.5 MATLAB参考程序和仿真内容
%*******************************************************************%
%mode: 模式（1--高通；2--低通；3--带通；4--带阻）
%n: 阶数，加窗的点数为阶数加1
%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率
%fs: 带通和带阻时指示上限频率
%window:加窗（1--矩形窗；2--三角窗；3--巴特利特窗；4--汉明窗；
% 5--汉宁窗；6--布莱克曼窗；7--凯泽窗；8--契比雪夫窗）
%r: 代表加chebyshev窗的r值和加kaiser窗时的beta值
%sample: 采样率
%h: 返回设计好的FIR滤波器系数
%*******************************************************************%
%mode: 模式（1--高通；2--低通；3--带通；4--带阻）
%n: 阶数，加窗的点数为阶数加1
%fp: 高通和低通时指示截止频率，带通和带阻时指示下限频率
%fs:
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）