数字信号处理实验报告一

数字信号处理

第一次实验报告

学 院：信息工程学院

班 级：2012级电子信息工程*班 姓 名：

学 号：********** 指导老师：

1

实验一：系统响应及系统稳定性

一、实验原理简述

1、在时域求系统响应的方法有两种，第一种是通过解差分方程求得系统输出；第二种是已知系统的单位脉冲响应，通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。

2、实验中检查系统是否稳定，不可能检查系统对所有有界的输入信号，输出是否都是有界输出，或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列，如果系统的输出趋于一个常数（包括0），就可断定系统是稳定的。系统的稳态输出时指当n—>无穷大时，系统的输出。如果系统稳定，信号加入系统后，系统输出的开始一段称为暂态效应。随n的增大，达到稳态输出。 二、程序与图形

1、%-----（2）调用filter解差分方程以及单位脉冲响应------ close all;clear all

A=[1,-0.9];B=[0.05,0.05]; %系统差分方程系数向量B和A

x1n=[ones(1,8) zeros(1,25)]; %产生信号x1(n)=R8(n)，用zeros用来加点的个数 x2n=ones(1,30); %产生信号x2(n)=u(n)

hn=impz(B,A,25); %求系统单位脉冲响应h(n) subplot(3,1,1);stem(hn); %调用函数stem绘图 title('(a) 系统单位脉冲响应h(n)');

y1n=filter(B,A,x1n); %求系统对x1(n)的响应y1(n) subplot(3,1,2);stem(y1n);

title('(b) 系统对R8(n)的响应y1(n)');

y2n=filter(B,A,x2n); %求系统对x2(n)的响应y2(n) subplot(3,1,3);stem(y2n);

title('(c) 系统对u(n)的响应y2(n)');

(a) 系统单位脉冲响应h(n)0.10.05010.5010.500510152025(b) 系统对R8(n)的响应y1(n)05101520253035(c) 系统对u(n)的响应y2(n)0510152025302

%-----（3）调用conv函数计算卷积------- x1n=[ones(1,8)]; %产生信号x1(n)=R8(n) h1n=[ones(1,10) zeros(1,10)]; h2n=[1 2.5 2.5 1 zeros(1,10)]; y21n=conv(h1n,x1n); y22n=conv(h2n,x1n); subplot(2,2,1);stem(h1n);

title('(d) 系统单位脉冲响应h1(n)'); subplot(2,2,3);stem(y21n);

title('(e) h1(n)与R8(n)的卷积y21(n)'); subplot(2,2,2);stem(h2n);

title('(f) 系统单位脉冲响应h2(n)'); subplot(2,2,4);stem(y22n);

title('(g) h2(n)与R8(n)的卷积y22(n)');

(d) 系统单位脉冲响应h1(n)1320.51(f) 系统单位脉冲响应h2(n)0051015200051015(e) h1(n)与R8(n)的卷积y21(n)86420010203086420(g) h2(n)与R8(n)的卷积y22(n)0102030

%-----（4）实验方法检查系统是否稳定------- close all;clear all

un=ones(1,256); %产生信号u(n) n=0:255;

xsin=sin(0.014*n)+sin(0.4*n); %产生正弦信号

A=[1,-1.8237,0.9801];B=[1/100.49,0,-1/100.49]; %系统差分方程系数向量B和A y1n=filter(B,A,un); %谐振器对u(n)的响应y31(n) y2n=filter(B,A,xsin); %谐振器对u(n)的响应y31(n) subplot(2,1,1);stem(y1n);

title('(h) 谐振器对u(n)的响应y31(n)');

3

subplot(2,1,2);stem(y2n);

title('(i) 谐振器对正弦信号的响应y32(n)');

(h) 谐振器对u(n)的响应y31(n)0.050-0.05050100150200250300(i) 谐振器对正弦信号的响应y32(n)10.50-0.5-1050100150200250300

三、实验结果分析

（a）25个点数和程序所写一致。

Filter函数实现线性常系数差分方程的递推求解。 调用格式如下：

