1、- 1 -实验 1 矩阵、序列的运算与显示一、实验目的掌握用 MATLAB 表示离散时间信号，以及它们的运算和显示。二、实验内容与要求1 用 MATLAB 产生并画出下列序列的样本。a 10()21,025mxnnmnb 2256.410u unc 3()0.9cos(./3),02nx nd ，式中 是在【-1，1】之间均匀24181w w分布的随机序列，你如何表征这个序列？e 为周期的，画出 5 个周期。 5.,3,.xn三、实验所用部分函数如下1.单位冲激序列（信号）生成函数 impseqx,n = impseq(n0,n1,n2)2.阶跃序列（信号）生成函数 stepseqx,n =
2、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 的正态分布的随机变
3、量。- 2 -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用
4、法与 stem 类似，具体用法可查看 MATLAB 帮助以上为 MATLAB 内置函数（在此仅为同学复习 MATLAB 提供）四、在 MATLAB 程序中变量赋值注意问题在 MATLAB 中，对变量赋值时其维数可以按需要动态地改变，这样虽然方便程序设计但同时容易出错。另外，频繁分配变量空间会大大降低程序的执行速度，因而应该尽量避免不必要的矩阵、向量维数的改变。通常先用 zeros()函数给变量分配足够大小的空间，再对变量进行赋值。例：依次执行下面的语句tic %开始计时for i=1:10000c (i) = i; %每次都重新分配空间endtoc %读取计时时间tic %开始计时d=zero
5、s(1,10000); %预先分配空间for i=1:10000d (i) = i; %直接赋值，不必重新分配空间endtoc %读取计时时间运行结果如下：elapsed_time =- 3 -11560elapsed_time =00470从结果可以看出，第 2 种赋值方法所用的时间比第 1 种方法所用时间少得多（以上是在主频为 266GHZ 的机器上运行的结果） 。- 4 -实验 2 序列与系统的时域表示、运算和卷积一、实验目的1 用矩阵乘法实现有限长序列卷积的计算；2 序列的抽取(decimation, dilation)。二、实验原理1 有限长序列 x(n)和 h(n)的卷积也可以用矩
6、阵乘实现。把卷积和 y(n)与输入序列 x(n)用列向量 y 和 x 表示，则有y = HxH 称为 Toeplitz 矩阵，由单位样本响应 h(n)的元素构成。2 序列的抽取(decimation, dilation)对序列的抽取操作定义为y(n) = x(nM)M 是一整数。三、实验内容与要求同学们以后在用 MATLAB 写程序时，应该加入适当的注释说明。1 用 MATLAB 解以下两题。（1）已知序列 1,2343,21xnhn和a 构造 Toeplitz 矩阵，它的维数与 h(n)和 x(n)的长度有关：设 h 和 x 的长度分别为Nh、N x，那么 y 的长度是多少？H 的维数呢？b
7、 利用给出的序列验证你的函数。（2）在已知第一行和第一列下，MATLAB 提供一个称之为 toeplitz 矩阵的函数用于产生Toeplitz 矩阵。a 观察矩阵 H 的第一列和第一行有什么特点，利用 toeplitz 函数，建立另外一个MATLAB 函数用于实现线性卷积。写成 x、h 是行或列向量都可以调用。这个函数的格式应该是Function =conv_tp(h,x),y%Linear Convolution using Toeplitz Matrix%-% =conv_tp(h,x),yH%y= output sequence - 5 -%H=Toeplitz matrix corre
8、sponding to sequence h so that y = Hx%h =Impulse response sequence %x = input sequence b对题（1）给出的序列验证你的函数。2 若 ,243,651,8xn 那么减采样因子 2 的序列为,y a 建立 一 MATLAB 函数 dnsample，它有形式为function y =dnsample（x，M）用于实现上述运算。细心关注在时间轴上的原点 n=0，利用 MATLAB 的编序机理。b 产生 以 4 对 抽取得到 。利用sin0.125,50,nxnynsubplot 画出 和 ，并对结果作评述。xyc 用
9、 重做（b） 。定性讨论信号在减采样后的效果。si.,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 均为整数）的速率调整。我们的实验只是抽取、内插。四、参考文件及程序
10、2 个音频波形文件供有兴趣的同学们做抽取、内插实验用，其中 pronounciation.wav 是 C+之父 Bjarne Stroustrup 告诉我们他自己的名字应该怎么念。注意：这个文件是立体声的（双声道） 。MATLAB Notebook 读音频文件.doc ，告诉你如何读音频文件，如何播放音频序列。附：音频文件的读入与音频信号播放1、音频文件的读入x = wavread (orig.wav);x, Fs, bits = wavread (orig.wav);- 6 -. = wavread (orig.wav, N)；. = wavread (orig.wav, N1 N2)siz
11、 = wavread (orig.wav, size)第 1 句从当前子目录下的音频波形文件 orig.wav 把数字音频数据读入矩阵 x；如果波形文件不在当前目录，应该在波形文件名参数加入相应的路径。当读入单声道数据时，x 是一列向量，读入立体声时，x 为 2 列的矩阵，分别对应 2 个声道。数据 x 每个样点值的范围在-1， 1之间。第 2 句除了读入波形数据，还把数字音频信号的抽样频率读入变量 Fs，每个样本占用的比特数放在变量 bits。如果不播放声音或把音频数据写入文件，不用改动后两个参数；你应该看出来，第 1 句其实是第 2 句的简化，只读入了音频数据，没有读其它参数；第 3、第
12、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 ()总
13、是返回矩阵最大的维数，所以使用它不能获得信号的信道数信息（mono or stereo） 。2、音频信号播放sound (y, Fs)sound (y, Fs, bits)sound (y)y 是音频信号向量/矩阵Fs 是抽样频率, 缺省值是 8192Hzbits 数字音频每个样本的比特数，8 或 16一般播放从音频文件读出的数据用第一种形式就可以，Fs 是数字音频自己的抽样频率。3、音频信号播放实例- 7 -如果你有耳机可以插在计算机面板上那个浅黄色的插孔。如果要直接执行下面的 MATLAB 语句，要先设置好 MATLAB 的 notebook，具体步骤参见为实验 1提供的文件。下面的语句打
14、开并播放 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(pronouncia
15、tion.wav, size) siz =265847 2 可以看到，pronounciation.wav 有 2 个声道，如果读入数据x, Fs, bits = wavread (pronounciation.wav);这时 x 应该视为 265647x2 的矩阵，即 2 列。注意：在一次播放音频信号未结束时，不要做另一次播放。- 8 -实验 3 离散信号、系统的频域表示一、实验目的1. 考察抽样间隔对信号频谱的影响；2. 考察信号观察时间对信号频谱的影响。二、实验原理根据对连续时间信号抽样与重建的原理，考察抽样间隔对信号频谱的影响。三、实验内容与要求在本次实验中，同学们要阅读、运行一些 M
16、ATLAB 程序，观察显示结果，然后对一些实验参数作调整后进一步观察，而后得出自己的结论。1. Lab3_1.m 是对正弦信号的抽样改变抽样间隔（哪个参数？） ，取若干不同的值作观察。这是为了下面的实验做准备，可以不反映在实验报告中；2. Lab3_2.m 是在时域观察抽样过程引起的混迭现象通过对抽样信号插值来恢复原信号。绿色是内插恢复的信号，触一个键后显示的蓝色是原信号。试着改变抽样间隔来改变恢复的效果。3. Lab3_3.m 在频域显示抽样引起的混迭现象本程序中的模拟信号的解析表达式是什么？观察改变抽样间隔引起的频域混迭现象，能得出什么结论？四、参考文件及程序清单Lab3_1.m，Lab3
17、_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); grid- 9 -xlabe 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 pr
18、oject, 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
19、% 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;
20、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 response- 10 -xa = cos (2*pi*f*t);hold onplot (t, xa)hold off% Program Lab3_3% Illustration of the Aliasing Effect % in the Frequency-D
21、omainclf;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);t
22、itle (|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)- 11 -wd = 0:pi/255:pi;hd = freqz (xs, 1, wd);plot (wd/(T*pi), T * abs (hd);
23、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
24、=linspace(a,b,n)生成有 n 个元素的行向量 x，它的元素值在 a 和 b 之间线性分布。- 12 -实验 4 DFT fft_time=zeros(1,Nmax);for n=1:Nmaxx=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，或更大些。- 13 -大致需
25、要这样一个循环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 文件内，执行语句loa
26、d updn要用的数据就会载入数组 siga 和 sigb，如何获取数组大小的信息？对 siga 做抽取，sigb 做内插实验，试用N = 2, 3, 4做抽取，内插，观察它们频谱的变化。提示：做谱分析时补一定的 0 在序列尾部。- 14 -实验 5 利用 FFT 计算线性卷积一、实验目的掌握用 FFT 计算线性卷积的方法。一、 实验原理对于输入信号序列 x(n)，一个单位样本响应为 h(n)的 LSI 系统的输出 y(n)为两序列的卷积和y(n) = x (n) * h (n)如 h (n)长度有限，则可以考虑把输入 x (n)分成有限长子序列之和，然后用循环卷积（通过FFT）计算系统响应
27、y (n)。教科书中已经叙述了重叠-保留 (overlap-save)算法，本实验要求同学们实现重叠-相加(overlap-add)算法，要用 FFT 实现循环卷积。二、 实验内容与要求参考重叠-保留算法的 MATLAB 实现 ovrlpsav(x, h, N)和 hsolpsav(x, h, N)（在 E: DSP-ARefProgramPWSK_DSP 目录下） ，用 MATLAB 实现下题。1、块卷积的重叠相加法是重叠保留法的另一种替代方法。设 是一个很长的序列，长xn度为 ML，其 。现将 分割为 M 段 ，每段长为 L。,1ML?xn,1,mxnM 010,nxmL 使 有其 余设
28、是 L 点的脉冲响应，那么h1100;MMmmmynxnxhnynxhn?显然， 是（2L-1）点序列。在这个方法中必须要保留中间卷积结果，然后在相加之前恰m当地重叠这些结果以形成最后结果 。为了对这个运算应用 DFT，必须要选 。yn 21NLa 利用循环卷积运算建立一个 MATLAB 函数实现重叠相加法，这个函数的格式应该是- 15 -function =ovrlpadd(x,h,N)y% Overlap-Add method of block convolution% = ovrlpadd(x,h,N)% y = output sequence% x = input sequence %
29、 h = impulse response % N = block length =2*length(h)-1b 在上面函数中结合基-2FFT 实现求得一个高速重叠相加块卷积程序。记住选 。2vN函数格式如下：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
30、 length (must be a power of two)%c 对下面两个序列验证你的函数：40os/5,1,xnnh重点是要实现用 FFT 的算法，即相应于 hsolpsav 的函数 hsolpadd，要求用 2 的整数次幂的点数做 FFT。可以先写在时域实现的 ovrlpadd，再写用 FFT 实现的 hsolpadd，但不限定一定要先写 ovrlpadd。写好后，还应该写一个测试程序，检验一下，就是题中的 c 部分。如何检验？可以用函数 conv 或书中重叠-保留算法的 hsolpsav 函数来做同样序列的卷积，比较它们的计算结果。但要注意，因为各自的算法不同，运算产生的结果类型（
31、实型、复型）和计算误差可能不同，不能简单地比较是否相等。也可以绘图比较，绘图比较时，但是要注意，我们得到结果序列较长（4 千多点） ，在一屏全部显示是我们用的显示器显示精度达不到的。所以绘图显示是辅助的，主要靠做数值分析。为方便检验，还可以用别的信号、响应序列。三、 思考题- 16 -1. 设 x(n)的长度为 M，h(n)的长度为 N，问完成两序列的线性卷积需要多少次乘、加法运算？2. 设 L = 2k，则 L 点 FFT 需要 L/2log2(L)次乘法。试比较一下直接卷积和循环卷积（借助FFT）算法的乘法计算量。如果保持 L = M + N -1，即恰好等于线性卷积的长度，且是 2 的整
32、数次幂，讨论一下 M、N（必然是一增一减）变化对计算效率的影响。四、 指数和对数函数表下表列出了本次实验中的可能用到的指数和对数函数，这些函数是针对数组，按数组运算规则进行运算的，取数组内每个元素的函数值。函数名 功能 函数名 功能 函数名 功能log2 以 2 为底的对数 log 自然对数 pow2 以 2 为底的指数log10 常用对数 exp 自然指数 sqrt 平方根- 17 -实验 6 IIR 数字滤波器的设计一、实验目的掌握用 MATLAB 设计 IIR 数字滤波器的方法。二、实验原理用 MATLAB 提供的函数，设计 IIR 数字滤波器。传统的 4 种类型 IIR 数字滤波器的设
33、计，可以概括为这样 4 组 MATLAB 函数：Butterworth 滤波器：b, a = butter (n, Wn)b, a = butter (n, Wn, ftype)输出参数：b, a：n+1 阶行向量。相应滤波器的系统函数 H(z)为1()(2)()nBzbzbzHAaa输入参数：n: 滤波器阶数Wn: 频率参数标量 w：3 分贝频率点(0 Ws)带通(Ws (1) Wp (1) Wp (2) Ws (2)- 18 -带阻(Wp (1) Ws (1) Ws (2) Wp (2)通俗地说：带通是“阻带包住通带” ，带阻是“通带包住阻带” 。类似地，Chebyshev I 和 II
34、型滤波器设计用：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 =
35、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
36、 = 50 dB， 过渡带宽可以用模拟频率（例如 200Hz）也可以用数字频率指定。这里选的设计参数与上次 FIR 滤波器设计一样，是为了让同学们对 FIR 和 IIR 滤波器的特性作比较。注意：IIR 滤波器单位样本响应、滤波与零、极点的观察。单位样本响应（unit sample response）：IIR 滤波器的单位样本响应可以用 MATLAB 函数impz（）得到h, t = impz (b, a)；- 19 -h, t = impz (b, a, n)；h：单位样本响应，列向量；t：h 对应的时间支集，t = 0 : n-1。在第一种形式，响应长度点数由 impz 自行确定，用户可以
37、用第二种形式指定返回的响应长度。如果不返回结果，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
38、的长度等于输入信号 x 的长度。IIR 滤波器的零、极点观察：用函数 zplane (b, a)。2. 选作：通过内插、抽取与滤波改变信号的数据速率。在实验 4，曾经对抽取、内插过程作过富里叶分析，下图是抽取率为 2 时的幅度谱：- 20 -下图是抽取率为 3 时的幅度谱，混叠现象（aliasing）很明显。实践中往往结合内插来改变信号的数据率。下图是内插率为 3 的幅度谱- 21 -实践中应该先做内插，经过低通滤波，然后才做抽取。试着实现这一过程。内插和抽取率可以作为参数。低通滤波器的截止频率应该如何确定？使用 IIR 还是 FIR 滤波器？可以先用实验 4 的数据文件 updn.mat，用
39、它做谱分析效果应该比较清楚，然后可以试着处理音频数据文件 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 读音
40、频文件.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
41、, Rs)% 现在得到滤波器的系统函数参数b, a = cheby1(n, Rp, Wn);- 22 -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: 数组，各 “关键”
42、频率点% 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),
43、 k*, w/pi, 20*log10(abs(h), r)subplot (2, 1, 2)plot (f, m, k, w/pi, abs(h), r) - 23 -实验 7 FIR 数字滤波器的设计一、实验目的掌握用 MATLAB 设计 FIR 数字滤波器的方法，重点要掌握窗函数和最优化设计方法。二、实验原理用 MATLAB 提供的函数，设计 FIR 数字滤波器三、 实验内容与要求1. 参考示范程序窗法： Examples 7.8-11（在 E: DSP-ARefProgramCHAP_07 目录下）最优化设计（Parks-McClellan-Remez）： Examples 7.23-
44、25（在 E: DSP-ARefProgramCHAP_07 目录下）在阅读这些例题时，应该注意：A) 如何从滤波器指标要求导出其它设计参数，如确定窗口类型、滤波器的阶数等；B) 例题中验证设计得到的滤波器是否满足设计指标的语句；C) 特别是在最优化设计的例中的迭代设计过程；D) 把滤波器的设计和其特性的显示分开，你自己的显示不一定要照搬书上的。特别最优化设计的那几个例子程序，显示部分需要函数 ampl_res ( )，需要自己写。在这基础上，试着直接用我们讲的 MATLAB 函数 fir1()、fir2()进行设计2. 用 Kaiser 窗设计一个 FIR 数字带阻滤波器，对模拟信号xa(t
45、) = 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)- 24 -就能观察有限长序列的 DTF
46、T，对对数显示不习惯的，可以用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% 用窗函
47、数设计线性相位 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);% 接下来用一个循环完成用各窗口函数设计的滤波器- 25 -% 并得到相应的传输函数向量，按列放在矩阵 Hs 中% 在这个循环中，win 依次取矩阵 Wi
