数字信号处理实验指导书（带源程序）

实验一 离散时间系统与MATLAB

一. 实验目的

1. 进一步加深对离散时间系统的理解。

2. 学习在MATLAB中怎样表示离散时间信号。 3. 熟悉离散时间信号的作图。 二. 实验步骤

1. 复习离散时间系统的有关内容。 2. 复习MATLAB的基本语法。 3. 按实验内容熟悉 stem。 4. 编写程序。

5. 输出结果，总结结论，按要求写出实验报告。 三. 实验内容

1．掌握stem函数

STEM(Y) plots the data sequence Y as stems from the x axis terminated with circles for the data value.

STEM(X,Y) plots the data sequence Y at the values specified in X. 例：t=[0:0.1:2]; x=cos(pi*t+0.6); stem(t,x); xn=[4,2,2,3,6,7]; stem(xn);

思考：STEM(Y)与STEM(X,Y)有什么不同？STEM与PLOT函数有什么不同？

2．掌握subplot函数

H = SUBPLOT(m,n,p), or SUBPLOT(mnp), breaks the Figure window into an m-by-n matrix of small axes, selects the p-th axes for the current plot, and returns the axis handle. The axes are counted along the top row of the Figure window, then the second row, etc. 例：

n1=0:3;x1=[1,1,1,1];subplot(221);stem(n1,x1);title('x1序列');

n2=0:7;x2=[1,2,3,4,4,3,2,1];subplot(222);stem(n2,x2);title('x2序列'); n3=0:7;x3=[4,3,2,1,1,2,3,4];subplot(223);stem(n3,x3);title('x3序列'); n4=0:7;x41=cos((pi/4)*n4);subplot(224);stem(n4,x41);title('x4序列'); 思考：subplot是怎样分配各个作图分区的顺序号的？

1

3．信号的运算

x1(n)?[1,0.7,0.4,0.1,0]，x2(n)?[0.1,0.3,0.5,0.7,0.9]，请作出x1(n)?x2(n)，x1(n)x2(n)的图形。

思考：假如x1(n)与x2(n)长度不同，序列的求和和乘积运算能否执行，结果怎样？ 程序：x1=[1,0.7,0.4,0.1,0] x2=[0.1,0.3,0.5,0.7,0.9]

subplot(121);stem(x1+x2);title('x1+x2'); subplot(122);stem(x1.*x2);title('x1*x2');

4．掌握freqs与freqz函数

freqs用于s域的频率响应的计算，f reqz用于z域的频率响应的计算。请参看两者的help文件。

例：建立2个函数文件，分别输入为系统函数的分子，分母，和角频率的最大值，输出为系统频率响应的db，幅度，相位，和角频率。参看下面的程序。 function[db,mag,pha,Omega]=freqs_m(b,a,Omega_Max) % s域频率响应的计算

Omega=[0:1:500]*Omega_Max/500; H=freqs(b,a,Omega); mag=abs(H);

db=20*log10((mag+eps)/max(mag));

2

pha=angle(H);

function[db,mag,pha,w]=freqz_m(b,a,w_Max); % z域频率响应的计算 w=[0:1:499]*w_Max/500; H=freqz(b,a,w); mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=angle(H); （1）若Ha?s??1，请画出该系统幅频和相频特性。

s2?5s?6程序： b=1

a=[1 5 6]

Omega_Max=2*pi

[db,mag,pha,Omega]=freqs_m(b,a,Omega_Max)

subplot(221);plot(Omega,mag);title('AF的幅度响应'); subplot(222);plot(Omega,db);title('AF的幅度响应db'); subplot(223);plot(Omega,pha);title('AF的相位响应');

1?z?1（2）若H?z??，请画出该的幅频和相频特性。 ?1?21?0.8z?0.64z程序： b=[1 1]

a=[1 0.8 0.64] w_Max=2*pi

[db,mag,pha,w]=freqz_m(b,a,w_Max);

subplot(221);plot(w,mag);title('AF的幅度响应'); subplot(222);plot(w,db);title('AF的幅度响应db'); subplot(223);plot(w,pha);title('AF的相位响应');

3

思考：freqs，freqz的输入的各个参数与系统函数怎么联系起来，各种调用形式之间有什么不同点？观察上述系统的幅频和相频曲线，能得出什么结论？ 四. 实验报告要求

1. 简述实验目的及原理。

2. 概括各函数的常用调用方式，记录程序运行的各种结果。 3. 回答思考题。

4

实验二 用FFT作谱分析

一. 实验目的

1. 进一步加深DFT算法原理和基本性质的理解（FFT是DFT的一种快速算法，FFT的运算结果必然满足DFT的基本性质）。

2. 熟悉FFT算法原理和FFT 函数的应用。

3. 学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因。 二. 实验步骤

1. 复习DFT的定义、性质和用 DFT作谱分析的有关内容。

2. 复习FFT算法原理与编程思想。熟悉FFT算法的MATLAB实现。

注：MATLAB提供fft函数来计算x?n?的DFT，fft函数是用机器语言，而不是以MATLAB指令写成的，因此执行速度快。

