常见离散信号产生和实现

数字信号处理实验报告

实验1 常见离散信号产生和实现

一、实验目的

1、加深对常用离散信号的理解；

2、熟悉使用MATLAB 在时域中产生一些基本的离散时间信号。

二、实验原理

MATLAB 语言提供了一系列函数用来产生信号，如exp, sin, cos, square, sawtooth ，ones, zeros 等函数。

1.

基本信号序列

1) 单位抽样序列

?

??=01)(n δ

00

≠=n n 在MATLAB 中可以利用zeros()函数实现。x=[1 zeros(1, n-1)]

示范程序: % Program P1_1

% Generation of a Unit Sample Sequence clf;

% Generate a vector from -10 to 20 n = -10:20;

% Generate the unit sample sequence u = [zeros(1,10) 1 zeros(1,20)]; % Plot the unit sample sequence stem(n,u);

xlabel('Time index n');ylabel('Amplitude'); title('Unit Sample Sequence'); axis([-10 20 0 1.2]);

如果)(n δ在时间轴上延迟了k 个单位，得到)(k n -δ即：

??

?=-01)(k n δ 0

≠=n k n 2) 单位阶跃序列

??

?0

1)(n u 00<≥n n 在MATLAB 中可以利用ones()函数实现。

);,1(N ones x =

3) 实指数序列

R a n a n x n

∈?=,)(

MATLAB 实现：

n

a x N n .^1:0=-=

4) 复指数序列

n e n x n

j ?=+)()(ωσ

MATLAB 实现：

)

)0exp((1

:0n w j lu x N n **+=-=

5) 随机序列

MATLAB 提供了两种随机信号：

rand(1,N) 产生[0,1]上均匀分布的随机矢量。

randn(1,N) 产生均值为0，方差为1的高斯随机序列，即白噪声序列。

2.

基本周期波形

1) 方波

MATLAB 工具箱函数square 可以产生方波;

t=0:0.1*pi:6*pi; y=square(t);

axis([0 7*pi -1.5 1.5]);

plot(t,y); xlabel(‘时间 t ’); ylabel(‘幅度y ’); 2) 正弦波

)/2sin()(?π+=Fs fn A n x

在MATLAB 中

)

