matlab中怎么定义tstem函数,数字信号处理实验报告一系统响应及系统稳定性

数字信号实验报告

实验名称 : 数字信号处理第一次实验 指导老师 : 学院 : 信息工程学院 班级 : 12电子信息工程2班 姓名 :

学号 :

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

实验目的

1. 掌握求系统响应的方法

2. 掌握时域离散系统的时域特性

3. 分析观察及检验系统的稳定性

实验原理

在时域中,描写系统特性的方法是差分方程和单位脉冲响应,在频域可以用系统函数描述系统特性。已知输入信号可以由差分方程/单位脉冲响应或系统函数求出系统对于该输入信号的响应。本实验仅在时域求解。在计算机上适合用递推法求差分方程的解,最简单的方法是采用MA TLAB 语言工具箱函数filter 函数。也可以用MA TLAB 语言的工具箱函数conv 函数计算输入信号和系统的单位脉冲响应的线性卷积,求出系统的响应。系统的时域特性指的是系统的线性时不变性质/因果性和稳定性。重点分析实验系统的稳定性,包括观察系统的暂态响应和稳定响应。系统的稳定性是指对任意有界的输入信号,系统都能得到有界的系统响应。或者系统的单位脉冲响应满足绝对可和的条件。系统的稳定性由其差分方程的系数决定。实际中检查系统是否稳定,不可能检查系统对所有有界的输入信号,输出是否都是有界输出,或者检查系统的单位脉冲响应满足绝对可和的条件。可行的方法是在系统的输入端加入单位阶跃序列,如果系统的输出趋近于一个常数(包括零),就可以断定系统是稳定的。系统的稳态输出是指当n 趋向于无穷时,系统的输出。如果系统稳定,则信号加入系统后,系统输出的开始一段称为暂态效应,随着n 的加大,幅度趋于稳定,达到稳态输出。

实验内容及步骤

1. 编制程序,包括产生输入信号、单位脉冲响应序列的子程序,用filter 函数或conv 函数求解系统输出响应的主程序。程序中要有绘制信号波形的功能。 2. 给定一个地通滤波器的差分方程为y (n )=0. 05x (n )+0. 05x (n -1) +0. 9y (n -1) 输入信号:

x 1(n ) =R 8(n ) x 2(n ) =μ(n )

①分别求出x 1(n ) =R 8(n ) 和x 2(n ) =μ(n ) 的系统响应y 1(n ) 和y 2(n ) ,并画出其

波形

②求出系统的单位脉冲响应,画出其波形。

(1) 给定系统的单位脉冲响应为 h 1(n ) =R 10(n )

h 2(n ) =δ(n ) +0. 25δ(n -1) +2. 5δ(n -2) +δ(n -3)

用线性卷积法求x 1(n ) =R 8(n ) 分别对系统h 1(n ) 和h 2(n ) 的输出响应y 21(n ) 和y 22(n ) 并画出波形。

(2) 给定一谐振器的差分方程为

谐振器y (n ) =1. 8237y (n -1) -0. 9801y (n -2) +b 0x (n ) -b 0x (n -2) 令b 0=1/100. 49,

的谐振频率为0.4rad

①用实验方法检查系统是否稳定,输入信号为μ(n ) 时,画出系统输出波形y 31(n ) 。 ②给定输入信号为x (n ) =sin(0. 014n ) +sin(0. 4n ) ,求出系统的输出响应y 32(n ) ,并画出其波形。

实验程序

实验1—1程序:

close all ; clear all ; A=[1,-0.9]; B=[0.05,0.05];

x1n=[1 1 1 1 1 1 1 1 zeros(1,50)]; x2n=ones(1,128); hn=impz(B,A,58); subplot(2,2,1); y='h(n)'; tstem(hn,y);

title('(a)系统单位脉冲响应h(n)'); box on

y1n=filter(B,A,x1n); subplot(2,2,2);y='y1(n)'; tstem(y1n,y);

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

y2n=filter(B,A,x2n); subplot(2,2,4);y='y2(n)'; tstem(y2n,y);

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

波形图:

(a) 系统单位脉冲响应

h(n)(b) 系统对R8(n)的响应y1(n)

y n

n

y n

n

(c) 系统对u(n)的响应y2(n)y n

n

实验1—2程序:

close all ; clear all ;

x1n=[1 1 1 1 1 1 1 1 ]; 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); figure(2) subplot(2,2,1); y='h1(n)'; tstem(h1n,y);

