实验一、数字信号处理基础.doc
《实验一、数字信号处理基础.doc》由会员分享,可在线阅读,更多相关《实验一、数字信号处理基础.doc(14页珍藏版)》请在冰点文库上搜索。
实验一:
数字信号处理基础
一、实验目的
1、熟悉并掌握离散系统的差分方程表示法;
2、加深对冲激响应和卷积分析方法的理解;
3、加深对离散信号的DFT的理解;
4、熟悉离散系统的频率响应分析方法、加深对零、极点分布的概念理解;
5、掌握Matlab常用函数的使用方法。
二、实验原理
1.LTI系统的表示
在时域中,离散时间系统对输入信号或者延迟信号进行运算处理,生成具有所需特性的输出信号,具体框图如下:
其输入、输出关系可用以下差分方程描述:
输入信号分解为冲激信号,
记系统单位冲激响应,则系统响应为如下的卷积计算式:
当时,h[n]是有限长度的(n=1:
M),称系统为FIR系统;反之,称系统为IIR系统。
2.离散傅里叶变换(DFT)的定义
N点序列的DFT和IDFT变换定义式如下:
,
利用旋转因子具有周期性,可以得到快速算法(FFT)。
3.离散时间系统的变换域分析方法
离散系统的时域方程为
其变换域分析方法如下:
频域:
系统的频率响应为:
Z域:
系统的转移函数为:
分解因式:
,
其中和称为零、极点。
三、预习要求
1、在MATLAB中,熟悉利用函数y=filter(b,a,x)实现差分方程的使用方法;
2、在MATLAB中,熟悉用函数y=conv(x,h)计算卷积的使用方法;
3、在MATLAB中,熟悉用函数y=impz(b,a)求系统冲激响应的使用方法。
4、在MATLAB中,熟悉用函数y=stepz(b,a)求系统冲激响应的使用方法。
5、在MATLAB中,函数和计算N点序列的DFT正、反变换的使用方法。
6、在MATLAB中,熟悉函数tf2zp、zplane、freqz、residuez、zp2sos的使用。
函数的使用方法:
1、filter:
这是一个一维数字滤波器;也可以用作卷积运算,结果与conv函数一样。
y=filter(A,B,X),x为待处理序列,H=B/A为一个系统函数,y的结果是H与X相乘后作傅氏反变换的结果,如果A=1,则y为两个序列的卷积。
2、conv:
求卷积函数,y=conv(x,h),y是x与h卷积的结果。
3、impz:
求系统的冲激响应。
[h,t]=impz(b,a):
b、a分别为系统的传递函数的分子和分母系数向量。
返回系统(b,a)的冲激响应h和相应的时间轴向量。
4、stepz:
求系统的单位阶跃响应。
[s,n]=stepz(b,a):
求解离散系统的单位阶跃响应,其中b、a为向量,n为时间序列
[s,n]=stepz(b,a,N):
求解离散系统的单位阶跃响应,N为采样点数
stepz(b,a):
在当前窗口,用stem(n,s)绘出图形。
5、求离散时间序列的傅立叶变换和反变换:
X=fft(x,N),x=ifft(X,N)
其中,X为x的DFT变换,而x则是X的傅立叶反变换,N为计算的点数,可以能过上面两个函数实现DFT。
6、tf2ss:
由传递函数到状态空间的转换。
[A,B,C,D]=tf2ss(NUM,DEN),NUM,DEM为向量
H(s)=NUM(s)/DEN(s)到x=Ax+Bu,y=Cx+Du
7、zplane(b,a):
绘制由行向量b和a构成的系统函数的零极点分布图;
zplane(z,p):
绘制由列向量z确定的零点、列向量p确定的极点构成的零极点分布图。
8、freqz:
求离散时间系统的频率响应函数;
[h,w]=freqz(b,a,n)。
可以得到数字滤波器的n点复频响应值,这n个点均匀地分布这[0,pi]上,并将这n个频点的频率记录在w中,相应的频响值记录在h中,n缺省时取512点。
[h,f]=freqz(b,a,n,Fs),用于对H(jw)在[0,Fs/2]上等间隔采样n点,采样点频率及相应频响值分别记录在f和h中,由用户指定Fs(以Hz为单位)的值。
9、[r,p,c]=residuez(b,a):
把b(z)/a(z)展开成部分分式;
[b,a]=residuez(r,p,c):
根据部分分式的r、p、c数组,返回有理多项式。
10、zp2sos:
变系统零极点形式为二阶分割形式
[SOS,G]=ZP2SOS(Z,P,K);
SOS为矩阵形式,G为增益。
………………………………………………………………………………………………..
四、实验内容
1.Matlab函数conv和filter的使用
1.1实验要求
以下程序中分别使用conv和filter函数计算h和x的卷积y和y1,运行程序,并分析y和y1是否有差别,为什么要使用x[n]补零后的x1来产生y1;具体分析当h[n]有M个值,x[n]有M个值,使用filter完成卷积功能,需要如何补零?
%ProgramP2_7
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;
x1=[xzeros(1,8)];
y1=filter(h,1,x1);
subplot(2,1,2);
stem(n,y1);
xlabel('Timeindexn');ylabel('Amplitude');
title('OutputGeneratedbyFiltering');grid;
1.2实验结果与分析
上面程序的运行结果如下:
将程序作如下改动:
h=[321-210-403];%impulseresponse
x=[1-23-4321];%inputsequence
y=conv(h,x)
y1=filter(h,1,x)
y2=filter(x,1,h)
结果输出如下:
y=
Columns1through10
3-46-1090172-910
Columns11through15
-2-20563
y1=
3-46-109017
y2=
3-46-1090172-9
从结果分析可知,conv(h,x)是专门用来做卷积运算的,得出的结果长度为(nh+nx-1),而对于filter(h,1,x),得出的结果长度与第三个参数x相同,而对于前面的运算结果与conv的结果一致,所以,用filter函数进行卷积运算时要对x进行补0操作。
当h[n]有M个值,x[n]有M个值时,需要令x=[xzeros(M-1,1)]。
………………………………………………………………………………………………..
2.Matlab函数impz和stepz的使用
2.1实验要求
编制程序用impz和stepz分别求解下列两个系统的单位冲激响应和阶跃响应:
y[n]+0.75y[n-1]+0.125y[n-2]=x[n]–x[n-1]
y[n]=0.25*(x[n-1]+x[n-2]+x[n-3]+x[n-4])
给出理论计算结果和程序计算结果并讨论。
2.2程序代码
N=16;
n=0:
N-1;
y1=[1,.75,0.125];
x1=[1,-1];
y2=1;
x2=[0.25,0.25,0.25,0.25];
h1=6*(-0.5).^n.*(n>=0)-5*(-0.25).^n.*(n>=0);%式1的单位冲激响应理论计算值
h2=2*(-0.5).^n.*(n>=0)-(-0.25).^n.*(n>=0);%式1的单位阶跃响应理论计算值
h3=[0,0.25,0.25,0.25,0.25,zeros(1,11)];
h4=0.25*((n>=1)+(n>=2)+(n>=3)+(n>=4));
subplot(4,2,1);
impz(x1,y1,N);%求系统的冲激响应
title('式1的冲激响应matlab结果');
subplot(4,2,2);
stem(n,h1);
title('式1的冲激响应理论计算结果');
xlabel('n(samples)');
ylabel('Amplitude');
subplot(4,2,3);
stepz(x1,y1,N);%求系统的阶跃响应
title('式1的阶跃响应matlab结果');
subplot(4,2,4);
stem(n,h2);
title('式1的阶跃响应理论计算结果');
xlabel('n(samples)');
ylabel('Amplitude');
subplot(4,2,5);
impz(x2,y2,N);
title('式2的冲激响应matlab结果');
subplot(4,2,6);
stem(n,h3);
title('式2的冲激响应理论计算结果');
xlabel('n(samples)');
ylabel('Amplitude');
subplot(4,2,7);
stepz(x2,y2,N);
title('式2的阶跃响应matlab结果');
subplot(4,2,8);
stem(n,h4);
title('式2的阶跃响应理论计算结果');
xlabel('n(samples)');
ylabel('Amplitude');
………………………………………………………………………………………………..
2.3实验结果与分析
理论计算:
对式1:
y[n]+0.75y[n-1]+0.125y[n-2]=x[n]–x[n-1]
单位冲击响应为:
阶跃响应为:
对式2:
y[n]=0.25*(x[n-1]+x[n-2]+x[n-3]+x[n-4])
单位冲击响应为:
阶跃响应为:
由结果图像可得,通过impz和stepz计算的单位冲激响应和阶跃响应与理论计算结果基本一致。
………………………………………………………………………………………………..
3.用N点DFT计算2N点实数的DFT
3.1实验要求
2N点实数序列
N=64。
用一个64点的复数FFT程序,一次算出,并绘出的图形。
3.2实现原理
v[n]是长度为2N的实序列,V[k]表示该实序列的2N点DFT。
定义两个长度均为N的实序列g[n]和h[n]为g[n]=v[2n],h[n]=v[2n+1],。
G[k]和H[k]表示它们的N点DFT。
通过下式可以算出v[n]的2N点傅立叶变换。
在matlab中,可以运用函数X=fft(x,N)直接计算出x[n]的N点傅立叶变换。
………………………………………………………………………………………………..
3.3实现代码
N=64;
n=0:
N-1;
k=0:
2*N-1;
x=cos(14*pi*k/N)+0.5*cos(38*pi*k/N);
g=cos(14*pi*2*n/N)+0.5*cos(38*pi*2*n/N);%g[n]=v[2n]
h=cos(14*pi*(2*n+1)/N)+0.5*cos(38*pi*(2*n+1)/N);%h[n]=v[2n+1]
G=fft(g,N);
H=fft(h,N);
form=1:
N%利用G,K求X,2N点的傅氏变换
X(m)=G(m)+exp(-j*2*pi*(m-1)/(2*N))*H(m);
end
form=(N+1):
2*N
X(m)=G(m-N)+exp(-j*2*pi*(m-1)/(2*N))*H(m-N);
end
X1=fft(x,2*N);%直接计算2N点傅区变换
dX=X1-X;%作差看两种算法的区别
subplot(3,1,1);
stem(k,abs(X));
title('利用N点DFT求2N点DFT');
xlabel('k');
ylabel('|X|');
subplot(3,1,2);
stem(k,abs(X1));
title('直接求2N点DFT');
xlabel('k');
ylabel('|X1|');
subplot(3,1,3);
stem(k,abs(dX));
title('两种算法的差值');
xlabel('k');
ylabel('|cX|');
………………………………………………………………………………………………..
3.4实验结果与分析
从结果图像可以看出,运用N点DFT求2N点DFT的方法求得的结果与直接求2N点DFT的结果基乎一样,虽然其中存在一定的差别,但是从工程上讲,这种差别可以忽略不计。
………………………………………………………………………………………………..
4.用DFT验证Z变换
4.1实验要求
已知某序列x(n)的Z变换在单位圆上的N=64等分样点为:
。
用N点IFFT程序计算出和。
4.2实现原理
在matlab中,可以运用函数x=ifft(X,N)来计算频域函数X的N点傅立叶反变换IDFT。
而对于理论计算,可以运用傅立叶变换对进行变换:
由于DFT相当于无限长序列的傅立叶变换的一个抽样,所以对于理论值的DFT只要取DTFT的有限个点即可。
………………………………………………………………………………………………..
4.3实现代码
N=64;
k=0:
N-1;
X=1./(1-0.8*exp(-j*2*pi*k/N));%原系统函数
x1=ifft(X,N);%运用matlab求解傅立叶反变换
x2=0.8.^k;%理论计算结果
dx=x2-x1;%理论结果与数值计算结果作差对比
subplot(3,1,1);
stem(k,x1);
axis([07001.01]);%对图像的x轴修正
title('matlab计算结果');
xlabel('n');
ylabel('x1');
subplot(3,1,2);
stem(k,x2);
title('理论计算结果');
xlabel('n');
ylabel('x2');
subplot(3,1,3);
stem(k,dx);
title('matlab计算结果与理论结果的差值');
xlabel('n');
ylabel('dx');
………………………………………………………………………………………………..
4.4实验结果与分析
理论计算结果为:
。
从结果图中可以看出,matlab计算结果与理论计算结果基本一致,但是,从误差分析来看,matlab计算结果比理论计算结果要大一些,不过,这个误差非常小,所以,从工程上说,这种误差应该是可以忽略不计的。
………………………………………………………………………………………………..
5离散时间系统的变换域分析
5.1实验要求
用matlab函数,求系统
的零、极点和幅度频率响应和相位响应。
5.2实现代码
b=[0.05280.07970.12950.12950.7970.0528];
a=[1-1.81072.4947-1.88010.9537-0.2336];
[z,p]=tf2zp(b,a);%求系统的零、极点
disp('零点');disp(z);%列出系统的零点计算结果
disp('极点');disp(p);%列出系统的极点计算结果
[h,w]=freqz(b,a);%求系统函数
t=w/pi;%将w对pi进行归一化
subplot(3,1,1);
zplane(b,a);
title('系统的零极点图');
subplot(3,1,2);
plot(t,abs(h));
title('系统的幅频特性曲线');
xlabel('\omega/\pi');
ylabel('Magnitude');
subplot(3,1,3);
plot(t,angle(h));
title('系统的相频特性曲线');
xlabel('\omega/\pi');
ylabel('Phase');
………………………………………………………………………………………………..
5.3实验结果与分析
零点
-1.5870+1.4470i
-1.5870-1.4470i
0.8657+1.5779i
0.8657-1.5779i
-0.0669
极点
0.2788+0.8973i
0.2788-0.8973i
0.3811+0.6274i
0.3811-0.6274i
0.4910
从实验结果可以看出,通过函数tf2zp计算的系统零、极点与用函数zplane画出的系统的零、极点图一到。
通过函数freqz可以得出系统函数h,并可以求出系统幅频特性曲线和相频特性曲线。
………………………………………………………………………………………………..
13