Y=[filter(B,A,x)] 计算系统对输入信号x的零状态响应输出信号向量Y，B,A是差分方程的系数向量。即B=[a1,a2……am] A=[b1,b2……bn] （d）（f）单位脉冲响应点数与程序要求一致 （e）（g）卷积点数满足M+N-1的要求，图形也满足要求。 Conv函数用于计算两个有限长序列的卷积

C=conv（A,B）计算两个有限长序列向量A和B的卷积

（h）中由：在系统的输入端加入单位阶跃序列，如果系统的输出趋近一个常数（包括零），就可以断定系统是稳定的输出显然趋近于零，所以是稳定的

（i）中谐振器具有对某个频率进行谐振的性质，本实验中的谐振器的谐振频率是0.4 rad,因此稳定波形为sin(0.4n)。 四、思考题简答

（1） 如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可否用线性卷积法求系统的响应? 如何求？

答：如果输入信号为无限长序列，系统的单位脉冲响应是有限长序列，可用分段线性卷积法求系统的响应 （2）如果信号经过低通滤波器，把信号的高频分量滤掉，时域信号会有何变化，用前面第一个实验结果进行分析说明。

答：如果信号经过低通滤波器， 则信号的高频分量被滤掉， 时域信号的变化减缓，在有阶跃处附近产生过渡带。因此，当输入矩形序列时，输出序列的开始和

4

终了都产生了明显的过渡带。

实验二：时域采样与频域采样

一、实验原理简述 1、时域采样定理：

a) 对模拟信号xa(t)以间隔T进行时域等间隔理想采样，形成的采样信号的

ˆ(j)是原模拟信号频谱X(j)以采样角频率（2/T）频谱Xssa为周期进行周期延拓。公式为：

1ˆ(j)FT[xˆa(t)]Xa(jjns) XaTnb) 采样频率s必须大于等于模拟信号最高频率的两倍以上，才能使采样信号的频谱不产生频谱混叠。 2、 频域采样定理：

公式为：xN(n)IDFT[XN(k)]N[x(niN)]RN(n)

i由公式可知，频域采样点数N必须大于等于时域离散信号的长度M(即N≥M)，才能使时域不产生混叠，则N点IDFT[XN(k)]得到的序列xN(n)就是原序列x(n),即xN(n)=x(n)。 二、程序与图形

% ---------（1）时域采样理论验证程序--------- Tp=64/1000; %观察时间Tp=64毫秒 %产生M长采样序列x(n) Fs=1000; T=1/Fs; M=Tp*Fs; n=0:M-1;

A=444.128; alph=pi*50*2^0.5 ;omega=pi*50*2^0.5; xnt=A*exp(-alph*n*T).*sin(omega*n*T); Xk=T*fft(xnt,M); %M点FFT[xnt)] yn='xa(nT)';subplot(3,2,1);

stem(xnt); %调用自编绘图函数stem绘制序列图 box on;title('(a) Fs=1000Hz'); k=0:M-1;fk=k/Tp;

subplot(3,2,2);stem(fk,abs(Xk));title('(a) T*FT[xa(nT)],Fs=1000Hz'); xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);

% Fs=300Hz和 Fs=200Hz的程序与上面Fs=1000Hz完全相同。

Tp=64/1000; %观察时间Tp=64毫秒 %产生M长采样序列x(n)

5

Fs=300; T=1/Fs; M=Tp*Fs; n=0:M-1;

A=444.128; alph=pi*50*2^0.5 ;omega=pi*50*2^0.5; xnt=A*exp(-alph*n*T).*sin(omega*n*T);

Xk=T*fft(xnt,M); %M点FFT[xnt)] yn='xa(nT)';subplot(3,2,3); stem(xnt); %调用自编绘图函数stem绘制序列图 box on;title('(b) Fs=300Hz'); k=0:M-1;fk=k/Tp;

subplot(3,2,4);stem(fk,abs(Xk));title('(b) T*FT[xa(nT)],Fs=300Hz'); xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);

Tp=64/1000; %观察时间Tp=64毫秒 %产生M长采样序列x(n) Fs=200; T=1/Fs; M=Tp*Fs; n=0:M-1;

