华工-数字信号处理实验报告.doc
《华工-数字信号处理实验报告.doc》由会员分享，可在线阅读，更多相关《华工-数字信号处理实验报告.doc（48页珍藏版）》请在冰点文库上搜索。
数字信号处理第一次实验报告
数字信号处理实验报告
班级：
姓名：
组号：
第九组
日期：
二零一四年十一月
实验1常见离散信号产生和实现
一、实验目的
1、加深对常用离散信号的理解；
2、熟悉使用MATLAB在时域中产生一些基本的离散时间信号。
二、实验原理
1、单位抽样序列
在MATLAB中可以利用函数实现。
2、单位阶越序列
在MATLAB中可以利用函数实现:
3、正弦序列
在MATLAB中实现过程如下：
4、复指数序列
在MATLAB中实现过程如下：
5、指数序列
在MATLAB中实现过程如下：
三、预习要求
1、预先阅读实验讲义（MATLAB基础介绍）；
2、讨论正弦序列、复指数序列的性质。
A．绘出信号，当、时、、时的信号实部和虚部图；当时呢？
此时信号周期为多少？
程序dsp1.m如下：
titlez1=-1/12+j*pi/6;
titlez2=1/12+j*pi/6;
z3=1/12;
z4=2+j*pi/6;
z5=j*pi/6;
n=0:
20;
x1=exp（titlez1*n）;
x2=exp（titlez2*n）;
x3=exp（z3*n）;
x4=exp（z4*n）;
x5=exp（z5*n）;
subplot（5,2,1）;
stem（n,real（x1））;
xlabel（'n'）;ylabel（'real（x1）'）;
title（'z1=-1/12+j*pi/6时'）
subplot（5,2,2）;
stem（n,imag（x1））;
xlabel（'n'）;ylabel（'imag（x1）'）;
title（'z1=-1/12+j*pi/6时'）
subplot（5,2,3）;
stem（n,real（x2））;
xlabel（'n'）;ylabel（'real（x2）'）;
title（'z2=1/12+j*pi/6时'）
subplot（5,2,4）;
stem（n,imag（x2））;
xlabel（'n'）;ylabel（'image（x2）'）;
title（'z2=1/12+j*pi/6时'）
subplot（5,2,5）;
stem（n,real（x3））;
xlabel（'n'）;ylabel（'real（x3）'）;
title（'z3=1/12时'）
subplot（5,2,6）;
stem（n,imag（x3））;
xlabel（'n'）;ylabel（'image（x3）'）;
title（'z3=1/12时'）
subplot（5,2,7）;
stem（n,real（x4））;
xlabel（'n'）;ylabel（'real（x4）'）;
title（'z4=2+j*pi/6时'）
subplot（5,2,8）;
stem（n,imag（x4））;
xlabel（'n'）;ylabel（'image（x4）'）;
title（'z4=2+j*pi/6时'）
subplot（5,2,9）;
stem（n,real（x5））;
xlabel（'n'）;ylabel（'real（x5）'）;
title（'z5=j*pi/6时'）
subplot（5,2,10）;
stem（n,imag（x5））;
xlabel（'n'）;ylabel（'image（x5）'）;
title（'z5=j*pi/6时'）
运行结果如下：
结论：
当Z=pi/6时，序列周期为12。
B．绘出信号的频率是多少？
周期是多少？
产生一个数字频率为0.9的正弦序列，并显示该信号，说明其周期。
（1）当x1=1.5*sin（2*pi*0.1*n）时
程序dsp2.m如下：
n=-20:
20;
x=1.5*sin（2*pi*0.1*n）;
stem（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'x1=1.5*sin（2*pi*0.1*n）'）
程序运行结果如下：
由图可知，x1的周期为10，频率为0.1.
（2）当x2=1.5*sin（0.9*pi*0.1*n）时
程序dsp3.m如下：
n=-20:
20;
x2=1.5*sin（0.9*n）;
stem（n,x2）
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'x2=1.5*sin（0.9*pi*0.1*n）'）
程序运行结果如下：
由图可知，x2为非周期函数，由理论分析可得，0.9不是π的倍数，所以不是周期函数。
3、使用帮助功能学习square（方波）,sawtooth（锯齿波）和sinc函数，并绘图。
程序如下：
（1）方波程序dsp4.m：
%%%%%%%%%%%%%%%%%%%%
n=-2*pi:
0.0001:
2*pi
y=square（n）;
plot（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'square'）
方波结果：
%%%%%%%%%%%%%%%%%%%
（2）锯齿波程序dsp5.m：
n=-10:
0.0001:
10;
y=sawtooth（n）;
plot（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'sawtooth'）
锯齿波结果：
%%%%%%%%%%%%%%%%%%%%%
（3）sinc函数程序dsp6.m：
n=-10:
0.1:
10;
y=sinc（n）;
plot（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'sinc'）
Sinc函数结果：
四、实验内容
编制程序产生上述5种信号，长度可输入确定，函数需要的参数可输入确定，并绘出其图形。
（1）单位抽样序列
函数部分fun1.m：
functiony=fun1（n）
y=[zeros（1,n-1）1zeros（1,n-1）];
调用部分dsp7.m：
n=3;
y=fun1（n）;
stem（（-n+1）:
（n-1）,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
（2）阶跃序列
函数部分fun2.m：
functiony=fun2（n）
y=[zeros（1,n-1）ones（1,n）];
调用部分dsp8.m：
n=3;
y=fun2（n）;
stem（（-n+1）:
（n-1）,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
（3）正弦序列dsp9.m
n=-20:
20;
A=2;
fai=pi/6;
f=10;
fs=100;
y=A*sin（2*pi*f*n/fs+fai）;
clf;
stem（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
图形如下：
（4）复指数序列
函数部分fun4.m：
functiony=fun4（n,r,w）
y=r*exp（j*w*n）;
subplot（2,1,1）;
stem（n,real（y））
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'realpart'）
subplot（2,1,2）;
stem（n,imag（y））;
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'imaginarypart'）;
调用部分dsp10.m：
n=-20:
20;
r=2;
w=5;
y=fun4（n,r,w）;
图形如下：
（5）指数序列
函数部分fun5.m：
functiony=fun5（a,n）;
y=a.^n;
stem（n,y）
xlabel（'时间序号'）;
ylabel（'幅度'）;
xlabel（'时间序号'）;
ylabel（'幅度'）;
title（'指数序号'）
调用部分dsp11.m：
n=-5:
10;
a=2;
y=fun5（a,n）;
图形如下：
实验2离散系统的时域分析
一、实验目的
1、熟悉并掌握离散系统的差分方程表示法；
2、加深对冲激响应和卷积分析方法的理解。
二、实验原理
在时域中，离散时间系统对输入信号或者延迟信号进行运算处理，生成具有所需特性的输出信号，具体框图如下:
其输入、输出关系可用以下差分方程描述：
输入信号分解为冲激信号,
记系统单位冲激响应，则系统响应为如下的卷积计算式：
当时，h[n]是有限长度的（），称系统为FIR系统；反之，称系统为IIR系统。
三、预习要求
1、在MATLAB中，熟悉利用函数实现差分方程的仿真；
2、在MATLAB中，熟悉用函数计算卷积，用求系统冲激响应的过程。
四、实验内容
1、以下程序中分别使用conv和filter函数计算h和x的卷积y和y1，运行程序，并分析y和y1是否有差别，为什么要使用x[n]补零后的x1来产生y1；具体分析当h[n]有i个值，x[n]有j个值，使用filter完成卷积功能，需要如何补零？
%Programdsp12.m
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;
x1=[xzeros（1,8）];
y1=filter（h,1,x1）;
subplot（2,1,2）;
stem（n,y1）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）;
title（'OutputGeneratedbyFiltering'）;grid;
答：
由于filter的使用方法中，所得矩阵的长度为max{length（A），length（B）}-1，但真正卷积所得矩阵长度为A+B-1，所以在使用filter时,先算得两个要卷积的数组的最大长度A+B-1,然后若A卷积B,则用filter时要将B补（A-1）个0,若B卷积A,则对A补（B-1）个0.而conv函数可以直接补零。
2、编制程序求解下列两个系统的单位冲激响应和阶跃响应，并绘出其图形。
要求分别用filter、conv、impz三种函数完成。
，
给出理论计算结果和程序计算结果并讨论。
（1）
①系统的单位冲激响应：
A.使用Filter函dsp13.m：
a=[1,0.75,-.125];
b=[1,-1];
n=0:
20;
x=[1,zeros（1,20）];
yfilter=filter（b,a,x）;
stem（n,yfilter）;
title（'filter'）;
xlabel（'x'）;
ylabel（'y'）;
B.使用Conv函数dsp14.m:
a=[1,0.75,0.125];
b=[1,-1];
x=[1zeros（1,10）];
[h]=impz（b,a,10）;
yconv=conv（h,x）;
n=0:
19;
stem（n,yconv）;
title（'conv'）;
xlabel（'x'）;
ylabel（'y'）;
C.使用Impz函数dsp15.m:
a=[1,0.75,0.125];
b=[1,-1];
c=impz（b,a,21）;
stem（c）
②系统的阶跃响应：
A.使用Filter函数dsp16.m:
a1=[1,0.75,0.125];
b1=[1,-1];
n=0:
20;
x2=ones（1,21）;
y1filter=filter（b1,a1,x2）;
stem（n,y1filter）;
title（'阶跃响应'）;
xlabel（'x'）;
ylabel（'y'）;
B.使用Impz函数dsp17.m:
a=[1,0.75,0.125];
b=1;
impz（b,a）
（2）
①冲击响应：
A.使用filter函数dsp18.m：
a=[1];
b=[0,0.25,0.25,0.25,0.25];
n=0:
20;
x=[1,zeros（1,20）];
yfilter=filter（b,a,x）;
stem（n,yfilter）;
title（'filter'）;
xlabel（'x'）;
ylabel（'y'）
B.使用Impz函数dsp19.m:
a=[1];
b=[0,0.25,0.25,0.25,0.25];
c=impz（b,a,21）;
stem（c）
②阶跃响应：
A.使用Filter函数dsp20.m:
a1=[1];
b1=[0,0.25,0.25,0.25,0.25];
n=0:
20;
x2=ones（1,21）;
y1filter=filter（b1,a1,x2）;
stem（n,y1filter）;
title（'阶跃响应'）;
xlabel（'x'）;
ylabel（'y'）
B.使用Impz函数dsp21.m:
n=0:
20;
b=[0,0.25,0.5,0.75,ones（1,17）];
a=1;
impz（b,a,21）
若用y=filter（p,d,x）用来实现差分方程，d表示差分方程输出y的系数，p表示输入x的系数，而x表示输入序列。
输出结果函数的长度即为x的长度。
而y=conv（x,h）是用来实现卷积的，对x序列和h序列进行卷积，输出的结果个数等于x的长度与h的长度之和减去1。
y=impz（p,d,N）是用来实现冲击响应的，d表示输入x的系数，d表示输出y的长度，N是x的长度。
作业题
M2.7
程序M27.m如下：
t=0:
0.001:
1;
f0=cos（6*pi*t）; %三个信号频率分别为3、7、13Hz
f1=cos（14*pi*t）;
f2=cos（26*pi*t）;
plot（t,f0,t,f1,t,f2）;
hold
n=0:
1:
10;
fn=cos（0.6*pi*n）; %抽样函数频率为0.3Hz
plot（n/10,fn,'o'）;
holdoff
实验结果如下：
结论：
如图所示，蓝色虚线，绿色虚线和红色实线相交于抽样点（蓝色圆圈），所以对于这三个信号抽样结果是相同的。
M2.9
程序M29.m如下：
N=10;
n=0:
N-1;
y=cos（2*n）;
x=rand（1,N）+y-0.5; %产生随机函数rand从零到N来模拟信号x的噪声
r=conv（x,fliplr（x））; %自相关序列
n1=length（x）-1;
k=-n1:
n1
stem（k,r）; %画出自相关离散图像
xlabel（'Lagindex'）;ylabel（'Amplitude'）;
实验结果如下：
结论：
如图所示，当n=0的时候序列的自相关序列呈现尖峰。
M4.1
程序M41.m如下：
%Program4_1
%SignalSmoothingbyaMoving-AverageFilter
R=50;
d=rand（R,1）-0.5;
m=0:
1:
R-1;
s=2*m.*（0.9.^m）;
x=s+d';
plot（m,d,'r-',m,s,'b--',m,x,'g:
'）
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
legend（'d[n]','s[n]','x[n]'）;
pause
M=input（'Numberofinputsamples='）;
b=ones（M,1）/M;
y=filter（b,1,x）;
plot（m,s,'r-',m,y,'b--'）
legend（'s[n]','y[n]'）;
xlabel（'Timeindexn'）;ylabel（'Amplitude'）
实验结果如下：
M=5
M=7
M=9
结论：
如图，红色实线是平滑后的输出，蓝色虚线是含噪输出。
由图可知随着滤波器长度的增加，曲线的平滑度上升，但是平滑后输出与含噪输出之间的延迟也增加了。
M4.2
程序M42.m如下：
a=9; %设定x[n]的系数alpha
y0=1; %设定初始值y[-1]=-1
y1=0.5*（y0+（a/y0））; %得到递归方程
whileabs（y1-y0）>0.00001 %当相邻两个y值接近时（收敛）的时候
y2=0.5*（y1+（a/y1））; %停止运算与赋值
y0=y1;
y1=y2;
end
实验结果如下：
1）若a=9
此时，变量结果如图：
2）若a=100
此时，变量结果如图：
可证y[n]收敛于。
M4.5
程序M45.m如下：
h=[-4.87889.5631-4.8788];
[H,w]=freqz（h,1）; %求出冲激响应和频率范围
m=abs（H）; %求出频响幅值
plot（w/pi,m）;
grid；
实验结果如下：
PPT声音处理题
课前每人录制一段自己的拼音字母“a、o、e、b、p、t”、以及“a、o、e”四个声调的声音（采样率8kHz，单声道，16bit量化）
利用文献中的方法编程求解出各自声音信号的基音周期，观察不同字母、声调的基音变化
语音信号需分帧处理，20ms一帧（160个样点），每一帧求一个周期
程序：
举例：
a的一声a1.wav程序ayisheng.m如下：
a1=wavread（'a1.wav'）;
subplot（3,1,1）;
stem（a1（6000:
11200））; %查看波形
title（'声音波形'）;
beg=6000; %声音开始位置
n=20; %帧数
fl=160; %每帧长度
subplot（3,1,2）;
fori=1:
n
s=a1（（beg+fl*（i-1））:
（beg+fl*i））;%截取一帧
x=xcorr（s）;%计算自相关
stem（x）
title（'自相关'）;
[M,I1]=max（x）;%找最大值位置
[M,I2]=max（x（1:
I1-10））;%找最大值附近极大值位置
d（i）=I1-I2%计算周期
end
subplot（3,1,3）;
stem（d）
dm=mode（d） %求基音周期众数dm
fd=8000/dm %求声音频率fd
xlabel（'帧数n'）;ylabel（'d'）;
title（'基音周期'）
运行结果：
1.“a”
（1）一声时
d=[292929292929282828282829 2929292929292929];
基音周期众数：
dm=29
频率：
fd=275.8621
（2）二声时
d=[4140404079793978773838383838373737363636];
基音周期众数：
dm=38
频率fd=210.5263
（3）三声时
d=[3940414141424242434343444444454545454545]
基音周期众数：
dm=45
频率：
fd=177.7778
（4）四声时
d=[1010172121222222222222232323232324242424]
基音周期众数：
dm=22
频率：
fd=363.6364
（5）四个声调对比：
初步结论：
由对比图可得后三种音调的波形图不如a的一声平稳，因为二三四声都经过了变调。
而从基音周期看，四种声调基音周期和频率各不相同，其中，四声基音周期最小，频率最大。
2.清音和浊音的对比
（1）a音一声
d=[292929292929282828282829 2929292929292929];
基音周期众数：
dm=29
频率：
fd=275.8621
（2）o音一声：
d=[2627272727272727272726262626262626262626]
基音周期众数：
dm=26
频率：
fd=307.6923
升级会员 