title('(d) 系统单位脉冲响应h1(n)'); box on

subplot(2,2,2);y='y21(n)'; tstem(y21n,y);

title('(e) h1(n)与R8(n)的卷积y21(n)'); box on

subplot(2,2,3);y='h2(n)'; tstem(h2n,y);

title('(f) 系统单位脉冲响应h2(n)'); box on

subplot(2,2,4);

y='y22(n)'; tstem(y22n,y);

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

波形图:

(d) 系统单位脉冲响应h1(n)

(e) h1(n)与R8(n)的卷积y21(n)

y n

y n

n

(f) 系统单位脉冲响应

h2(n)

n

(g) h2(n)与R8(n)的卷积y22(n)y n

y n

n

n

实验1—3程序:

close all; clear all; un=ones(1,256); 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]; y31n=filter(B,A,un); y32n=filter(B,A,xsin); figure(3) subplot(2,1,1); y='y31(n)'; tstem(y31n,y);

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

subplot(2,1,2);

y='y32(n)'; tstem(y32n,y);

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

波形图:

(h) 谐振器对u(n)的响应

y31(n)

y n

n

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

y n

n

分析讨论:

(1) 综合起来,在时域求系统响应的方法有两种,第一是通过解差分方程求得系统输出,

注意合理地选择初始条件;第二种是已知系统的单位脉冲响应,通过求输入信号和系统单位脉冲响应的线性卷积求得系统输出。用计算机求解时最好使用MATLAB 语言进行。

(2) 实际中要检验系统的稳定性,其方法是在输入端加入单位阶跃序列,观察输出波形,

如果波形稳定在一个常数值上,系统稳定,否则不稳定。如第三个实验就是稳定的。

(3) 谢正奇具有对某个频率进行谐振的性质,本实验中的谐振器的谐振频率是0.4rad ,

因此稳定波形为sin(0.4n)。

思考题

思考题1-1

如果输入信号为无限长序列,系统的单位脉冲响应是有限长序列,可用分段线性卷积法求系统的响应。

思考题1-2

如果信号经过经过低通滤波器,则信号的高频分量被滤掉,时域信号的变化减缓,在有阶跃处附近产生过渡变化时间。因此,当输入矩形序列时,输出序列的开始和终了都产生了明显的上升和下降时间,详细可见第一个实验结果的波形。

实验报告要求

(1) (2) (3) (4) (5)

简述在时域求系统响应的方法

简述通过实验判断系统稳定性的方法。分析上面第三个实验的稳定输出的波形。 对各实验所得的结果进行简单的分析和解释。 简要回答思考题

打印程序清单和要求的各信号波形。

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

实验目的

时域采样理论与频域采样理论是数字信号处理中的重要理论。要求掌握模拟信号采样前后的频谱变化,以及如何选择采样频率才能使采样后的信号不丢失信息;要求掌握频域采样会引起时域周期化的概念,以及频率域采样定理及其对频域采样点数选择的指导作用。

实验原理

1. 时域采样定理的要点

(1)对模拟信号x a (t ) 以T 进行时域等间隔理想采样,形成的采样信号的频谱X (j Ω) 会以采样角频率Ωs (Ωs =2π/T ) 为周期进行周期延拓。公式为

1∞

X a (j Ω) =FT [x a (t ) ]=∑X a (j Ω-jn Ωs )

T n =-∞

(2)采样频率Ωs 必须大于等于模拟信号最高频率的两倍以上,才能使采样信号的频谱不产生频谱混叠。

利用计算机计算X a (j Ω) 并不方便,下面我们导出另外一个公式,以便在计算机上进行实验。

理想采样信号x a (t ) 和模拟信号x a (t ) 之间的关系为

x a (t ) =x a (t ) ∑δ(t -nT ) 对上式进行傅立叶变换,得到

n =-∞

X a (j Ω) =⎰[x a (t ) ∑δ(t -nT ) ]e

-∞

n =-∞

-j Ωt

dt =

n =-∞

∑⎰

-∞

x a (t ) δ(t -nT ) e -j Ωt dt

在上式的积分号内只有当t=nT时,才有非零值,因此

X a (j Ω) =

n =-∞

∑x (nT ) e

a

-j ΩnT

上式中,在数值上x a (nT ) =x (n ) ,再将ω=Ω代入,得到

X a (j Ω) =

n =-∞

∑x (n ) e

