2011信号处理A指导书正文.doc

实验1 矩阵、序列的运算与显示一、实验目的掌握用MATLAB表示离散时间信号，以及它们的运算和显示。二、实验内容与要求1 用MATLAB产生并画出下列序列的样本。abcd，式中是在【-1，1】之间均匀分布的随机序列，你如何表征这个序列？e为周期的，画出5个周期。 三、实验所用部分函数如下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、cosy = sin(x) ,y = cos(x) (注意:x以弧度为单位)5. 随机序列生成函数rand，用法如：Y = rand (n) 生成nn阶的均匀分布随机阵；Y = rand (m, n) 生成mn阶的随机阵；rand 返回在0，1区间上的一个随机数；将上面的rand写成randn则可以生成均值为0、方差为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; %每次都重新分配空间endtoc %读取计时时间tic %开始计时d=zeros(1,10000); %预先分配空间for i=1:10000d (i) = i; %直接赋值，不必重新分配空间endtoc %读取计时时间运行结果如下：elapsed_time = 11560elapsed_time =00470从结果可以看出，第2种赋值方法所用的时间比第1种方法所用时间少得多（以上是在主频为266GHZ的机器上运行的结果）。实验2 序列与系统的时域表示、运算和卷积一、实验目的1 用矩阵乘法实现有限长序列卷积的计算；2 序列的抽取(decimation, dilation)。二、实验原理1 有限长序列x(n)和h(n)的卷积也可以用矩阵乘实现。把卷积和y(n)与输入序列x(n)用列向量y和x表示，则有y = HxH称为Toeplitz矩阵，由单位样本响应h(n)的元素构成。2 序列的抽取(decimation, dilation)对序列的抽取操作定义为y(n) = x(nM)M是一整数。三、实验内容与要求同学们以后在用MATLAB写程序时，应该加入适当的注释说明。1 用MATLAB解以下两题。（1）已知序列a 构造Toeplitz矩阵，它的维数与h(n)和x(n)的长度有关：设h和x的长度分别为Nh、Nx，那么y的长度是多少？H的维数呢？b 利用给出的序列验证你的函数。（2）在已知第一行和第一列下，MATLAB提供一个称之为toeplitz矩阵的函数用于产生Toeplitz矩阵。a 观察矩阵H的第一列和第一行有什么特点，利用toeplitz函数，建立另外一个MATLAB函数用于实现线性卷积。写成x、h是行或列向量都可以调用。这个函数的格式应该是Function =conv_tp(h,x)%Linear Convolution using Toeplitz Matrix%-%=conv_tp(h,x)%y= output sequence %H=Toeplitz matrix corresponding to sequence h so that y = Hx%h =Impulse response sequence %x = input sequence b对题（1）给出的序列验证你的函数。2 若那么减采样因子2的序列为a 建立 一MATLAB函数dnsample，它有形式为function y =dnsample（x，M）用于实现上述运算。细心关注在时间轴上的原点n=0，利用MATLAB的编序机理。b 产生以4对抽取得到。利用subplot画出和，并对结果作评述。c 用重做（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)；. = wavread (orig.wav, N1 N2)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 or stereo）。2、音频信号播放sound (y, Fs)sound (y, Fs, bits)sound (y)y 是音频信号向量/矩阵Fs 是抽样频率, 缺省值是8192Hzbits 数字音频每个样本的比特数，8或16一般播放从音频文件读出的数据用第一种形式就可以，Fs是数字音频自己的抽样频率。3、音频信号播放实例如果你有耳机可以插在计算机面板上那个浅黄色的插孔。如果要直接执行下面的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 = 265847 2 可以看到，pronounciation.wav 有2个声道，如果读入数据x, Fs, bits = wavread (pronounciation.wav);这时x应该视为265647x2 的矩阵，即2列。注意：在一次播放音频信号未结束时，不要做另一次播放。实验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 of the Sampling Process in the Time-Domainclf;t = 0:0.0005:1;f = 13;xa = cos (2*pi*f*t);subplot (2,1,1)plot (t,xa); gridxlabe l(Time, msec); ylabel (Amplitude);title (Continuous-time signal x_a(t);axis (0 1 -1.2 1.2)subplot (2, 1, 2);% in the lab project, should run for 4 different sampling periods T = 0.1;n = 0:T:1;xs = cos (2*pi*f*n);k = 0:length(n)-1;stem (k, xs); grid;xlabel (Time index n); ylabel(Amplitude);title (Discrete-time signal xn);axis (0 (length(n)-1) -1.2 1.2)% Program Lab3_2.m% Illustration of Aliasing Effect in the Time-Domain% Program adapted from Kra94 with permission from% 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 signal y_a(t);axis (0 1 -1.2 1.2);% plot the original analog signalpause% waiting user responsexa = cos (2*pi*f*t);hold onplot (t, xa)hold off% Program Lab3_3% Illustration of the Aliasing Effect % in the Frequency-Domainclf;t = 0:0.005:10;xa = 2 * t .* exp (-t);subplot (2, 2, 1)plot (t, xa); gridxlabel (Time, msec); ylabel (Amplitude);title (Continuous-time signal x_a(t);subplot (2, 2, 2)wa = 0 : 10/511 : 10;ha = freqs (2, 1 2 1, wa);plot (wa/(2*pi), abs (ha); grid;xlabel (Frequency, kHz); ylabel (Amplitude);title (|X_a(jOmega)|);axis (0 5/pi 0 2);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 index n);ylabel (Amplitude);title (Discrete-time signal xn);subplot (2, 2, 4)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 (0 min (1/T, 5/pi) 0 2)五、实验中所用部分函数介绍1线性等分向量生成函数linspacelinspace函数生成线性等分向量，其功能类似冒号算子“：”，但它不象冒号那样，根据给定的起始值、增量和终止值控制生成向量元素的个数，而是直接给出元素个数，从而给出各个元素的值。linspace函数的格式为：（1）x=linspace(a,b) 生成有100个元素的行向量x，它的元素值在a和b之间线性分布。（2）x=linspace(a,b,n)生成有n个元素的行向量x，它的元素值在a和b之间线性分布。实验4 DFT & FFT一、实验目的掌握用FFT做谱分析的方法。二、 实验内容与要求1、用DFT/FFT对模拟信号做傅里叶分析参考上次的实验，以频率fs 对以下信号抽样N点xa(t) = cos (a t) + cos (b t) + cos (c t)相应的参数是a = 2*pi*6500, b = 2*pi *7000, c = 2*pi*9000fs = 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，观察幅频特性。通过这些观察，你能作出什么判断？试着改变fs ，令其分别为 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);endn=1:Nmax;plot(n,fft_time,.)xlabel(N);ylabel(Time in Sec.)title(FFT execution times)比如，从N1点到N2点内的效率比较，设N1 = 8，N2可以考虑取256，或更大些。大致需要这样一个循环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在序列尾部。实验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、块卷积的重叠相加法是重叠保留法的另一种替代方法。设是一个很长的序列，长度为ML，其。现将分割为M段，每段长为L。设是L点的脉冲响应，那么显然，是（2L-1）点序列。在这个方法中必须要保留中间卷积结果，然后在相加之前恰当地重叠这些结果以形成最后结果。为了对这个运算应用DFT，必须要选。a 利用循环卷积运算建立一个MATLAB函数实现重叠相加法，这个函数的格式应该是function =ovrlpadd(x,h,N)% Overlap-Add method of block convolution% = ovrlpadd(x,h,N)% y = output sequence% x = input sequence % h = impulse response % N = block length =2*length(h)-1b 在上面函数中结合基-2FFT实现求得一个高速重叠相加块卷积程序。记住选。函数格式如下：function y = hsolpadd (x, h, N)% High-speed Overlap-Add method of block convolutions using FFT% -% y = hsolpadd (x, h, N)% y = output sequence% x = input sequence% h = impulse response% N = block length (must be a power of two)%c 对下面两个序列验证你的函数：重点是要实现用FFT的算法，即相应于hsolpsav的函数hsolpadd，要求用2的整数次幂的点数做FFT。可以先写在时域实现的ovrlpadd，再写用FFT实现的hsolpadd，但不限定一定要先写ovrlpadd。写好后，还应该写一个测试程序，检验一下，就是题中的c部分。如何检验？可以用函数conv或书中重叠-保留算法的hsolpsav函数来做同样序列的卷积，比较它们的计算结果。但要注意，因为各自的算法不同，运算产生的结果类型（实型、复型）和计算误差可能不同，不能简单地比较是否相等。也可以绘图比较，绘图比较时，但是要注意，我们得到结果序列较长（4千多点），在一屏全部显示是我们用的显示器显示精度达不到的。所以绘图显示是辅助的，主要靠做数值分析。为方便检验，还可以用别的信号、响应序列。三、 思考题1. 设x(n)的长度为M，h(n)的长度为N，问完成两序列的线性卷积需要多少次乘、加法运算？2. 设L = 2k，则L点FFT需要L/2log2(L)次乘法。试比较一下直接卷积和循环卷积（借助FFT）算法的乘法计算量。如果保持L = M + N -1，即恰好等于线性卷积的长度，且是2的整数次幂，讨论一下M、N（必然是一增一减）变化对计算效率的影响。四、 指数和对数函数表下表列出了本次实验中的可能用到的指数和对数函数，这些函数是针对数组，按数组运算规则进行运算的，取数组内每个元素的函数值。函数名功能函数名功能函数名功能log2以2为底的对数log自然对数pow2以2为底的指数log10常用对数exp自然指数sqrt平方根实验6 IIR 数字滤波器的设计一、实验目的掌握用MATLAB设计IIR数字滤波器的方法。二、实验原理用 MATLAB 提供的函数，设计IIR数字滤波器。传统的4种类型IIR数字滤波器的设计，可以概括为这样4组MATLAB函数：Butterworth 滤波器：b, a = butter (n, Wn)b, a = butter (n, Wn, ftype)输出参数：b, a：n+1阶行向量。相应滤波器的系统函数H(z)为输入参数：n: 滤波器阶数Wn: 频率参数标量 w：3 分贝频率点(0w1)数组 w1, w2: (w1w2) band pass/stopftype: 滤波器类型。缺省时依频率参数Wn的形式分别表示低通或带通high: high passstop: band stopButterworth 滤波器的阶数n和频率参数Wnn, Wn = buttord (Wp, Ws, Rp, Rs) 输入参数：Rp: 通带允许的最大波纹(db，用分贝表示)Rs: 阻带衰减 (db)Wp, Ws: 通带、阻带的频率参数(01)低通(Wp Ws)带通(Ws (1) Wp (1) Wp (2) Ws (2)带阻(Wp (1) Ws (1) Ws (2) Wp (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 = 32000 Hz滤波器的通带波纹 Rp = 0.25 dB， 阻带衰减As = 50 dB， 过渡带宽可以用模拟频率（例如200Hz）也可以用数字频率指定。这里选的设计参数与上次FIR滤波器设计一样，是为了让同学们对FIR和IIR滤波器的特性作比较。注意：IIR滤波器单位样本响应、滤波与零、极点的观察。单位样本响应（unit sample response）：IIR滤波器的单位样本响应可以用MATLAB函数impz（）得到h, t = impz (b, a)；h, t = impz (b, a, n)；h：单位样本响应，列向量；t：h 对应的时间支集，t = 0 : n-1。在第一种形式，响应长度点数由impz自行确定，用户可以用第二种形式指定返回的响应长度。如果不返回结果，impz将直接在当前图形窗口显示单位样本响应，如impz (b, a)对所得滤波器的频率响应也不应是freqz (h, 1)，应为从理论上说h是无限长的。滤波：IIR滤波器设计得到的是系统函数H(z) = b(z)/a(z)的分子、分母多项式系数数组，不像FIR滤波器设计后直接得到单位样本响应，所以滤波一般不是经过h, t = impz（b, a）；得到单位样本响应，再经过函数conv（）得到输出响应y = conv(h, x)；而是直接用函数filter（）求滤波结果y = filter (b, a, x)；这时输出y的长度等于输入信号x的长度。IIR滤波器的零、极点观察：用函数zplane (b, a)。2. 选作：通过内插、抽取与滤波改变信号的数据速率。在实验4，曾经对抽取、内插过程作过富里叶分析，下图是抽取率为2时的幅度谱：下图是抽取率为3时的幅度谱，混叠现象（aliasing）很明显。实践中往往结合内插来改变信号的数据率。下图是内插率为3的幅度谱实践中应该先做内插，经过低通滤波，然后才做抽取。试着实现这一过程。内插和抽取率可以作为参数。低通滤波器的截止频率应该如何确定？使用IIR还是FIR滤波器？可以先用实验4的数据文件updn.mat，用它做谱分析效果应该比较清楚，然后可以试着处理音频数据文件 orig.wav，pronounciation.wav, global.wav。pronounciation.wav是C+之父Bjarne Stroustrup说他的名字应该如何发音的，你听得清吗？试着改变信号的速率看能不能容易听些。如何读音频文件，如何在MATLAB中播放声音，参见第2次实验的文件。四、 参考文件1. cheby1example.m 切比雪夫I型设计例；2. YulewalkExample.m Yule-Walker算法例，相应于FIR滤波器的最优化设计；3. MATLAB Notebook 读音频文件.doc；4. orig.wav，pronounciation.wav, global.wav：音频波形文件。cheby1example.m% 双线性变换法设计数字切贝雪夫低通滤波器% 显示用窗口号数。改变它以方便保留先前打开的窗口（比如与其它方法的结果做对照）figno = 1;ShowFlag = 1; % 控制是否显示零极点、频率响应等，设为 0 则不显示% 归一化的通带、阻带数字频率点Wp = 0.2; Ws = 0.3;% 响应的通带、阻带衰减参数（分贝）Rp = 0.1; Rs = 25;% 滤波器的阶数、“关键”频率n, Wn = cheb1ord (Wp, Ws, Rp, Rs)% 现在得到滤波器的系统函数参数b, a = cheby1(n, Rp, Wn);if ShowFlag = 0 % 显示零极点 figure (figno) zplane (b, a) % 显示幅度和相位响应 figure (figno+1) h,w = freqz (b, a); freqzplot(h, w, linear) figure (figno+2) freqz (b, a) % in DB % 显示滤波器的单位冲激响应 figure (figno+3) impz (b,a) endYulewalkExample.m% yulewalk 滤波器设计例% f: 数组，各“关键”频率点% m：数组，维数和 f 相同。目标滤波器在上述频率点的“理想”幅度% n: 滤波器阶数n = 20;m = 0 0 1 1 0 0 1 1 0 0;f = 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 1;b,a = yulewalk (n, f, m);% 频率响应，在0到pi内抽256点，如果省略这个参数，将返回512点h,w = freqz (b, a, 256);% 用黑色显示设计目标，红色显示得到的滤波器幅度响应subplot (2, 1, 1)% m+eps 免得解释器抱怨对0取对数plot (f, 20*log10(m+eps), k*, w/pi, 20*log10(abs(h), r)subplot (2, 1, 2)plot (f, m, k, w/pi, abs(h), r) - 28 -实验7 FIR 数字滤波器的设计一、实验目的掌握用MATLAB设计FIR数字滤波器的方法，重点要掌握窗函数和最优化设计方法。二、实验原理用 MATLAB 提供的函数，设计FIR数字滤波器三、 实验内容与要求1. 参考示范程序窗法： Examples 7.8-11（在E: DSP-ARefProgramCHAP_07目录下）最优化设计（Parks-McClellan-Remez）： Examples 7.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 = 32000 Hz这同我们第四次实验。但采样点数应该比较大，可以用 N = 4096。滤波器的 Rp = 0.25 dB， As = 50 dB， 过渡带宽可以用模拟频率（例如200Hz）也可以用数字频率指定。还可以改变As（比如30dB）观察滤波效果。这个实验一般不须在时域观察，主要应在频域作观察。用freqz(x, 1)就能观察有限长序列的DTFT，对对数显示不习惯的，可以用H, w = freqz(x, 1)；得到数据，再用freqzplot(H, w, linear)等方式显示。具体用法参见参考程序。四、 参考程序：Work35DSPLAB61、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);% 接下来用一个循环完成用各窗口函数设计的滤波器% 并得到相应的传输函数向量，按列放在矩阵 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.35 0.65; % 通带或阻带% Haming window by 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)FIR2Demo3.m% 用窗函数设计线性相位FIR滤波器% 设计 2 个带通/带阻滤波器演示