离散时间系统及离散卷积.docx
《离散时间系统及离散卷积.docx》由会员分享，可在线阅读，更多相关《离散时间系统及离散卷积.docx（21页珍藏版）》请在冰豆网上搜索。
离散时间系统及离散卷积
实验一、离散时间系统及离散卷积
1、单位脉冲响应
源程序：
functionpr1（）%定义函数pr1
a=[1,-1,0.9];%定义差分方程y（n）-y（n-1）+0.9y（n-2）=x（n）
b=1;
x=impseq（0,-20,120）;%调用impseq函数（matlab软件的函数库）
n=[-20:
120];%定义n的范围，从-20到120
h=filter（b,a,x）;%调用函数给纵坐标赋值
figure
（1）%绘图figure1（冲激响应）
stem（n,h）;%在图中绘出冲激
title（'单位冲激响应（耿海锋）'）;%定义标题为:
'冲激响应（耿海锋）'
xlabel（'n'）;%绘图横座标为n
ylabel（'h（n）'）;%绘图纵座标为h（n）
figure
（2）%绘图figure2
[z,p,g]=tf2zp（b,a）;%绘出零极点图
zplane（z,p）
function[x,n]=impseq（n0,n1,n2）%声明impseq函数
n=[n1:
n2];
x=[（n-n0）==0];
结果：
Figure1:
Figure2:
2、离散系统的幅频、相频的分析
源程序：
functionpr2（）
b=[0.0181,0.0543,0.0543,0.0181];
a=[1.000,-1.76,1.1829,-0.2781];
m=0:
length（b）-1;%m的范围，从0到3
l=0:
length（a）-1;%l的范围，从0到3
K=5000;
k=1:
K;
w=pi*k/K;%角频率w
H=（b*exp（-j*m'*w））./（a*exp（-j*l'*w））;%对系统函数的定义
figure
（1）
magH=abs（H）;%magH为幅度
angH=angle（H）;%angH为相位
plot（w/pi,magH-耿海锋）;%绘制w（pi）-magH-耿海锋的图形
figure
（2）
axis（[0,1,0,1]）;%限制横纵座标从0到1
xlabel（'w（pi）'）;%x座标为w（pi）
ylabel（'|H|'）;%y座标为angle（H）-耿海锋
title（'幅度，相位响应（耿海锋）'）;%图的标题为:
'幅度，相位响应（耿海锋）'
plot（w/pi,angH）;%绘制w（pi）-angH的图形
grid;%为座标添加名称
xlabel（'w（pi）'）;%x座标为w（pi）
ylabel（'angle（H）'）;%y座标为angle（H）
结果：
Figure1
Figure2
3、卷积计算
源程序：
functionpr3（）
n=-5:
50;%声明n的范围，从-5到50
u1=stepseq（0,-5,50）;%调用stepseq函数声用明u1=u（n）
u2=stepseq（10,-5,50）;%调用stepseq函数声用明u2=u（n-10）
%输入x（n）和冲激响应h（n）
x=u1-u2;%x（n）=u（n）-u（n-10）
h=（（0.9）.^n）.*u1;%h（n）=0.9^n*u（n）
figure
（1）
subplot（3,1,1）;%绘制第一个子图
stem（n,x）;%绘制图中的冲激
axis（[-5,50,0,2]）;%限定横纵座标的范围
title（'输入序列-52101702-耿海锋'）;%规定标题为:
'输入序列-52101702-耿海锋'
xlabel（'n'）;%横轴为n
ylabel（'x（n）'）;%纵轴为x（n）
subplot（3,1,2）;%绘制第二个子图
stem（n,h）;%绘制图中的冲激
axis（[-5,50,0,2]）;%限定横纵座标的范围
title（'冲激响应序列-52101702-耿海锋'）;%规定标题为:
'冲激响应序列-52101702-耿海锋'
xlabel（'n'）;%横轴为n
ylabel（'h（n）'）;%纵轴为h（n）
%输出响应
[y,ny]=conv_m（x,n,h,n）;%调用conv_m函数
subplot（3,1,3）;%绘制第三个子图
stem（ny,y）;
axis（[-5,50,0,8]）;
title（'输出响应-52101702-耿海锋'）;%规定标题为:
'输出响应-52101702-耿海锋'
xlabel（'n'）;
ylabel（'y（n）'）;%纵轴为y（n）
%stepseq.m子程序
%实现当n>=n0时x（n）的值为1
function[x,n]=stepseq（n0,n1,n2）
n=n1:
n2;
x=[（n-n0）>=0];
%con_m的子程序
%实现卷积的计算
function[y,ny]=conv_m（x,nx,h,nh）
nyb=nx
（1）+nh
（1）;
nye=nx（length（x））+nh（length（h））;
ny=[nyb:
nye];
y=conv（x,h）;
结果：
实验二、离散傅立叶变换与快速傅立叶变换
1、离散傅立叶变换（DFT）
源程序：
functionpr4（）
F=50;
N=64;
T=0.000625;
n=1:
N;
x=cos（2*pi*F*n*T）;%x（n）=cos（pi*n/16）
subplot（2,1,1）;%绘制第一个子图x（n）
stem（n,x）;%绘制冲激
title（'x（n）'）;%标题为x（n）
xlabel（'n'）;%横座标为n
X=dft（x,N）;%调用dft函数计算x（n）的傅里叶变换
magX=abs（X）;%取变换的幅值
subplot（2,1,2）;%绘制第二个子图DFT|X|
stem（n,X）;
title（'DFT|X|'）;
xlabel（'f（pi）'）;%横座标为f（pi）
%dft的子程序
%实现离散傅里叶变换
function[Xk]=dft（xn,N）
n=0:
N-1;
k=0:
N-1;
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
结果：
F=50,N=64,T=0.000625时的波形
F=50,N=32,T=0.000625时的波形:
2、快速傅立叶变换（FFT）
源程序：
%functionpr5（）
F=50;N=64;T=0.000625;
n=1:
N;
x=cos（2*pi*F*n*T）;%x（n）=cos（pi*n/16）
subplot（2,1,1）;plot（n,x）;
title（'x（n）'）;xlabel（'n'）;%在第一个子窗中绘图x（n）
X=fft（x）;magX=abs（X）;
subplot（2,1,2）;plot（n,X）;
title（'DTFT|X|'）;xlabel（'f（pi）'）;%在第二个子图中绘图x（n）的快速傅%里叶变换
结果：
3、卷积的快速算法
源程序：
functionpr6（）
n=0:
15;
x=1.^n;
h=（4/5）.^n;
x（16:
32）=0;
h（16:
32）=0;
%到此x（n）=1,n=0~15;x（n）=0,n=16~32
%h（n）=（4/5）^n,n=0~15;h（n）=0,n=16~32
subplot（3,1,1）;
stem（x）;
title（'x（n）'）;
axis（[1,32,0,1.5]）;%在第一个子窗绘图x（n）横轴从1到32,纵轴从0到1.5（边框范围）
subplot（3,1,2）;
stem（h）;
title（'h（n）'）;
axis（[1,32,0,1.5]）;%在第二个子窗绘图h（n）横轴从1到32,纵轴从0到1.5
X=fft（x）;%X（n）为x（n）的快速傅里叶变换
H=fft（h）;%H（n）为h（n）的快速傅里叶变换
Y=X.*H;%Y（n）=X（n）*H（n）
%Y=conv（x,h）;
y=ifft（Y）;%y（n）为Y（n）的傅里叶反变换
subplot（3,1,3）%在第三个子窗绘图y（n）横轴从1到32,纵轴从0到6
stem（abs（y））;
title（'y（n=x（n）*h（n））'）;
axis（[1,32,0,6]）;
结果：
实验总结与思考
1、在较短的傅里叶变换中，FFT的计算速度与DFT相比不是很明显，序列计算长度越长，计算时间差距越大，FFT较快；
2、对于不同序列的较小长度的频谱分析可能会得到相同的频谱，适当加倍长度会避免这种情况的发生；
3、对同一序列的不同间隔的FFT变换，在满足奈奎斯特定律的情况下也会产生栅栏效应、频谱泄露、旁瓣效应等，采取适当的方法可以减弱这些不利效应；
4、在计算两个序列的离散卷积的时候要注意序列的长度L>=M+N-1。
实验三、IIR数字滤波器设计
源程序：
functionpr7（）
wp=0.2*pi;%频率转换（WP:
通带截止频率）
ws=0.3*pi;%频率转换（WS:
阻带截止频率）
Rp=1;%对于参数的说明
As=15;
T=1;
Fs=1/T;
OmegaP=（2/T）*tan（wp/2）;%OmegaP（w）=2*tan（0.1*pi）
OmegaS=（2/T）*tan（ws/2）;%OmegaS（w）=2*tan（0.15*pi）
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep.^2））;
Attn=1/10^（As/20）;
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;
OmegaC=OmegaP/（（10.^（Rp/10）-1）.^（1/（2*N）））;
[cs,ds]=u_buttap（N,OmegaC）;
[b,a]=bilinear（cs,ds,Fs）;
[mag,db,pha,w]=freqz_m（b,a）;
subplot（3,1,1）;%统共要绘制三个图，把此图放在第一个子窗绘制
%幅度响应的图形
plot（w/pi,mag）;
title（'幅度响应'）
xlabel（'w（pi）'）;
ylabel（'H'）;
axis（[0,1,0,1.1]）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]）;
set（gca,'YTickmode','manual','YTick',[0,Attn,Ripple,1]）;
grid;
subplot（3,1,2）;%在第二个子窗以分贝为单位绘制幅度响应的图形
plot（w/pi,db）;
title（'幅度响应（dB）'）;
xlabel（'w（pi）'）;
ylabel（'H'）;
axis（[0,1,-40,5]）;
set（gca,'XTickmode','manual','XTick',[0,0.2,0.35,1.1]）;
set（gca,'YTickmode','manual','YTick',[-50,-15,-1,0]）;
grid;
subplot（3,1,3）;%在第三个子
升级会员 