常见离散信号产生和实现.docx
《常见离散信号产生和实现.docx》由会员分享，可在线阅读，更多相关《常见离散信号产生和实现.docx（46页珍藏版）》请在冰点文库上搜索。
常见离散信号产生和实现
数字信号处理实验报告
实验1常见离散信号产生和实现
一、实验目的
1、加深对常用离散信号的理解；
2、熟悉使用MATLAB在时域中产生一些基本的离散时间信号。
二、实验原理
MATLAB语言提供了一系列函数用来产生信号，如exp,sin,cos,square,sawtooth，ones,zeros等函数。
1.基本信号序列
1）单位抽样序列
在MATLAB中可以利用zeros（）函数实现。
x=[1zeros（1,n-1）]
示范程序:
%ProgramP1_1
%GenerationofaUnitSampleSequence
clf;
%Generateavectorfrom-10to20
n=-10:
20;
%Generatetheunitsamplesequence
u=[zeros（1,10）1zeros（1,20）];
%Plottheunitsamplesequence
stem（n,u）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'UnitSampleSequence'）;
axis（[-102001.2]）;
如果
在时间轴上延迟了k个单位，得到
即：
2）单位阶跃序列
在MATLAB中可以利用ones（）函数实现。
3）实指数序列
MATLAB实现：
4）复指数序列
MATLAB实现：
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（[07*pi-1.51.5]）;
plot（t,y）;
xlabel（‘时间t’）;
ylabel（‘幅度y’）;
2）正弦波
在MATLAB中
3）锯齿波
工具箱函数sawtooth函数可以产生锯齿波
Fs=10000;
t=0:
1/Fs:
1.5;%抽样长度1.5s,抽样频率为10kHz
x=sawtooth（2*pi*50*t）;%信号频率为50Hz
plot（t,x）;
axis（）[00.2-11];%画出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）抽样和
MATLAB实现：
y=sum（x（n1:
n2））;
6）抽样积
MATLAB实现：
y=prod（x（n1:
n2））;
7）信号能量
MATLAB实现：
Ex=sum（abs（x）.^2）;
8）信号功率
MATLAB实现：
Px=sum（abs（x）.^2）/N;
实验和程序：
1）单位抽样序列
functionuss（n）%构造函数
N=0:
n-1
u=[1,zeros（1,n-1）];%构造矩阵
stem（N,u）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'UnitSampleSequence'）;
axis（[-10n01.2]）;
右图为取n为10的图形
2）单位阶跃序列
functionjieyue（n）%构造函数
N=-10:
n-1;
u=[zeros（1,10）ones（1,n）];
stem（N,u）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'UnitjumpSequence'）;
axis（[-10n01.2]）;
右图为取n为10的图形
3）实指数序列
functionindex（z1,N）%构造函数
%n1=0:
N;
n=0:
N/2;
x=z1.^n;
stem（n,x）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
右图为取z1=2,N=30的图形
4）复指数序列
functionfindex（z1,z2）%构造函数
N=10;
n=0:
N/2;
x=exp（（z1+j*z2）.*n）;
stem（n,abs（x））;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
右图为z1=2，z2=3的图形
5）方波
functionsquare1（g）%G为方波的占空比
t=0:
0.1:
6*pi;
y=square（t,g）;
plot（t,y）;
xlabel（'时间t'）;
ylabel（'幅度y'）;
axis（[07*pi-1.51.5]）;%要放在最后才可以正确显示
右图分别为G=70，G=30的图形
6）正弦波
functionsin1（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（'Timeindexn'）;ylabel（'Amplitude'）;
axis（[0N/2-2.52.5]）;
右图为sin1（60,2,1000,50,10）的图形
7）锯齿波
functionsaw（f,F）%f为锯齿波的频率，F为间隔%其中F不能取太大，不易观察图形
t=0:
1/F:
2;
y=sawtooth（2*pi*f*t）;
plot（t,y）;
axis（[02-1.21.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信号
functionsinc1（t1,t2,A）%[t1,t2]为定义域A为幅度
t=linspace（t1,t2）;
x=A*sinc（t）;
plot（t,x）;
右图为[t1,t2]=[-10，10]，A=2的图形
10）序列相加
functionsequenceplus（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的长度必须相同
x=x1+x2;
subplot（223）
stem（n,x）;
axis（[-n1n102]）
elseerror（'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）序列相乘
functionsequencenmultiply（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的长度必须相同
x=x1.*x2;
subplot（223）
stem（n,x）;
axis（[-n1n102]）
elsedisp（'errorx1与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）改变比例
functionchangerate（x,k）
n1=length（x）;%x为需要改变的序列
n=-n1/2:
（n1/2-1）;，k为改变的比例
x1=k.*x;
stem（n,x1）;
axis（[-n1n10k+1]）;
右图为x=[ones（1,20）,ones（1,10）,ones（1,10）];K=1.5的图形
13）折叠
functionpucker（x）
n1=length（x）;
n=-n1/2:
（n1/2-1）;%这里的n必须关于0两边对称
x1=fliplr（x）;，不然就会出错
subplot（211）
stem（n,x）%画出原序列
axis（[-n1n101.2]）;
subplot（212）
stem（n,x1）;%画出折叠后的序列
axis（[-n1n101.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=
Columns1through6
0.95010.23110.60680.48600.89130.7621
Columns7through12
0.45650.01850.82140.44470.61540.7919
Columns13through18
0.92180.73820.17630.40570.93550.9169
Columns19through24
0.41030.89360.05790.35290.81320.0099
Columns25through26
0.13890.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）
（2）
（3）
（4）
（5）
讨论：
由上面的图形可以得出当z的实部为负数时，经过复指数运算后的实部和虚部都是振幅按指数衰减的正弦振荡。
而z的实部为正数时，经过复指数运算后的实部和虚部都是振幅按
指数增长的正弦振荡。
对于实数经过复指数运算后只有实部，而当z为纯虚数时，经过复指数运算后的实部与虚部都是正弦振荡的，且振荡周期为虚部除以PI再2倍，上题中
时，周期为12。
实验2离散系统的时域分析
一、实验目的
1、熟悉并掌握离散系统的差分方程表示法；
2、加深对冲激响应和卷积分析方法的理解。
二、实验原理
在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:
其输入、输出关系可用以下差分方程描述：
输入信号分解为冲激信号,
记系统单位冲激响应
，则系统响应为如下的卷积计算式：
当
时，h[n]是有限长度的（
），称系统为FIR系统；反之，称系统为IIR系统
四、实验内容
1、用MATLAB求系统响应
1）卷积的实现
线性移不变系统可由它的单位脉冲响应来表征。
若已知了单位脉冲响应和系统激励就可通过卷积运算来求取系统响应，即
程序：
functionconv1（x,h）
y=conv（x,h）;
N=length（y）;
n=0:
1:
N-1;
disp（'outputsequence='）;
disp（y）;
stem（n,y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
输入为：
x=[-201-13]
h=[120-1]
functionconv1（x,h）
输出为：
outputsequence=
-2-413151-3
2）单位脉冲响应的求取
线性时不变因果系统可用MATLAB的函数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）
程序实现：
functionfilter1（N,b,a）%a=input（'Typeinthevectora='）;N=input（'Desiredimpuseresponse
x=[1zeros（1,N-1）];length='）;%b=input（'Typeinthevectorb='）;
y=filter（b,a,x）;
disp（'outputsequence'）
disp（y）
k=0:
1:
N-1;
stem（k,y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
输入数据：
N=41
b=[0.8-0.440.360.02]
a=[10.7-0.45-0.6]
输出为：
outputsequence
Columns1through6
0.8000-1.00001.4200-0.94400.6998-0.0627
Columns7through12
-0.20760.5370-0.50690.4719-0.23630.0736
Columns13through18
0.1253-0.19640.2380-0.17980.1151-0.0187
Columns19through24
-0.04300.0908-0.09410.0809-0.04450.0111
Columns25through30
0.0207-0.03620.0414-0.03280.0198-0.0038
Columns31through36
-0.00810.0158-0.01700.0142-0.00810.0018
Columns37through41
0.0036-0.00650.0073-0.00590.0035
2、以下程序中分别使用conv和filter函数计算h和x的卷积y和y1，运行程序，并分析y和y1是否有差别，为什么要使用x[n]补零后的x1来产生y1；具体分析当h[n]有i个值，x[n]有j个值，使用filter完成卷积功能，需要如何补零？
clf;
h=[321-210-403];%impulseresponse
x=[1-23-4321];%inputsequence
y=conv（h,x）;
n=0:
14;
subplot（2,1,1）;
stem（n,y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'OutputObtainedbyConvolution'）;grid;
n1=length（h）;
x1=[xzeros（1,n1-1）]%要补零,因为输入x1与输出y1要等长的,
y1=filter（h,1,x1）
subplot（2,1,2）;
stem（n,y1）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'OutputGeneratedbyFiltering'）;grid;
max（y-y1）%判断两种方法求出的有没有差别
以上程序段为有判断y与y1的差别的,用conv与filter求的结果一样,只是两个的原理不同.一个是卷积,一个是循环卷积.而对于循环卷积来说,它的两个对象必须长度相等,所以就需要进行补零
当h[n]有i个值，x[n]有j个值，使用filter完成卷积功能，需要补零,首先判断哪个序列的长度更大,再将短的一个序列后面补零,使其与长序列相等长度.
3、编制程序求解下列两个系统的单位冲激响应，分别用filter和impz实现，并绘出其图形。
给出理论计算结果和程序计算结果并讨论。
下面为程序段:
functionzzz
%%%%%%%%%%%%%%%%%%%%
（1）y[n]+0.75y[n-1]+0.125y[n-2]=x[n]-x[n-1]
%%%用filter求
t1=[10.750.125];%系统函数的分母
t2=[1-1];%系统函数的分子
subplot（2,1,1）
filter（40,t2,t1）
title（'filter'）
%%%%用impz求
t1=[10.750.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.50.50.50.5]%系统函数的分子
subplot（211）
filter1（40,t2,t1）
title（'filter'）
%%%%用impz求
t1=[1]
t2=[0.50.50.50.5]
subplot（212）
filter1（40,t2,t1）
title（'impz'）
functionfilter1（N,b,a）%a=input（'Typeinthevectora='）;N=input（'Desiredimpuseresponselength='）;%b=input（'Typeinthevectorb='）;
x=[1zeros（1,N-1）];
y=filter（b,a,x）;
disp（'outputsequence'）;
disp（y）;
k=0:
1:
N-1;
stem（k,y）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
%%%%%%%%%%%%%%%%%%%%%%%
functionimpz1（N,b,a）%a=input（'Typeinthevectora='）;N=input（'Desiredimpuseresponselength='）;%b=input（'Typeinthevectorb='）;
%x=[1zeros（1,N-1）];
y1=impz（b,a,N）;
disp（'outputsequence'）
disp（y1）;
k=0:
1:
N-1;
stem（k,y1）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
第一,二题的图形如下图所示:
理论计算的结果为
（1）h（n）=-5*（-0.25）^n*u（n）+6*（-0.5）^n*u（n）
（2）h（n）=0.25*[d（n-1）+d（n-2）+d（n-3）+d（n-4）]d代表单位冲激响应
下面是用理论计算的结果会出的第一题的图形,与程序的结果大致相同.
分析以上结果可以知道,用两种函数都可以求出差分方程的单位冲激响应.所求结果也大致相同.
实验3离散系统的变换域分析
一、实验目的
1、熟悉对离散系统的频率响应分析方法；
2、加深对零、极点分布的概念理解。
二、实验原理
离散系统的时域方程为
其变换域分析方法如下：
频域:
系统的频率响应为:
Z域:
系统函数为:
分解因式:
，
其中
和
称为零、极点。
四、实验内容
求系统
的零、极点和幅度频率响应和相位响应。
程序段如下;
num=[0.05280.07970.12950.12950.7970.05
升级会员 