/***2sin(*1

:0fai Fs n f pi A x N n +=-=

3) 锯齿波

工具箱函数sawtooth 函数可以产生锯齿波

Fs=10000;

t=0:1/Fs:1.5; %抽样长度1.5s, 抽样频率为10kHz x=sawtooth(2*pi*50*t); %信号频率为50Hz plot(t, x);

axis()[0 0.2 -1 1]; %画出0.2秒的波形

3.

基本非周期波形

工具箱函数chirp 能产生一种扫射频率信号，其特点是信号的瞬时频率随时间按照一定规律变化

t=0:1/1000:2 %抽样频率1kHz, 抽样时间2s.

x=chirp(t, 0.1, 150) %0时刻为DC 信号，1s 时频率为150Hz 。 specgram(x, 256, 1000, 256, 250); 4.

sinc 信号

MATLAB 实现： t=linspace(-5,5); x=sinc(t); plot(t,x); 5.

序列的操作

1) 信号加 x(n)=x1(n)+x2(n) MATLAB 实现：

x=x1+x2;

注意：x1和x2序列应该具有相同的长度，位置对应，才能相加。 2) 信号乘 x(n)=x1(n)*x2(n) MATLAB 实现：

x=x1.*x2; %数组乘法 3) 改变比例 y(n)=k*x(n) MATLAB 实现： y=k*x;

4) 折叠 y(n)=x(-n); MATLAB 实现： y=fliplr(x);

5) 抽样和 ∑==

2

1

)(n n n n x y

MATLAB 实现：

y=sum(x(n1:n2));

6) 抽样积 ∏==2

1

)(n n n n x y

MATLAB 实现：

y=prod(x(n1:n2));

7) 信号能量 2

|)(|∑∞

-∞

==

n n x Ex

MATLAB 实现：

Ex=sum(abs(x).^2);

8) 信号功率 2

1

|

)(|1∑-==N n n x N Px

MATLAB 实现：

Px=sum(abs(x).^2)/ N;

实验和程序： 1) 单位抽样序列

function uss(n)%构造函数

N=0:n-1

u = [1, zeros(1,n-1)];%构造矩阵

stem(N,u);

xlabel('Time index n');ylabel('Amplitude'); title('Unit Sample Sequence');

axis([-10 n 0 1.2]);

右图为取n为10的图形

2) 单位阶跃序列

function jieyue(n)%构造函数

N = -10:n-1;

u = [zeros(1,10) ones(1,n)];

stem(N,u);

xlabel('Time index

n');ylabel('Amplitude');

title('Unit jump Sequence');

axis([-10 n 0 1.2]);

右图为取n为10的图形

3) 实指数序列

function index(z1,N)%构造函数

% n1=0:N;

n=0:N/2;

x=z1.^n;

stem(n,x);

xlabel('Time index

n');ylabel('Amplitude');

右图为取z1=2,N=30的图形

4) 复指数序列

function findex(z1,z2)%构造函数

Time index n

A

m

p

l

i

t

u

d

e

Unit Sample Sequence

Time index n

A

m

p

l

i

t

u

d

e

Unit jump Sequence

4

Time index n A

m

p

l

i

t

u

d

e

5

Time index n

A

m

p

l

i

t

u

d

e

N=10;

n=0:N/2;

x=exp((z1+j*z2).*n);

stem(n,abs(x));

xlabel('Time index n');ylabel('Amplitude'); 右图为z1=2，z2=3的图形

5) 方波

function square1(g)%G为方波的占空比

t=0:0.1:6*pi;

y=square(t,g);

plot(t,y);

xlabel('时间 t');

ylabel('幅度y');

axis([0 7*pi -1.5 1.5]); %要放在最后才可以正确显示

右图分别为G=70，G=30的图形

6）正弦波

function sin1(N,A,Fs,f,Q)%N为显示

的范围，A为幅度，

%FS为抽样频率，f为模拟信号的频率，Q为初相位

n=0:N/2;

x=A*sin(2*pi*f*n/Fs+Q)

stem(n,x);

xlabel('Time index

n');ylabel('Amplitude');

axis([0 N/2 -2.5 2.5]);

右图为sin1(60,2,1000,50,10)的图形

7）锯齿波

时间 t

幅

度

y

时间 t

Time index n

A

m

p

l

i

t

u

d

e

function saw(f,F)%f 为锯齿波的频率，F 为间隔 %其中F 不能取太大，不易观 察图形 t=0:1/F:2;

y=sawtooth(2*pi*f*t); plot(t,y);

axis([0 2 -1.2 1.2]); 右图为saw(10,100)的图形 8）基本非周期波形

%利用chirp 来实现

t=0:1/500:2 %抽样频率500Hz, 抽样时间2s.

x=chirp(t, 0.1, 150) %0时刻为DC 信号，1s 时频

率为150Hz 。

specgram(x, 256, 500, 256, 250); 右图分别为抽样频率为500与10000的图形 9） sinc 信号

function sinc1(t1,t2,A)%[t1,t2]为定义域A

为幅度

t=linspace(t1,t2);

x=A*sinc(t);

plot(t,x);

右图为[t1,t2]=[-10，10]，A=2的图形

10）序列相加 function sequenceplus(x1,x2) n1=length(x1); n2=length(x2);

Time

F r e q u e n c

y

50100

150200250

Time

F r e q u e n c

y

50

100

150

200

250

-20

-10

10

20

-40

-2002040

n=-n1/2:(n1/2-1);

subplot(221)

stem(n,x1)

subplot(222)

stem(n,x2)

if(n1==n2) % x1，x2 的长度必须相同

x=x1+x2;

subplot(223)

stem(n,x);

axis([-n1 n1 0 2])

else error('x1与x2不一样长')

end

右图为x3=[ones(1,20),ones(1,10),ones(1,10)];

x4=[zeros(1,20),ones(1,10),zeros(1,10)];

sequenceplus(x3,x4)的图形

如果输入

x3=[ones(1,10),ones(1,10),ones(1,10)];

x4=[zeros(1,20),ones(1,10),zeros(1,10)];

则会显示x1与x2不一样长

11）序列相乘

function sequencenmultiply(x1,x2)

n1=length(x1);

n2=length(x2);

n=-n1/2:(n1/2-1);

subplot(221)

stem(n,x1)

subplot(222)

stem(n,x2)

if(n1==n2) % x1，x2的长度必须相-20-1001020

同

x=x1.*x2; subplot(223) stem(n,x);

axis([-n1 n1 0 2])

else disp('error x1与x2不一样长') end

右图为x3=[ones(1,20),ones(1,10),ones(1,10)];

x4=[zeros(1,20),ones(1,10),zeros(1,10)]; sequenceplus(x3,x4)的图形 如果输入

x3=[ones(1,10),ones(1,10),ones(1,10)]; x4=[zeros(1,20),ones(1,10),zeros(1,10)]; 则会显示x1与x2不一样长 12）改变比例

function changerate(x,k) n1=length(x); %x 为需要改变的序列 n=-n1/2:(n1/2-1); ，k 为改变的比例

x1=k.*x; stem(n,x1);

axis([-n1 n1 0 k+1]); 右图为

x=[ones(1,20),ones(1,10),ones(1,10)];K=1.5的图形 13)折叠

function pucker(x)

n1=length(x);

n=-n1/2:(n1/2-1);%这里的n 必须关于0两边对称 x1=fliplr(x); ，不然就会出错 subplot(211)

-40

-30

-20

-10

10

20

30

40

stem(n,x)%画出原序列

axis([-n1 n1 0 1.2]);

subplot(212)

stem(n,x1);%画出折叠后的序列

axis([-n1 n1 0 1.2]);

14)抽样和函数

n=0:24;

x1=[zeros(1,5),ones(1,3),zeros(1,5),ones(1,6),zeros(1,6)];

x=sum(x1)

计算结果为x =

9

15)序列积

n=0:25;

x1=rand(1,26)

x=prod(x1)

计算结果为x1 =

Columns 1 through 6

0.9501 0.2311 0.6068 0.4860 0.8913 0.7621

Columns 7 through 12

0.4565 0.0185 0.8214 0.4447 0.6154 0.7919

Columns 13 through 18

0.9218 0.7382 0.1763 0.4057 0.9355 0.9169

Columns 19 through 24

0.4103 0.8936 0.0579 0.3529 0.8132 0.0099

Columns 25 through 26

0.1389 0.2028

x =

4.6704e-012

16) 序列能量

x3=[ones(1,15),ones(1,20),ones(1,16)];

p=sum(abs(x3).^2)

计算结果为p =

51

17)信号功率

x1=[zeros(1,5),-2*ones(1,3),zeros(1,5),5*ones(1,6),zeros(1,6)];

N=length(x1)

p=sum(abs(x1).^2)/N

计算结果为N =

25

p =

6.4800

四、讨论复指数序列的性质。

绘出信号，当、时、、

时的信号实部和虚部图；当时呢？此时信号周期为多少？

调用复指数函数，绘出上面的z的图形如下所示

（1）

05101520253035404550

051015

202530354045

50

Time index n

A m p l i t u d e

（2）

5

101520253035404550

051015

2025303540

4550

Time index n

A m p l i t u d e

（3）

051015

20253035404550

Time index n

A m p l i t u d e

（4）

051015

x 10

11

05

1015

12

Time index n

A m p l i t u d e

（5）

Time index n

A m p l i t u d e

讨论：由上面的图形可以得出当z 的实部为负数时，经过复指数运算后的实部和虚部都是振幅按指数衰减的正弦振荡。而z 的实部为正数时，经过复指数运算后的实部和虚部都是振幅按 指数增长的正弦振荡。对于实数经过复指数运算后只有实部，而当z 为纯虚数时，经过复指数运

算后的实部与虚部都是正弦振荡的，且振荡周期为虚部除以PI 再2倍，上题中时，周

期为12。

实验2 离散系统的时域分析

一、实验目的

1、熟悉并掌握离散系统的差分方程表示法；

2、加深对冲激响应和卷积分析方法的理解。

二、实验原理

在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:

其输入、输出关系可用以下差分方程描述：

输入信号分解为冲激信号,

记系统单位冲激响应，则系统响应为如下的卷积计算式：

当

时，h[n]是有限长度的（），称系统为FIR 系统；反之，

称系统为IIR 系统

四、实验内容

1、用MATLAB 求系统响应

1） 卷积的实现

线性移不变系统可由它的单位脉冲响应来表征。若已知了单位脉冲响应和系统激励就可通过卷积运算来求取系统响应，即)(*)()(n h n x n y 程序：

function conv1(x,h)

y=conv(x,h); N=length(y); n=0:1:N-1;

disp('output sequence=');

disp(y);

stem(n,y);

xlabel('Time index n');ylabel('Amplitude');

输入为：

x=[-2 0 1 -1 3]

h=[1 2 0 -1]

function conv1(x,h)

输出为：output sequence=

-2 -4 1 3 1 5 1 -3

2）单位脉冲响应的求取

线性时不变因果系统可用MA TLAB的函数filter来仿真

y=filter(b,a,x);

其中，x和y是长度相等的两个矢量。矢量x表示激励，矢量a，b表示系统函数形式

滤波器的分母和分子系数，得到的响应为矢量y。例如计算以下系统的单位脉冲响应y(n)+0.7y(n-1)-0.45y(y-2)-0.6y(y-3)=0.8x(n)-0.44x(n-1)+0.36x(n-2)+0.02x(n-3) 程序实现：

function filter1(N,b,a) % a=input('Type in the vector a='); N=input('Desired impuse response

x=[1 zeros(1,N-1)]; length=');%b=input('Type in the vector b=');

y=filter(b,a,x);

disp('output sequence')

disp(y)

k=0:1:N-1;

stem(k,y);

xlabel('Time index n'); ylabel('Amplitude');

输入数据：

N=41

b=[0.8 -0.44 0.36 0.02]

a=[1 0.7 -0.45 -0.6]

输出为：output sequence

Columns 1 through 6

0.8000 -1.0000 1.4200 -0.9440 0.6998 -0.0627 Columns 7 through 12

-0.2076 0.5370 -0.5069 0.4719 -0.2363 0.0736 Columns 13 through 18

0.1253 -0.1964 0.2380 -0.1798 0.1151 -0.0187 Columns 19 through 24

-0.0430 0.0908 -0.0941 0.0809 -0.0445 0.0111 Columns 25 through 30

0.0207 -0.0362 0.0414 -0.0328 0.0198 -0.0038 Columns 31 through 36

-0.0081 0.0158 -0.0170 0.0142 -0.0081 0.0018 Columns 37 through 41

0.0036 -0.0065 0.0073 -0.0059 0.0035

Time index n

A m p l i t u d e

2、以下程序中分别使用conv和filter函数计算h和x的卷积y和y1，运行程序，并分析y和y1是否有差别，为什么要使用x[n]补零后的x1来产生y1；具体分析当h[n]有i个值，x[n]有j个值，使用filter完成卷积功能，需要如何补零？

clf;

h = [3 2 1 -2 1 0 -4 0 3]; %impulse response

x = [1 -2 3 -4 3 2 1]; %input sequence

y = conv(h,x);

n = 0:14;

subplot(2,1,1);

stem(n,y);

xlabel('Time index n'); ylabel('Amplitude');

title('Output Obtained by Convolution'); grid;

n1=length(h);

x1 = [x zeros(1,n1-1)] %要补零,因为输入x1与输出y1要等长的,

y1 = filter(h,1,x1)

subplot(2,1,2);

stem(n,y1);

xlabel('Time index n'); ylabel('Amplitude');

title('Output Generated by Filtering'); grid;

max(y-y1) %判断两种方法求出的有没有差别

以上程序段为有判断y与y1的差别的,用 conv 与filter求的结果一样,只是两个的原理不同.一个是卷积,一个是循环卷积.而对于循环卷积来说,它的两个对象必须长度相等,所以就需要进行补零

当h[n]有i个值，x[n]有j个值，使用filter完成卷积功能，需要补零,首先判断哪个序列的长度更大,再将短的一个序列后面补零,使其与长序列相等长度.

3、编制程序求解下列两个系统的单位冲激响应，分别用filter 和 impz实现，并绘出其图形。

给出理论计算结果和程序计算结果并讨论。

下面为程序段:

function zzz

%%%%%%%%%%%%%%%%%%%%(1)y[n]+0.75y[n-1]+0.125y[n-2]=x[n]-x[n-1]

%%%用filter求

t1=[1 0.75 0.125];%系统函数的分母

t2=[1 -1];%系统函数的分子

subplot(2,1,1)

filter(40,t2,t1)

title('filter')

%%%%用impz求

t1=[1 0.75 0.125];%系统函数的分母

t2=[1 -1];%系统函数的分子

subplot(212)

impz(20,t2,t1)

title('impz') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%(2)y[n]=0.5{x[n-1]+x[n-2]+x[n-3]+x[n-4]} %%%%%用filter求

t1=[1]%系统函数的分母

t2=[0.5 0.5 0.5 0.5]%系统函数的分子

subplot(211)

filter1(40,t2,t1)

title('filter')

%%%%用impz求

t1=[1]

t2=[0.5 0.5 0.5 0.5]

subplot(212)

filter1(40,t2,t1)

title('impz')

function filter1(N,b,a)% a=input('Type in the vector a=');