DSP-B指导书正文

1、实验1 矩阵、序列及系统的时域表示与运算一、实验目的掌握用MATLAB表示离散时间信号和系统，以及它们的运算和显示。二、实验内容与要求 （在预习报告中写出生成如下序列的程序）请在实验报告中记录以下程序和结果1 用MATLAB产生并画出下列序列的样本。abcd，式中是在-1，1之间均匀分布的随机序列，你如何表征这个序列？在实验报告回答用matlab生成随机序列的函数e为周期的，画出5个周期。 2若线性时不变系统的单位样值响应为，输入序列，求系统的输出，并画出其波形图。(思考：你可以用几种方法来实现？) 在实验报告回答表示LTI系统的方法和求解输出响应的方法在实验报告回答：函数plot、stem、

2、subplot、xlabel、ylabel、title、hold on/off的作用三、实验所用部分函数如下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 = r

3、and (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指定杆顶为实心，若无此选项则默认空心。 st

4、em(.,LineSpec) 参数LineSpec指定杆形图的线形、数据标志符号及颜色，具体用法可查看MATLAB帮助9线性坐标平面绘图函数plot 用法与stem类似，具体用法可查看MATLAB帮助以上为MATLAB内置函数（在此仅为同学复习MATLAB提供） 四、在MATLAB程序中变量赋值注意问题在MATLAB中，对变量赋值时其维数可以按需要动态地改变，这样虽然方便程序设计但同时容易出错。另外，频繁分配变量空间会大大降低程序的执行速度，因而应该尽量避免不必要的矩阵、向量维数的改变。通常先用zeros()函数给变量分配足够大小的空间，再对变量进行赋值。例：依次执行下面的语句tic %开始计

5、时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. 掌握用FFT做谱分析的方法。二、实验

6、内容与要求 （在预习报告中写出对信号进行傅立叶分析的程序，并思考下面的问题）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), *)来显示，然后确定用哪种显示方式。

7、注意安排信号的时域、频域的显示。1）对N=16点的序列作16点DFT（FFT），观察其幅频特性；2）对N=16点的序列作M=256点DFT（FFT），这意味着在x后补了M-N个0，再观察幅频特性；3）对N=256点的序列作256点DFT（FFT），观察其幅频特性；4）改变fs ，令其分别为 24000，19000，18000，17000，16000，你看到了什么，做怎样的解释？1）在实验报告中回答：离散信号如何通过对连续信号抽样来获得，并按照以上要求写出对xa(t)抽样后的结果（即xn的表达式）；2）在实验报告中回答：如何对信号进行频谱分析，fft函数的用法？3）在实验报告中回答：对原时域信号

8、补零后，其幅频特性有何变化？4）在实验报告中回答：高密度谱和高分辨率谱之间有何区别？5）在实验报告中回答：改变时域信号的采样率fs对信号的幅频特性会产生什么影响？为什么？在实际中选取fs的依据是什么？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

9、 execution times)比如，从N1点到N2点内的效率比较，设N1 = 8，N2可以考虑取256，或更大些。大致需要这样一个循环for n = N1:N2产生数据（可以用固定数据，或是随机数据）计时开始DFT计时结束，统计计时开始FFT计时结束，统计end显示统计数据计时函数clock, tic等参阅联机帮助。在实验报告中回答：函数clock、etime、tic、toc的作用和使用方法。实验3 FIR数字滤波器的设计一、实验目的1）掌握MATLAB中设计FIR数字滤波器的常用函数，并编写简单程序。2）熟悉用FDATool进行FIR滤波器设计与分析的方法。二、实验原理用 MATLAB

10、提供的函数，设计FIR数字滤波器三、实验内容与要求1. 用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）观察滤波效果

11、。这个实验一般不须在时域观察，主要应在频域作观察。用freqz(x, 1)就能观察有限长序列的DTFT，对对数显示不习惯的，可以用 H, w = freqz(x, 1)；得到数据，再用 freqzplot(H, w, linear)注意：A) 如何从滤波器指标要求导出其它设计参数，如确定窗口类型、滤波器的阶数等；B) 如何验证设计得到的滤波器是否满足设计指标的语句；试着直接用我们讲的MATLAB函数fir1()、fir2()进行设计 2. 利用FDATool设计一个如上要求的FIR滤波器，对比使用不同的方法设计的滤波器。四、参考程序：1、FIR1Demo1.m 用 fir1 函数设计高通滤波器

12、例，各种窗口都试了一下；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); % K

13、aiser 窗参数，假设阻带要求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);% 接下来用一个循环完成

14、用各窗口函数设计的滤波器% 并得到相应的传输函数向量，按列放在矩阵 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)freqzp

15、lot (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 =

16、 H, HH;freqzplot (TH, w, squared)FIR2Demo3.m% 用窗函数设计线性相位FIR滤波器% 设计 2 个带通/带阻滤波器演示n = 48f = 0 .2 0.4 0.6 0.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);

17、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 = 0 0.4 0.5 1; %

18、Frequency band edges%m = 1 1 0 0; % Desired amplitudes% 下面的例子是用这两种方法设计的两个通带的滤波器% 再加上 FIR2 设计的作为对照% useplot 是控制显示方式的变量，参看下面显示部分useplot = 0;% 试着改动各参数观察一下，比如改变过渡带的宽度、阶数等n = 40;m = 1 1 0 0 1 1 0 0;f = 0 0.3 0.4 0.5 0.6 0.7 0.8 1;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); % 这里把分母项省略了，返回的点数和上行一样，所以不再需要 wHb3 = freqz (b3);if useplot % 用一句 plot 画3条曲线，第2条用绿色，虚线，第3条红色；grid 画出网格线 subplot (2, 1, 1) p