-jwn

上式的右边就是序列的傅立叶变换X (e ) ,即

jw

X a (j Ω) =X (e jw ) |w =ΩT

上式说明理想 的傅立叶变换可用相应的采样序列的傅立叶变换得到,只要将自变量ω用ΩT 代替即可。 2. 频域采样定理的要点

(1)对信号x(n)的频谱函数X (e ) 在[0,2π]上等间隔采样N 点,得到

jw

X N (k ) =X (e jw ) |

w =

2k πN

( k=0,1,2…N-1),则N 点IDFT [X N (k )]得到的序列就是原序列x(n)以N 为周期进行周期延拓后的主值区序列,公式x N (n ) =IDFT [X N (k )]N =[

i =-∞

∑x (n +iN ) ]R

N

(n )

(2)由上式可知,频率采样点数N 必须大于等于时域离散信号的长度M (即N ≥M ),才能使时域不产生混叠,这时N 点

[IDFT (X N k )]得到的序列x N (n ) 就是原序列x (n ) ,即

x N (n ) =x (n ) 。如果N>M,则x N (n ) 比原序列尾部多了N-M 个零点;如果N

而且N

因此

x N (n ) 与x (n ) 不相同。

实验程序

实验2—1程序: close all; clear all; Tp=64/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); yn='xa(nT)'; subplot(3,2,1); tstem(xnt,yn); box on;

title('(a) Fs=1000Hz'); k=0:M-1; fk=k/Tp;

subplot(3,2,2); plot(fk,abs(Xk));

title('(a) T*FT[xa(nT)],Fs=1000Hz'); xlabel('f(Hz)');ylabel('幅度'); axis([0,Fs,0,1.2*max(abs(Xk))]) Tp=64/300; 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); yn='xa(nT)'; subplot(3,2,3); tstem(xnt,yn); box on;

title('(a) Fs=300Hz'); k=0:M-1; fk=k/Tp;

subplot(3,2,4); plot(fk,abs(Xk));

title('(a) T*FT[xa(nT)],Fs=300Hz'); xlabel('f(Hz)');ylabel('幅度'); axis([0,Fs,0,1.2*max(abs(Xk))]) Tp=64/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); yn='xa(nT)'; subplot(3,2,5); tstem(xnt,yn); box on;

title('(a) Fs=200Hz'); k=0:M-1; fk=k/Tp;

subplot(3,2,6); plot(fk,abs(Xk));

title('(a) T*FT[xa(nT)],Fs=200Hz'); xlabel('f(Hz)');

ylabel('幅度');

axis([0,Fs,0,1.2*max(abs(Xk))]) 波形图:

(a) Fs=1000Hz

n

(a) Fs=300Hz

n

(a) Fs=200Hz

11

(a) T*FT[xa(nT)],Fs=1000Hz

幅度

y n

0.500

5001000f(Hz)

(a) T*FT[xa(nT)],Fs=300Hz

幅度

y n

0.500

100

200

300

f(Hz)

(a) T*FT[xa(nT)],Fs=200Hz

幅度

n

y n

0.500

50

100f(Hz)

150

200

实验2—2程序:

M=27; N=32; n=0:M;

xa=0:floor(M/2); xb= ceil(M/2)-1:-1:0; xn=[xa,xb]; Xk=fft(xn,1024); X32k=fft(xn,32); x32n=ifft(X32k); X16k=X32k(1:2:N); x16n=ifft(X16k,N/2); 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)]

200

(b) 三角波序列x(n)

|X (e j ω) |

10000

0.5ω/π

(c) 16点频

域采样

1

x (n )

n

(d) 16点

IDFT[X16(k)]

|X 16(k ) |

k

x 16(n )

10

20

30

(e) 32点频域采样

n

(f) 32点IDFT[X32(k)]

|X 32(k ) |

k

x 32(n )

n

分析讨论:

(1) 时域采样定理的验证程序的运行结果如图实验2-1的波形图所示。由图可见,当采

样频率为1000Hz 时,频谱混叠很小;当采样频率为300Hz 时,频谱混叠很严重;当采样频率为200Hz 时,频谱混叠更严重。

(2) 频域采样定理的验证程序的运行结果如图实验2-2的波形图所示。该图验证了频域

采样定理和频域采样理论。对信号x(n)的频谱函数X (e ) 在[0,2π]上等间隔采

jw