A=444.128; alph=pi*50*2^0.5 ;omega=pi*50*2^0.5; xnt=A*exp(-alph*n*T).*sin(omega*n*T); Xk=T*fft(xnt,M); %M点FFT[xnt)] yn='xa(nT)';subplot(3,2,5);

stem(xnt); %调用自编绘图函数stem绘制序列图 box on;title('(c) Fs=200Hz'); k=0:M-1;fk=k/Tp;

subplot(3,2,6);stem(fk,abs(Xk));title('(c) T*FT[xa(nT)],Fs=200Hz'); xlabel('f(Hz)');ylabel('幅度');axis([0,Fs,0,1.2*max(abs(Xk))]);

(a) Fs=1000Hz20011000-1002000-2000510(c) Fs=200Hz2001520020406080(a) T*FT[xa(nT)],Fs=1000Hz幅度0.500(b) Fs=300Hz15001000f(Hz)(b) T*FT[xa(nT)],Fs=300Hz幅度0.500100200300f(Hz)(c) T*FT[xa(nT)],Fs=200Hz幅度0510150-2000.50050100f(Hz)150200

6

%-------------(2)频域采样理论验证----------- M=27;N=32;n=0:M;

%产生M长三角波序列x(n)

xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb];

Xk=fft(xn,1024); %1024点FFT[x(n)], 用于近似序列x(n)的TF X32k=fft(xn,32) ;%32点FFT[x(n)]

x32n=ifft(X32k); %32点IFFT[X32(k)]得到x32(n) X16k=X32k(1:2:N); %隔点抽取X32k得到X16(K)

x16n=ifft(X16k,N/2); %16点IFFT[X16(k)]得到x16(n) subplot(3,2,2);stem(n,xn);box on

title('(b) 三角波序列x(n)');xlabel('n');ylabel('x(n)');axis([0,32,0,20]) k=0:1023;wk=2*k/1024; %

subplot(3,2,1);plot(wk,abs(Xk));title('(a)FT[x(n)]');

xlabel('\\omega/\\pi');ylabel('|X(e^j^\\omega)|');axis([0,1,0,200]) k=0:N/2-1;

subplot(3,2,3);stem(k,abs(X16k));box on

title('(c) 16点频域采样');xlabel('k');ylabel('|X_1_6(k)|');axis([0,8,0,200]) n1=0:N/2-1;

subplot(3,2,4);stem(n1,x16n);box on title('(d)16点IDFT[X_1_6(k)]');

xlabel('n');ylabel('x_1_6(n)');axis([0,32,0,20]) k=0:N-1;

subplot(3,2,5);stem(k,abs(X32k));box on

title('(e) 32点频域采样');xlabel('k');ylabel('|X_3_2(k)|');axis([0,16,0,200]) n1=0:N-1;

subplot(3,2,6);stem(n1,x32n);box on title('(f)32点IDFT[X_3_2(k)]');

xlabel('n');ylabel('x_3_2(n)');axis([0,32,0,20])

(a)FT[x(n)]20020(b) 三角波序列x(n)|X(ej)|10002000-2000510(c) Fs=200Hz2001520x(n)00.5/(b) Fs=300Hz110002030n(b) T*FT[xa(nT)],Fs=300Hz101幅度0.500100200300f(Hz)(c) T*FT[xa(nT)],Fs=200Hz幅度0510150-2000.50050100f(Hz)1502007

三、实验结果分析

第一个程序中当采样频率为1000 Hz时，频谱混叠很小：当采样频率为300 Hz时，频谱混叠很严重；当采样频率为200 Hz时，频谱混叠更很严重。 fft函数的调用格式：Xk=fft(xn,N)

调用参数xn为被交换的时域序列向量，N是DFT变换的区间长度，当N大于xn的长度时，fft函数自动在xn后面补零。当N小于xn的长度时，fft函数计算xn的前面N个元素构成的N长序列的N点DFT,忽略xn后面的元素。

第二个程序验证了频域采样理论和频域采样定理。 对信号x(n)的频谱函数X(ejω)在［0， 2π］上等间隔采样N=16时， N点IDFT［XN(k)］得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列 当N=16时，由于N当N=32时，由于N>M，频域采样定理，所以不存在时域混叠失真ifft函数用法同fft函数 四、思考题简答

