DSPA指导书正文.doc

DSPA指导书正文 实验1矩阵、序列的运算与显示 一、实验目的掌握用MATLAB表示离散时间信号，以及它们的运算和显示。 二、实验内容与要求1用MATLAB产生并画出下列序列的样本。 ax1(n)?m?0?m?1?n?2m?n?2m?1?,0?n?25n102bx2(n)?nu?n?5?u?n?6?10?n?20?0.5?u?n?4?u?n?10?cx3(n)?0.9?cos(0.2?n?/3),0?n?20dx4(n)?10cos(0.0008?n)?w?n?,0?n?100，式中w?n?是在【-1，1】之间均匀分2n布的随机序列，你如何表征这个序列？?5?n?.,1,2,3,2,1,2,3,2,.为周期的，画出5个周期。 ex? 三、实验所用部分函数如下1.单位冲激序列（信号）生成函数impseqx,n=impseq(n0,n1,n2)2.阶跃序列（信号）生成函数stepseqx,n=stepseq(n0,n1,n2)3.序列（信号）相加函数sigaddy,n=sigadd(x1,n1,x2,n2)以上为MATLAB没有，需外加入的函数（将相应函数拷贝到自己当前目录下）4.正（余）弦生成函数sin、cos y=sin(x),y=cos(x)(注意:x以弧度为单位)5.随机序列生成函数rand，用法如Y=rand(n)生成nn阶的均匀分布随机阵；Y=rand(m,n)生成mn阶的随机阵；rand返回在0，1区间上的一个随机数；将上面的rand写成randn则可以生成均值为 0、方差为1的正态分布的随机变量。 -1-6全1矩阵生成函数ones(m,n)生成mn阶全1矩阵7全0矩阵生成函数zeros(m,n)生成mn阶全0矩阵8离散序列绘图函数stem stem(y)以 1、 2、3为横坐标，y为纵坐标画杆形图；stem(x,y)以x为横坐标，y为纵坐标画杆形图(x与y数据个数必须一致)；stem(,fill)选项fill指定杆顶为实心，若无此选项则默认空心。 stem(.,LineSpec)参数LineSpec指定杆形图的线形、数据标志符号及颜色，具体用法可查看MATLAB帮助9线性坐标平面绘图函数plot用法与stem类似，具体用法可查看MATLAB帮助以上为MATLAB内置函数（在此仅为同学复习MATLAB提供） 四、在MATLAB程序中变量赋值注意问题在MATLAB中，对变量赋值时其维数可以按需要动态地改变，这样虽然方便程序设计但同时容易出错。 另外，频繁分配变量空间会大大降低程序的执行速度，因而应该尽量避免不必要的矩阵、向量维数的改变。 通常先用zeros()函数给变量分配足够大小的空间，再对变量进行赋值。 例依次执行下面的语句tic%开始计时for i=1:10000c(i)=i;%每次都重新分配空间end toc%读取计时时间tic%开始计时d=zeros(1,10000);%预先分配空间for i=1:10000d(i)=i;%直接赋值，不必重新分配空间end toc%读取计时时间运行结果如下elapsed_time=11560-2-elapsed_time=00470从结果可以看出，第2种赋值方法所用的时间比第1种方法所用时间少得多（以上是在主频为266GHZ的机器上运行的结果）。 -3-实验2序列与系统的时域表示、运算和卷积 一、实验目的1用矩阵乘法实现有限长序列卷积的计算；2序列的抽取(decimation,dilation)。 二、实验原理1有限长序列x(n)和h(n)的卷积也可以用矩阵乘实现。 把卷积和y(n)与输入序列x(n)用列向量y和x表示，则有y=Hx H称为Toeplitz矩阵，由单位样本响应h(n)的元素构成。 2序列的抽取(decimation,dilation)对序列的抽取操作定义为y(n)=x(nM)M是一整数。 三、实验内容与要求同学们以后在用MATLAB写程序时，应该加入适当的注释说明。 1用MATLAB解以下两题。 （1）已知序列x?n?1,2,3,4和h?n?3,2,1?a构造Toeplitz矩阵，它的维数与h(n)和x(n)的长度有关设h和x的长度分别为N h、N x，那么y的长度是多少？H的维数呢？b利用给出的序列验证你的函数。 （2）在已知第一行和第一列下，MATLAB提供一个称之为toeplitz矩阵的函数用于产生Toeplitz矩阵。 a观察矩阵H的第一列和第一行有什么特点，利用toeplitz函数，建立另外一个MATLAB函数用于实现线性卷积。 写成x、h是行或列向量都可以调用。 这个函数的格式应该是Function?y,H?=conv_tp(h,x)%Linear Convolutionusing ToeplitzMatrix%-%?y,H?=conv_tp(h,x)%y=output sequence-4-%H=Toeplitz matrixcorresponding tosequence hso thaty=Hx%h=Impulse responsesequence%x=input sequenceb对题 （1）给出的序列验证你的函数。 2若x?n?,?2,4,3,?6,5,?1,8,?那么减采样因子2的序列为y?n?,?2,3,5,8,?a建立一MATLAB函数dnsample，它有形式为function y=dnsample（x，M）用于实现上述运算。 细心关注在时间轴上的原点n=0，利用MATLAB的编序机理。 b产生x?n?sin?0.125?n?,?50?n?50,以4对x?n?抽取得到y?n?。 利用subplot画出x?n?和y?n?，并对结果作评述。 c用x?n?sin?0.5?n?,?50?n?50,重做（b）。 定性讨论信号在减采样后的效果。 3选作序列的内插运算(interpolation,or up-sampling)y(n)=x(n/L),n=kL,k为整数=0，n不满足以上条件时用MATLAB函数表示，应该是y=upsampling(x,L)注意upsample是MATLAB信号处理工具箱中的函数，所以我们用upsampling。 MATLAB还提供了一个对信号作内插、抽取、滤波的函数y=resample(x,p,q)对信号做p/q（p,q均为整数）的速率调整。 我们的实验只是抽取、内插。 四、参考文件及程序2个音频波形文件供有兴趣的同学们做抽取、内插实验用，其中pronounciation.wav是C+之父Bjarne Stroustrup告诉我们他自己的名字应该怎么念。 注意这个文件是立体声的（双声道）。 MATLAB Notebook读音频文件.doc，告诉你如何读音频文件，如何播放音频序列。 附音频文件的读入与音频信号播放 1、音频文件的读入x=wavread(orig.wav);x,Fs,bits=wavread(orig.wav);.=wavread(orig.wav,N)；-5-.=wavread(orig.wav,N1N2)siz=wavread(orig.wav,size)第1句从当前子目录下的音频波形文件orig.wav把数字音频数据读入矩阵x；如果波形文件不在当前目录，应该在波形文件名参数加入相应的路径。 当读入单声道数据时，x是一列向量，读入立体声时，x为2列的矩阵，分别对应2个声道。 数据x每个样点值的范围在-1，1之间。 第2句除了读入波形数据，还把数字音频信号的抽样频率读入变量Fs，每个样本占用的比特数放在变量bits。 如果不播放声音或把音频数据写入文件，不用改动后两个参数；你应该看出来，第1句其实是第2句的简化，只读入了音频数据，没有读其它参数；第 3、第4句写法表示赋值等号左边写法同上2句，但第3句只读入每个声道的前N点数据，而第4句读出每个声道从N1到N2点的数据。 第5句读入音频文件的相应信息，size是这种用法的固定参数，这时将不读入音频数据。 返回siz是一个2元数组，组成成分如下siz=samples channels第一个元素siz (1)是音频文件中所含的样本数，也就是第一句读入的x的长度，第二个元素siz (2)是该文件所含数字音频的声道数，单声道是1，立体声为2。 求取向量x的长度也可以用Len=length(x)这样就不用用上面的第2句。 但函数length()总是返回矩阵最大的维数，所以使用它不能获得信号的信道数信息（mono orstereo）。 2、音频信号播放sound(y,Fs)sound(y,Fs,bits)sound(y)y是音频信号向量/矩阵Fs是抽样频率,缺省值是8192Hz bits数字音频每个样本的比特数，8或16一般播放从音频文件读出的数据用第一种形式就可以，Fs是数字音频自己的抽样频率。 3、音频信号播放实例如果你有耳机可以插在计算机面板上那个浅黄色的插孔。 -6-如果要直接执行下面的MATLAB语句，要先设置好MATLAB的notebook，具体步骤参见为实验1提供的文件。 下面的语句打开并播放orig.wav（这个文件所在目录应该是当前目录，否则应在文件名字符串中加入文件所在路径）x,Fs,bits=wavread(orig.wav);sound(x,Fs)再听这句呢？sound(x)是不是不一样？这是因为播放的缺省速率和原来信号的速率不一样。 如果只播放信号向量的一部分，对单声道一维向量，播放N1:N2区间的样本，可以用sound(x(N1:N2),Fs)对二维立体声，则用sound(x(N1:N2,:),Fs)下面这句返回pronounciation.wav数据大小的参数siz=wavread(pronounciation.wav,size)siz=2658472可以看到，pronounciation.wav有2个声道，如果读入数据x,Fs,bits=wavread(pronounciation.wav);这时x应该视为265647x2的矩阵，即2列。 注意在一次播放音频信号未结束时，不要做另一次播放。 -7-实验3离散信号、系统的频域表示 一、实验目的1.考察抽样间隔对信号频谱的影响；2.考察信号观察时间对信号频谱的影响。 二、实验原理根据对连续时间信号抽样与重建的原理，考察抽样间隔对信号频谱的影响。 三、实验内容与要求在本次实验中，同学们要阅读、运行一些MATLAB程序，观察显示结果，然后对一些实验参数作调整后进一步观察，而后得出自己的结论。 1.Lab3_1.m是对正弦信号的抽样改变抽样间隔（哪个参数？），取若干不同的值作观察。 这是为了下面的实验做准备，可以不反映在实验报告中；2.Lab3_2.m是在时域观察抽样过程引起的混迭现象通过对抽样信号插值来恢复原信号。 绿色是内插恢复的信号，触一个键后显示的蓝色是原信号。 试着改变抽样间隔来改变恢复的效果。 3.Lab3_3.m在频域显示抽样引起的混迭现象本程序中的模拟信号的解析表达式是什么？观察改变抽样间隔引起的频域混迭现象，能得出什么结论？ 四、参考文件及程序清单Lab3_1.m，Lab3_2.m，Lab3_3.m%Program Lab3_1.m%Illustration ofthe SamplingProcess in the Time-Domain clf;t=0:0.0005:1;f=13;xa=cos(2*pi*f*t);subplot(2,1,1)plot(t,xa);grid-8-xlabe l(Time,msec);ylabel(Amplitude);title(Continuous-time signalx_a(t);axis(01-1.21.2)subplot(2,1,2);%in thelab project,should runfor4different samplingperiods T=0.1;n=0:T:1;xs=cos(2*pi*f*n);k=0:length(n)-1;stem(k,xs);grid;xlabel(Time indexn);ylabel(Amplitude);title(Discrete-time signalxn);axis(0(length(n)-1)-1.21.2)%Program Lab3_2.m%Illustration ofAliasing Effectin theTime-Domain%Program adaptedfromKra94with permissionfrom%The Mathworks,Inc.,Natick,MA.clf;T=0.1;f=13;n=(0:T:1);xs=cos(2*pi*f*n);t=linspace(-0.5,1.5,500);%下句参考p.67(3-33)和p.69(3-35)ya=sinc(1/T)*t(:,ones(size(n)-(1/T)*n(:,ones(size(t)*xs;plot(n,xs,o,t,ya);grid;xlabel(Time,msec);ylabel(Amplitude);title(Reconstructed continuous-time signaly_a(t);axis(01-1.21.2);%plot theoriginal analogsignal pause%waiting userresponse-9-xa=cos(2*pi*f*t);hold onplot(t,xa)hold off%Program Lab3_3%Illustration ofthe AliasingEffect%intheFrequency-Domain clf;t=0:0.005:10;xa=2*t.*exp(-t);subplot(2,2,1)plot(t,xa);grid xlabel(Time,msec);ylabel(Amplitude);title(Continuous-time signalx_a(t);subplot(2,2,2)wa=0:10/511:10;ha=freqs(2,121,wa);plot(wa/(2*pi),abs(ha);grid;xlabel(Frequency,kHz);ylabel(Amplitude);title(|X_a(jOmega)|);axis(05/pi02);subplot(2,2,3)T=1;n=0:T:10;xs=2*n.*exp(-n);k=0:length(n)-1;stem(k,xs);grid;xlabel(Time indexn);ylabel(Amplitude);title(Discrete-time signalxn);subplot(2,2,4)-10-wd=0:pi/255:pi;hd=freqz(xs,1,wd);plot(wd/(T*pi),T*abs(hd);grid;xlabel(Frequency,kHz);ylabel(Amplitude);title(|X(ejomega)|);axis(0min(1/T,5/pi)02) 五、实验中所用部分函数介绍1线性等分向量生成函数linspace linspace函数生成线性等分向量，其功能类似冒号算子“”，但它不象冒号那样，根据给定的起始值、增量和终止值控制生成向量元素的个数，而是直接给出元素个数，从而给出各个元素的值。 linspace函数的格式为 （1）x=linspace(a,b)生成有100个元素的行向量x，它的元素值在a和b之间线性分布。 （2）x=linspace(a,b,n)生成有n个元素的行向量x，它的元素值在a和b之间线性分布。 -11-实验4DFT&FFT 一、实验目的掌握用FFT做谱分析的方法。 二、实验内容与要求 1、用DFT/FFT对模拟信号做傅里叶分析参考上次的实验，以频率f s对以下信号抽样N点xa(t)=cos(a t)+cos(b t)+cos(c t)相应的参数是a=2*pi*6500,b=2*pi*7000,c=2*pi*9000f s=32000，N=16对这N点序列作N点DFT，观察其幅频特性，如果X=fft(x)w是频率坐标向量，你可以考虑用stem(w,abs(X),plot(w,abs(X),plot(w,abs(X),*)来显示，然后确定用哪种显示方式。 接下来，对x作M=256点FFT，这意味着在x后补了M-N个0，再观察幅频特性；现在令抽样点数N=256，再对这个抽样序列作N点FFT，观察幅频特性。 通过这些观察，你能作出什么判断？试着改变f s，令其分别为24000，19000，18000，17000，16000，你看到了什么，做怎样的解释？注意安排你的时域、频域的显示。 2、参考以下程序段，对DFT和FFT算法做计算效率比较Nmax=2048;fft_time=zeros(1,Nmax);for n=1:Nmax x=rand(1,n);t=clock;fft(x);fft_time(n)=etime(clock,t);end n=1:Nmax;plot(n,fft_time,.)xlabel(N);ylabel(Time inSec.)title(FFT executiontimes)比如，从N1点到N2点内的效率比较，设N1=8，N2可以考虑取256，或更大些。 -12-大致需要这样一个循环for n=N1:N2产生数据（可以用固定数据，或是随机数据）计时开始DFT计时结束，统计计时开始FFT计时结束，统计end显示统计数据计时函数clock,tic等参阅联机帮助。 3、对抽取/内插前后的信号做傅里叶分析本次实验的信号均假定起始时间下标为0，也就是对信号x(n)作N:1抽取，只要y=x(1:N:end)即可，而1:N内插则为y=zeros(1,N*length(x);y(1:N:length(y)=x;我们要着重观察的是抽取、内插后频谱的改变。 本实验的数据放在updn.mat文件内，执行语句load updn要用的数据就会载入数组siga和sigb，如何获取数组大小的信息？对siga做抽取，sigb做内插实验，试用N=2,3,4做抽取，内插，观察它们频谱的变化。 提示做谱分析时补一定的0在序列尾部。 -13-实验5利用FFT计算线性卷积 一、实验目的掌握用FFT计算线性卷积的方法。 一、实验原理对于输入信号序列x(n)，一个单位样本响应为h(n)的LSI系统的输出y(n)为两序列的卷积和y(n)=x(n)*h(n)如h(n)长度有限，则可以考虑把输入x(n)分成有限长子序列之和，然后用循环卷积（通过FFT）计算系统响应y(n)。 教科书中已经叙述了重叠-保留(overlap-save)算法，本实验要求同学们实现重叠-相加(overlap-add)算法，要用FFT实现循环卷积。 二、实验内容与要求参考重叠-保留算法的MATLAB实现ovrlpsav(x,h,N)和hsolpsav(x,h,N)（在E:DSP-ARefProgramPWSK_DSP目录下），用MATLAB实现下题。 1、块卷积的重叠相加法是重叠保留法的另一种替代方法。 设x?n?是一个很长的序列，长度为ML，其M,L?1。 现将x?n?分割为M段x m?n?,m?1,?,M，每段长为L。 M?1?x?n?,mL?n?m?1?L?1x?n?使有x?n?x m?n?m0,其余nm?0?设h?n?是L点的脉冲响应，那么y?n?x?n?h?n?x m?n?h?n?y m?n?;y m?n?x m?n?h?n?m?0m?0M?1M?1显然，y m?n?是（2L-1）点序列。 在这个方法中必须要保留中间卷积结果，然后在相加之前恰当地重叠这些结果以形成最后结果y?n?。 为了对这个运算应用DFT，必须要选N?2L?1?。 a利用循环卷积运算建立一个MATLAB函数实现重叠相加法，这个函数的格式应该是function?y?=ovrlpadd(x,h,N)%Overlap-Add methodof blockconvolution-14-%?y?=ovrlpadd(x,h,N)%y=output sequence%x=input sequence%h=impulse response%N=block length=2*length(h)-1b在上面函数中结合基-2FFT实现求得一个高速重叠相加块卷积程序。 记住选N?2。 c对下面两个序列验证你的函数vx?n?cos?n/500?4000?n?,h?n?1,?1,1,?1?函数格式如下functiony=hsolpadd(x,h,N)%High-speed Overlap-Add methodof blockconvolutions usingFFT%-%y=hsolpadd(x,h,N)%y=output sequence%x=input sequence%h=impulse response%N=block length(must bea powerof two)%重点是要实现用FFT的算法，即相应于hsolpsav的函数hsolpadd，要求用2的整数次幂的点数做FFT。 可以先写在时域实现的ovrlpadd，再写用FFT实现的hsolpadd，但不限定一定要先写ovrlpadd。 写好后，还应该写一个测试程序，检验一下，就是题中的c部分。 如何检验？可以用函数conv或书中重叠-保留算法的hsolpsav函数来做同样序列的卷积，比较它们的计算结果。 但要注意，因为各自的算法不同，运算产生的结果类型（实型、复型）和计算误差可能不同，不能简单地比较是否相等。 也可以绘图比较，绘图比较时，但是要注意，我们得到结果序列较长（4千多点），在一屏全部显示是我们用的显示器显示精度达不到的。 所以绘图显示是辅助的，主要靠做数值分析。 为方便检验，还可以用别的信号、响应序列。 三、思考题1.设x(n)的长度为M，h(n)的长度为N，问完成两序列的线性卷积需要多少次乘、加法运算？-15-2.设L=2k，则L点FFT需要L/2log2(L)次乘法。 试比较一下直接卷积和循环卷积（借助FFT）算法的乘法计算量。 如果保持L=M+N-1，即恰好等于线性卷积的长度，且是2的整数次幂，讨论一下M、N（必然是一增一减）变化对计算效率的影响。 四、指数和对数函数表下表列出了本次实验中的可能用到的指数和对数函数，这些函数是针对数组，按数组运算规则进行运算的，取数组内每个元素的函数值。 函数名log2log10功能以2为底的对数常用对数函数名log exp功能自然对数自然指数函数名pow2sqrt功能以2为底的指数平方根-16-实验6FIR数字滤波器的设计 一、实验目的掌握用MATLAB设计FIR数字滤波器的方法，重点要掌握窗函数和最优化设计方法。 二、实验原理用MATLAB提供的函数，设计FIR数字滤波器 三、实验内容与要求1.参考示范程序窗法Examples7.8-11（在E:DSP-ARefProgramCHAP_07目录下）最优化设计（Parks-McClellan-Remez）Examples7.23-25（在E:DSP-ARefProgramCHAP_07目录下）在阅读这些例题时，应该注意A)如何从滤波器指标要求导出其它设计参数，如确定窗口类型、滤波器的阶数等；B)例题中验证设计得到的滤波器是否满足设计指标的语句；C)特别是在最优化设计的例中的迭代设计过程；D)把滤波器的设计和其特性的显示分开，你自己的显示不一定要照搬书上的。 特别最优化设计的那几个例子程序，显示部分需要函数ampl_res()，需要自己写。 在这基础上，试着直接用我们讲的MATLAB函数fir1()、fir2()进行设计2.用Kaiser窗设计一个FIR数字带阻滤波器，对模拟信号xa(t)=cos(a t)+cos(b t)+cos(c t)a=2*pi*6500，b=2*pi*7000,c=2*pi*9000滤波，要求滤去7000Hz的频率成分。 系统采样率为fs=32000Hz这同我们第四次实验。 但采样点数应该比较大，可以用N=4096。 滤波器的Rp=0.25dB，As=50dB，过渡带宽可以用模拟频率（例如200Hz）也可以用数字频率指定。 还可以改变As（比如30dB）观察滤波效果。 这个实验一般不须在时域观察，主要应在频域作观察。 用freqz(x,1)-17-就能观察有限长序列的DTFT，对对数显示不习惯的，可以用H,w=freqz(x,1)；得到数据，再用freqzplot(H,w,linear)等方式显示。 具体用法参见参考程序。 四、参考程序Work35DSPLAB 61、FIR1Demo1.m用fir1函数设计高通滤波器例，各种窗口都试了一下； 2、FIR1Demo2.m用fir1函数设计多通带滤波器例； 3、FIR2Demo3.m用fir2函数设计多通带滤波器例； 4、FIRLsRemesDemo.m用firls和remez函数设计例，fir2也在这里做了对照。 1和4项所列的文件注释详细些。 FIR1Demo1.m%用窗函数设计线性相位FIR滤波器%教科书讲到的6种窗都用一下%本例未涉及如何确定滤波器的设计参数figNo=1;%显示用窗口号，接连用了3个n=48%滤波器阶数Wn=0.35;%截止频率beta=0.1102*(60-8.7);%Kaiser窗参数，假设阻带要求60分贝衰减(p.253)%生成窗口矩阵，各窗函数看教科书pp.243-53%先创建一个空矩阵，然后再把各窗函数列向量加进去Windows=;Windows=Windows,rectwin(n+1);Windows=Windows,bartlett(n+1);Windows=Windows,blackman(n+1);Windows=Windows,hann(n+1);Windows=Windows,hamming(n+1);Windows=Windows,kaiser(n+1,beta);%接下来用一个循环完成用各窗口函数设计的滤波器-18-%并得到相应的传输函数向量，按列放在矩阵Hs中%在这个循环中，win依次取矩阵Windows的各列Hs=;for win=Windows b=fir1(n,Wn,high,win);%把high参数去掉就是低通滤波器h,w=freqz(b,1);Hs=Hs,h;end%现在用3种尺度来显示,请观察对比各窗口figure(figNo)freqzplot(Hs,w,linear)figure(figNo+1)freqzplot(Hs,w,squared)figure(figNo+2)freqzplot(Hs,w)%用分贝数显示figure(figNo)%把第一个窗口推到前头FIR1Demo2.m%用窗函数设计线性相位FIR带通/带阻滤波器n=48;Wn=0.350.65;%通带或阻带%Haming windowby defaultb=fir1(n,Wn,stop);%再省去第三个参数stop就是带通win=rectwin(n+1);%矩形窗，宽度是滤波器阶数+1bb=fir1(n,Wn,stop,win);H,w=freqz(b,1,512);HH=freqz(bb,1,512);TH=H,HH;freqzplot(TH,w,squared)-19-FIR2Demo3.m%用窗函数设计线性相位FIR滤波器%设计2个带通/带阻滤波器演示n=48f=0.20.40.60.8;%f的元素一定不能包含0和1！bhm=fir1(n,f,DC-1);%DC-0表示频率0,0.2是阻带，DC-1是通带win=rectwin(n+1);%上面用缺省海明窗，下面用的矩形窗brec=fir1(n,f,DC-1,win);%显示Hm,w=freqz(bhm,1);Hrec=freqz(brec,1);H2=Hm,Hrec;freqzplot(H2,w,linear)FIRLsRemesDemo.m%最优化FIR线性相位滤波器设计例%均方误差最小准则%b=firls(n,f,m)%最大误差最小化准则，雷米兹交换算法%b=remez(n,f,m)%参数%b:返回的滤波器传输函数分子多项式的系数，也就是单位样本响应%n:滤波器的阶数%f:行向量，以0开始递增，各关键频率点，以1结束%m:行向量，维数与f相同，各关键频率点对应的响应幅度值%n=20;%Filter order%f=00.40.51;%Frequency bandedges%m=1100;%Desired amplitudes-20-%下面的例子是用这两种方法设计的两个通带的滤波器%再加上FIR2设计的作为对照%useplot是控制显示方式的变量，参看下面显示部分useplot=0;%试着改动各参数观察一下，比如改变过渡带的宽度、阶数等n=40;m=11001100;f=00.30.40.50.60.70.81;b1=firls(n,f,m);b2=remez(n,f,m);%用FIR2设计b3=fir2(n,f,m);%下面的代码显示比较它们的频率响应Hb1,w=freqz(b1,1);%分母多项式只有常数项1，没指定返回点数，缺省512点Hb2=freqz(b2);%这里把分母项省略了，返回的点数和上行一样，所以不再需要w Hb3=freqz(b3);if useplot%用一句plot画3条曲线，第2条用绿色，虚线，第3条红色；grid画出网格线subplot(2,1,1)plot(w/pi,abs(Hb1),w/pi,abs(Hb2),g-,w/pi,abs(Hb3),r),grid subplot(2,1,2)%plot(w/pi,angle(Hb),w/pi,angle(Hbb),r-),grid plot(w/pi,unwrap(angle(Hb1),w/pi,unwrap(angle(Hb2),g-,w/pi,unwrap(angle(Hb3),r)grid else%Hb,Hbb是和w同维数的列向量(从Wokspace里可以看到)%H是2列合成的矩阵H=Hb1,Hb2,Hb3;-21-%也可以用下面的来显示%用不同的方式显示，如果三种都要同时看，你还得加几句%例如用pause在时间上隔开，或是用3个窗口figure (1)freqzplot(H,w)figure (2)freqzplot(H,w,linear)figure (3)freqzplot(H,w,squared)end 六、实验所用数值处理函数表MATLAB中提供了一些数值处理函数，在编程时常会用到，它们也是按照数组运算规则进行运算的。 函数名fix功能向零取整函数名ceil功能向正无穷方向取整函数名mod功能除法求余（与除数同号）除法求余（与被除数同号）round四舍五入floor向负无穷方向取整rem-22-实验7IIR数字滤波器的设计 一、实验目的掌握用MATLAB设计IIR数字滤波器的方法。 二、实验原理用MATLAB提供的函数，设计IIR数字滤波器。 传统的4种类型IIR数字滤波器的设计，可以概括为这样4组MATLAB函数Butterworth滤波器b,a=butter(n,Wn)b,a=butter(n,Wn,ftype)输出参数b,an+1阶行向量。 相应滤波器的系统函数H(z)为B(z)b (1)?b (2)z?1?b(n?1)z?n H(z)?1?nA(z)1?a (2)z?a(n?1)z输入参数n:滤波器阶数Wn:频率参数标量w3分贝频率点(0 缺省时依频率参数Wn的形式分别表示低通或带通high:high passstop:band stopButterworth滤波器的阶数n和频率参数Wnn,Wn=buttord(Wp,Ws,Rp,Rs)输入参数Rp:通带允许的最大波纹(db，用分贝表示)Rs:阻带衰减(db)Wp,Ws:通带、阻带的频率参数(01)低通(WpWs)带通(Ws (1) (1) (2) (2)-23-带阻(Wp (1) (1) (2) (2)通俗地说带通是“阻带包住通带”，带阻是“通带包住阻带”。 类似地，Chebyshev I和II型滤波器设计用b,a=cheby1(n,Rp,Wn)b,a=cheby1(n,Rp,Wn,ftype)n,Wn=cheb1ord(Wp,Ws,Rp,Rs)b,a=cheby2(n,Rs,Wn)b,a=cheby2(n,Rs,Wn,ftype)n,Wn=cheb2ord(Wp,Ws,Rp,Rs)注意cheby1（）用Rp而cheby2（）用Rs。 而椭圆滤波器设计用b,a=ellip(n,Rp,Rs,Wn)b,a=ellip(n,Rp,Rs,Wn,ftype)n,Wn=ellipord(Wp,Ws,Rp,Rs)这些函数的参数与Butterworth滤波器设计函数相同。 三、实验内容与要求1.设计一个IIR带阻滤波器，可以用任何类型的，由于设计方式相似，改用另一类型引起的程序改动不大，所以鼓励同学们设计多于一种类型的滤波器。 仿照上次实验的场景，滤去模拟信号7000Hz的成分xa(t)=cos(a t)+cos(b t)+cos(c t)a=2*pi*6500,b=2*pi*7000,c=2*pi*9000系统采样率为fs=32000Hz滤波器的通带波纹Rp=0.25dB，阻带衰减As=50dB，过渡带宽可以用模拟频率（例如200Hz）也可以用数字频率指定。 这里选的设计参数与上次FIR滤波器设计一样，是为了让同学们对FIR和IIR滤波器的特性作比较。 注意IIR滤波器单位样本响应、滤波与零、极点的观察。 单位样本响应（unit sampleresponse）IIR滤波器的单位样本响应可以用MATLAB函数impz（）得到h,t=impz(b,a)；-24-h,t=impz(b,a,n)；h单位样本响应，列向量；th对应的时间支集，t=0:n-1。 在第一种形式，响应长度点数由impz自行确定，用户可以用第二种形式指定返回的响应长度。 如果不返回结果，impz将直接在当前图形窗口显示单位样本响应，如impz(b,a)对所得滤波器的频率响应也不应是freqz(h,1)，应为从理论上说h是无限长的。 滤波IIR