样N=16时,N 点IDFT [X N (k)]得到的序列正是原序列x(n)以16为周期进行周期延拓后的主值区序列,由于N

思考题

先对原序列x (n ) 以N 为周期进行周期延拓后去主值区序列,

x N (n ) =[∑x (n +Ni ) ]R N (n )

i =-∞

再计算N 点的DFT ,则得到N 点的频域采样:

X N (k ) =DFT [x N (n )]N =X (e jw ) |

w =

2πk N

k =0, 1, 2,..., N -1

实验报告要求

(1) (2) (3) (4)

运行程序,打印要求显示的图形。

对各实验所得的结果进行简单的分析和解释。 简要回答思考题。

打印程序清单和要求的各信号波形。

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

实验目的

学习用FFT 对连续信号的时域离散信号进行谱分析的方法,了解可能出现的分析误差及其原因,以便正确应用FFT 。

实验原理

用FFT 对信号作频谱分析是学习数字信号处理的重要内容。经常需要进行谱分析的信号是模拟信号和时域离散信号。对信号进行谱分析的重要问题是频谱分辨率D 和分析误差。频谱分辨率直接和FFT 的变换区间N 有关,因为FFT 能够实现的频率分辨率是2∏/N,因此要求2∏/N≤D 。可以根据此式选择FFT 的变换区间N 。误差主要来自于用FFT 频谱分析时,得到的是离散谱,而信号(周期信号除外)是连续谱的,只有当N 较大时,离散谱的包络才能逼近于连续谱,因此N 要适当选择大一些。

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

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

析进行。

实验步骤及内容

(1) 对以下序列进行谱分析:

x 1(n ) =R 4(n )

⎧n +1(0≤n ≤3) ⎪

x 2(n ) =⎨8-n (4≤n ≤7)

⎪0(其他n ) ⎩⎧4-n (0≤n ≤3) ⎪

x 3(n ) =⎨n -3(4≤n ≤7)

⎪0(其他n ) ⎩

选择FFT 的变换区间N 为8和16两种情况进行频谱分析。分别打印其幅频特性曲线,并进行对比,分析和讨论。

(2) 对以下周期序列进行谱分析:

x 4(n ) =cos x 5(n ) =cos

π

4

n

πn

4

+cos

πn

8

选择FFT 的变换区间N 为8和16两种情况分别对以上序列进行频谱分析。分别打印其幅频特性曲线,并进行对比,分析和讨论。 (3) 对模拟周期信号进行谱分析:

x 6(t ) =cos 8πt +cos 16πt +cos 20πt

选择采样频率F s =64Hz ,对变换区间N=16,32,64三种情况进行谱分析。分别打印其幅频特性曲线,并进行分析和讨论。

实验程序

实验3—1程序: clear all; close all;

x1n=[ones(1,4)]; M=8;xa=1:(M/2); xb=(M/2):-1:1; x2n=[xa,xb]; x3n=[xb,xa];

X1k8=fft(x1n,8); X1k16=fft(x1n,16);

X2k8=fft(x2n,8); X2k16=fft(x2n,16); X3k8=fft(x3n,8); X3k16=fft(x3n,16); subplot(2,2,1); mstem(X1k8);

title('(1a)8点DFT[x_1(n)]'); xlabel('ω/π');

ylabel('幅度');

axis([0,2,0,1.2*max(abs(X1k8))]) subplot(2,2,2); mstem(X1k16);

title('(1b)16点DFT[x_1(n)]'); xlabel('ω/π');

ylabel('幅度');

axis([0,2,0,1.2*max(abs(X1k16))]) subplot(2,2,3); mstem(X2k8);

title('(2a)8点DFT[x_2(n)]'); xlabel('ω/π');

ylabel('幅度');

axis([0,2,0,1.2*max(abs(X2k8))]) subplot(2,2,4); mstem(X2k16);

title('(2b)16点DFT[x_2(n)]'); xlabel('ω/π');

ylabel('幅度'); figure(2)

axis([0,2,0,1.2*max(abs(X2k16))]) subplot(2,2,1); mstem(X3k8);

title('(3a) 8点DFT[x_3(n)]'); xlabel('ω/π');

ylabel('幅度');

axis([0,2,0,1.2*max(abs(X3k8))]) subplot(2,2,2); mstem(X3k16);

title('(3b)16点DFT[x_3(n)]'); xlabel('ω/π'); ylabel('幅度');

axis([0,2,0,1.2*max(abs(X3k16))])