jX(e)在[0,2]上的N点等间隔1、 如果序列x(n)的长度为M，希望得到其频谱

采样，当N答：先对原序列x(n)以N为周期进行周期延拓后取主值区序列，

xN(n)[x(niN)]RN(n)i

再计算N点DFT则得到N点频域采样：

XN(k)DFT[xN(n)]N =X(ej)2kN , k0,1,2,,N1

实验三：用FFT对信号作频谱分析

一、 实验原理简述

1、周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。如果不知道信号周期，可以尽量选择信号的观察时间长一些。

2、对模拟信号进行谱分析时，首先要按照采样定理将其变成时域离散信号。如果是模拟周期信号，也应该选取整数倍周期的长度，经过采样后形成周期序列，按照周期序列的谱分析进行。

3、对信号进行谱分析的重要问题是频谱分辨率D和分析误差。频谱分辨率直接和FFT的变换区间N有关，因为FFT能够实现的频率分辨率是2/N，因此要求2/ND。可以根据此式选择FFT的变换区间N。误差主要来自于用FFT作频谱分析时，得到的是离散谱，而信号（周期信号除外）是连续谱，只有当N较大时离散谱的包络才能逼近于连续谱，因此N要适当选择大一些。 二、程序与图形

%-----子程序-----% function mstem(Xk) M=length(Xk);

k=0:M-1; wk=2*k/M;

8

stem(wk,abs(Xk),'.'); box on; xlabel('ω/π'); ylabel('幅度'); axis([0,2,0,1.2*max(abs(Xk))]);

clear all;close all

%---------------实验内容(1)----------------

x1n=[ones(1,4)]; %产生序列向量x1(n)=R4(n)

M=8; xa=1:(M/2); xb=(M/2):-1:1; x2n=[xa,xb]; %产生长度为8的三角波序列x2(n) x3n=[xb,xa];

X1k8=fft(x1n,8); %计算x1n的8点DFT X1k16=fft(x1n,16); %计算x1n的16点DFT X2k8=fft(x2n,8); %计算x1n的8点DFT X2k16=fft(x2n,16); %计算x1n的16点DFT X3k8=fft(x3n,8); %计算x1n的8点DFT X3k16=fft(x3n,16); %计算x1n的16点DFT %以下绘制幅频特性曲线 subplot(3,2,1);

mstem(X1k8); %绘制8点DFT的幅频特性图