格式（1）： y=fft(x) 计算x的FFT变换y。当x为矩阵，计算x中每一列信号的离散傅氏变换。当x的长度为2的幂时，采用基2算法，否则采用分裂基算法。

格式（2）： y=fft(x,n) 计算x的n点FFT，当x长度大于n时，截断x，否则补零。 Plot线性绘图函数。stem：绘制离散序列图。subplot：多坐标设置与定位当前坐标系。figure：创建新的图形窗口（用于输出图形的窗口）。 3. 产生下列信号并进行谱分析。 （1）x1?n??R4?n?

?n?1,0?n?3?（2）x2?n???8?n4?n?7

?0其他n??4?n,0?n?3?（3）x3?n???n?34?n?7

?0其他n?（4）x4?n??cos（5）x5?n??sin?4n n

?8（6）x6?t??cos8?t?cos16?t?cos20?t 4. 编写程序

5. 输出结果，总结结论，按要求写出实验报告。 三. 实验内容

1. 对上述6个信号，逐个进行谱分析

5

对于x1?n?，x2?n?，x3?n?，x4?n?，x5?n? ：N?8,16。对于x6?n? ：在一个周期内取N?8,16点抽样，做频谱分析。

（参考）每个信号谱分析的流程设计为：

开始

读入长度N

产生实验信号

绘制时间波形图

用fft函数进行谱分析

曲线绘制X?k?

结束

程序： % fft figure(1);

n1=0:3;x1=[1,1,1,1];subplot(221);stem(n1,x1);title('x1序列');

k1=0:7;y11=fft(x1,8);magy11=[abs(y11)];subplot(222);stem(k1,magy11);title('x1的8点FFT');

k2=0:15;y12=fft(x1,16);magy12=[abs(y12)];subplot(224);stem(k2,magy12);title('x1的16点FFT');

%x2 figure(2);

n2=0:7;x2=[1,2,3,4,4,3,2,1];subplot(221);stem(n2,x2);title('x2序列');

k1=0:7;y21=fft(x2,8);magy21=abs(y21);subplot(222);stem(k1,magy21);title('x2的8点FFT');

k2=0:15;y22=fft(x2,16);magy22=abs(y22);subplot(224);stem(k2,magy22);title('x2的16

6

点FFT');

%x3 figure(3);

n3=0:7;x3=[4,3,2,1,1,2,3,4];subplot(221);stem(n3,x3);title('x3序列');

k1=0:7;y31=fft(x3,8);magy31=abs(y31);subplot(222);stem(k1,magy31);title('x3的8点FFT');

k2=0:15;y32=fft(x3,16);magy32=abs(y32);subplot(224);stem(k2,magy32);title('x3的16点FFT');

%x4 figure(4);

n41=0:7;x41=cos((pi/4)*n41);subplot(221);stem(n41,x41);title('x4的8点序列');

k1=0:7;y41=fft(x41,8);magy41=abs(y41);subplot(222);stem(k1,magy41);title('x4的8点FFT');

n42=0:15;x42=cos((pi/4)*n42);subplot(223);stem(n42,x42);title('x4的16点序列');

k2=0:15;y42=fft(x42,16);magy42=abs(y42);subplot(224);stem(k2,magy42);title('x4的16点FFT');

7

%x5 figure(5);

n51=0:7;x51=sin((pi/8)*n51);subplot(221);stem(n51,x51);title('x5的8点序列');

k1=0:7;y51=fft(x51,8);magy51=abs(y51);subplot(222);stem(k1,magy51);title('x5的8点FFT');

n52=0:15;x52=sin((pi/8)*n52);subplot(223);stem(n52,x52);title('x5的16点序列');

k2=0:15;y52=fft(x52,16);magy52=abs(y52);subplot(224);stem(k2,magy52);title('x5的16点FFT');

%x6 figure(6);

f1=4;f2=8;f3=10;t=0:0.01:1;

x6=cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t);subplot(221);plot(t,x6);title('x6模拟信号');

%x61 8 T=0.5/8; t=0:T:0.5;

x6=cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t); k61=0:7;

y61=fft(x6,8);magy61=abs(y61);subplot(222);stem(k61,magy61);title('x6的8点FFT'); %x62 16 T=0.5/16; t=0:T:0.5;

8

x6=cos(f1*2*pi*t)+cos(f2*2*pi*t)+cos(f3*2*pi*t); k62=0:15;

y62=fft(x6,16);magy62=abs(y62);subplot(224);stem(k62,magy62);title('x6的16点FFT');

2. 令x7?n??x4?n??x5?n?，求出其N?16的FFTX7?k?。

当已求出X7?k?，根据DFT的对称性，求出X4?k??FFT?x4?n??和

X5?k??FFT?x5?n??，与起第1步求出的X4?k?、X5?k?比较。

程序： %x7 figure(7);

n7=0:15;x7=x42+x52;subplot(221);stem(n7,x7);title('x7的16点序列');

k7=n7;y7=fft(x7,16);magy7=abs(y7);subplot(222);stem(k7,magy7);title('x7的16点FFT');

X4=real(y7);X5=imag(y7)*i;subplot(223);stem(k7,abs(X4));title('x4的16点序列'); subplot(224);stem(k7,abs(X5));title('x5的16点序列');