title('(1a) 8点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X1k8))])

subplot(3,2,2);mstem(X1k16); %绘制16点DFT的幅频特性图 title('(1b)16点DFT[x_1(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X1k16))])

subplot(3,2,3);mstem(X2k8); %绘制8点DFT的幅频特性图 title('(2a) 8点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X2k8))])

subplot(3,2,4);mstem(X2k16); %绘制16点DFT的幅频特性图 title('(2b)16点DFT[x_2(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X2k16))])

subplot(3,2,5);mstem(X3k8); %绘制8点DFT的幅频特性图 title('(3a) 8点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X3k8))])

subplot(3,2,6);mstem(X3k16); %绘制16点DFT的幅频特性图 title('(3b)16点DFT[x_3(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X3k16))])

9

(1a) 8点DFT[x1(n)]44(1b)16点DFT[x1(n)]幅度2000.511.5ω/π(2a) 8点DFT[x2(n)]2幅度20011.5ω/π(2b)16点DFT[x2(n)]0.522020幅度10000.511.5ω/π(3a) 8点DFT[x3(n)]2幅度100011.5ω/π(3b)16点DFT[x3(n)]0.522020幅度10000.51ω/π1.52幅度10000.51ω/π1.52

%-------------实验内容（2）-----------

N=8;n=0:N-1; %FFT的变换区间N=8 x4n=cos(pi*n/4);

x5n=cos(pi*n/4)+cos(pi*n/8);

X4k8=fft(x4n); %计算x4n的8点DFT X5k8=fft(x5n); %计算x5n的8点DFT N=16;n=0:N-1; %FFT的变换区间N=16 x4n=cos(pi*n/4);

x5n=cos(pi*n/4)+cos(pi*n/8);

X4k16=fft(x4n); %计算x4n的16点DFT X5k16=fft(x5n); %计算x5n的16点DFT

subplot(2,2,1);mstem(X4k8); %绘制8点DFT的幅频特性图 title('(a) 8点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X4k8))])

subplot(2,2,3);mstem(X4k16); %绘制16点DFT的幅频特性图 title('(b)16点DFT[x_4(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X4k16))])

subplot(2,2,2);mstem(X5k8); %绘制8点DFT的幅频特性图 title('(a) 8点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X5k8))])

subplot(2,2,4);mstem(X5k16); %绘制16点DFT的幅频特性图 title('(b)16点DFT[x_5(n)]');xlabel('ω/π');ylabel('幅度'); axis([0,2,0,1.2*max(abs(X5k16))])

10

(a) 8点DFT[x4(n)]46(a) 8点DFT[x5(n)]幅度21000.511.5ω/π(b)16点DFT[x4(n)]2幅度342000.511.5ω/π(b)16点DFT[x5(n)]288幅度642000.51ω/π1.52幅度642000.51ω/π1.52

%----------实验内容（3）--------- Fs=64;T=1/Fs;

N=16;n=0:N-1; %FFT的变换区间N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)16点采样 X6k16=fft(x6nT); %计算x6nT的16点DFT

X6k16=fftshift(X6k16); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心） subplot(3,1,1);stem(fk,abs(X6k16),'.');box on %绘制8点DFT的幅频特性图 title('(6a) 16点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k16))]) N=32;n=0:N-1; %FFT的变换区间N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)32点采样 X6k32=fft(x6nT); %计算x6nT的32点DFT

X6k32=fftshift(X6k32); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心） subplot(3,1,2);stem(fk,abs(X6k32),'.');box on %绘制8点DFT的幅频特性图 title('(6b) 32点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k32))]) N=64;n=0:N-1; %FFT的变换区间N=16

x6nT=cos(8*pi*n*T)+cos(16*pi*n*T)+cos(20*pi*n*T); %对x6(t)64点采样 X6k64=fft(x6nT); %计算x6nT的64点DFT

X6k64=fftshift(X6k64); %将零频率移到频谱中心 Tp=N*T;F=1/Tp; %频率分辨率F

k=-N/2:N/2-1;fk=k*F; %产生16点DFT对应的采样点频率（以零频率为中心）

11

subplot(3,1,3);stem(fk,abs(X6k64),'.'); box on%绘制8点DFT的幅频特性图 title('(6a) 64点|DFT[x_6(nT)]|');xlabel('f(Hz)');ylabel('幅度'); axis([-N*F/2-1,N*F/2-1,0,1.2*max(abs(X6k64))])

(6a) 16点|DFT[x6(nT)]|10幅度50-30-20-10010(6b) 32点f(Hz)|DFT[x6(nT)]|2030幅度100-30-20-10010f(Hz)(6a) 64点|DFT[x6(nT)]|2030幅度200-30-20-100f(Hz)102030

二、 思考题简答

（1）对于周期序列，如果周期不知道，如何用FFT进行谱分析？

答：可先截取M点进行DFT,再将截取长度扩大1倍，比较两次的结果。如果二者的主谱差别满足分析误差要求，则以两者中的一个近似表示周期序列的频谱，否则，继续把截取长度加倍，重复上述步骤。

（2）如何选择FFT的变换区间？（包括非周期信号和周期信号）

答：对于非周期信号：有频谱分辨率F，而频谱分辨率直接和FFT的变换区间有关，因为FFT能够实现的频率分辨率是2π/N...因此有最小的N>2π/F。就可以根据此式选择FFT的变换区间。对于周期信号，周期信号的频谱是离散谱，只有用整数倍周期的长度作FFT，得到的离散谱才能代表周期信号的频谱。 （3）当N=8时，x2(n)和

x3(n)的幅频特性会相同吗？为什么？N=16 呢？

答 ：不同，因为这样会影响是不是周期的整数倍的问题，即影响了频谱的正确性。

12