3. 令x8?n??x4?n??jx5?n?，求出其N?16的FFTX8?k?。

当已求出X8?k?，根据DFT的对称性，求出X4?k??FFT?x4?n??和

9

X5?k??FFT?x5?n??，与起第1步求出的X4?k?、X5?k?比较。

程序： %x8 figure(8);

n8=0:15;x8=x42+x52*j;subplot(221);stem(n8,abs(x8));title('x8的16点序列');

k8=n8;y8=fft(x8,16);magy8=abs(y8);subplot(222);stem(k8,magy8);title('x8的16点FFT');

tempk=mod((16-k8),16); tempy8=y8(tempk+1);

convy8=real(tempy8)-imag(tempy8)*i;

X4=0.5*(y8+convy8);subplot(223);stem(k8,abs(X4));title('x4的16点序列'); X5=1/(2*j)*(y8-convy8);subplot(224);stem(k8,abs(X5));title('x5的16点序列');

四. 思考题

1. 在N?8时，x2?n?和x3?n?的幅频特性会相同吗？为什么？N?16呢？

2. 如果周期信号的周期预先不知道，如何使用FFT进行分析。 五. 实验报告要求

1. 简述实验的原理和目的。

2. 输出实验中所用的程序以及所得的时域频域性曲线，与理论结果比较，分析误差产生的原因，以及FFT做谱分析时有关参数的选择方法。 3. 总结实验所得的主要结论。 4. 简要回答思考题。

10

db=20*log10((mag+eps)/max(mag)); pha=angle(H);

subplot(222);plot(w/pi,mag);title('实际的幅度响应'); subplot(224);plot(w/pi,pha);title('实际的相位响应'); figure;

plot(w/pi,db);title('实际的幅度db响应');

db=20*log10((mag+eps)/max(mag)); figure; N=33;

wc=0.25*pi; n=0:1:N-1;

hd=ideal_lp(wc,N);

subplot(221);stem(n,hd);title('理想单位脉冲响应'); w_ham=(hamming(N))'; h=hd.*w_ham;

subplot(223);stem(n,h);title('实际单位脉冲响应'); [H,w]=freqz(h)

mag=abs(H);

db=20*log10((mag+eps)/max(mag)); pha=angle(H);

subplot(222);plot(w/pi,mag);title('实际的幅度响应'); subplot(224);plot(w/pi,pha);title('实际的相位响应'); figure;

plot(w/pi,db);title('实际的幅度db响应');

ideal_lp函数：

function hd=ideal_lp(wc,N)

%N为奇数，理想低通滤波器的脉冲响应h(n) alpha=(N-1)/2; n=0:1:N-1; m=n-alpha;

hd=sin(wc*m)./(pi*m); hd(alpha+1)=wc/pi;

16

（2）N=33，?c?0.25?，用四种窗函数(矩形窗、汉宁窗、海明窗、布拉克曼窗)设计线性相位低通滤波器，绘制相应的幅频特性曲线（dB），观察3dB和20dB带宽以及阻带的衰减，比较四种窗函数对滤波器特性的影响。

程序：

%不同的窗口对滤波器的影响 N=33;

wc=0.25*pi;

hd=ideal_lp(wc,N);

w_box=(boxcar(N))';%矩形窗 h_box=hd.*w_box;

w_han=(hanning(N))';%汉宁窗 h_han=hd.*w_han;

w_ham=(hamming(N))';%海明窗 h_ham=hd.*w_ham;

w_bla=(blackman(N))';%二阶升余弦窗 h_bla=hd.*w_bla; [H,w]=freqz(h_box);

subplot(221);plot(w/pi,20*log10(abs(H)));title('由矩形窗设计的filter幅度响应'); [H,w]=freqz(h_han);

subplot(222);plot(w/pi,20*log10(abs(H)));title('由汉宁窗设计的filter幅度响应'); [H,w]=freqz(h_ham);

17

subplot(223);plot(w/pi,20*log10(abs(H)));title('由海明窗设计的filter幅度响应'); [H,w]=freqz(h_bla);

subplot(224);plot(w/pi,20*log10(abs(H)));title('由二阶升余弦窗设计的filter幅度响应');

2. 实验步骤

参考下列程序流程图：

开始

读入窗口长度

计算hd?n?

调用窗函数求w?n? 计算h?n??hd?n?w?n?

计算频响Hdej?

绘制幅频和相频曲线

结束

四. 思考题

1. 如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器？写出设计步骤。

2. 如果要求用窗函数法设计带通滤波器，给定上下边带截止频率，求出理想带通的单位脉冲响应，并试写出设计数字滤波器的步骤。

3. 若用频率抽样法设计N=33的FIR滤波器，应如何设计？简要写出设计思路。 五. 实验报告要求

1. 简述实验目的和原理。

2. 按照实验步骤及要求，比较各种情况下的滤波性能，说明窗口长度N和窗函数类型对滤波器特性的影响。

3. 总结用窗函数法设计FIR滤波器的主要特点。

?? 18

本文来源：https://www.bwwdw.com/article/g3xo.html