数字信号处理编程作业机电.docx
《数字信号处理编程作业机电.docx》由会员分享，可在线阅读，更多相关《数字信号处理编程作业机电.docx（33页珍藏版）》请在冰豆网上搜索。
数字信号处理编程作业机电
数字信号处理编程作业
姓名：
学号：
单位：
专业：
导师：
实验一离散时间系统及离散卷积
1、离散时间系统的单位脉冲响应
（1）选择一个离散时间系统；
y（n）=x（n）+2y（n-1）
y（n）=0,n<0;h（n）=0,n<0;
（2）用笔进行差分方程的递推计算；
当x（n）=
，y（n）=h（n），得
h（0）=1+2y（-1）=1;
h
（1）=0+2h（0）=2;
h
（2）=0+2h
（1）=4;
h（n）=0+2h（n-1）=2n.
从而，系统函数的单位冲击响应为
h（n）=2nu（n）;
（3）编制差分方程的递推计算程序；
n=[0:
9];
y=zeros（1,length（n））;
a=2;
form=1:
10
if（m==1）
y（1,1）=1;
else
y（1,m）=a*y（1,m-1）;
end
end
subplot（1,1,1）,stem（n,y）;title（'³å¼¤ÏìÓ¦'）;xlabel（'n'）;ylabel（'h（n）'）;
（4）在计算机上实现递推运算；
得到的结果，如图所示
（5）将程序计算结果与笔算的计算结果进行比较，验证程序运行的正确性；
>>y
y=
1248163264128256512
与笔算结果完全符合。
1-2参考实例验证：
离散时间系统的单位脉冲响应
%单位取样序列函数
function[x,n]=impseq（n0,n1,n2）
n=[n1:
n2];
x=[（n-n0）==0];
%单位阶跃序列函数
function[x,n]=stepseq（n0,n1,n2）
n=[n1:
n2];
x=[（n-n0）>=0];
（1）参考实例验证
>>[x,n]=impseq（0,-20,120）;
>>a=[1,-1,0.9];b=1;
>>h=filter（b,a,x）;
>>subplot（2,1,1）,stem（n,h）;title（'冲激响应'）;xlabel（'n'）;ylabel（'h（n）'）;
subplot（2,1,2）,zplane（b,a）;title（'系统零极点图'）;
2、离散系统的幅频、相频的分析方法
（1）参考实例验证
b=[0.0181,0.0543,0.0543,0.0181];
a=[1.000,-1.76,1.1829,-0.2781];
functionwork1（h,N）
sum=0;
w=0:
0.1:
pi;
forn=1:
N
sum=sum+h（n）.*exp（-j*n*w）;
end
Hm=abs（sum）;
Ha=angle（sum）;
subplot（2,1,1）;plot（w/pi,Hm）;title（'幅度响应'）;grid;
ylabel（'幅度'）;xlabel（'以\pi为单位的频率'）;
subplot（2,1,2）;plot（w/pi,Ha）;title（'相位响应'）;grid;
ylabel（'相位/pi'）;xlabel（'以\pi为单位的频率'）;
>>b=[0.0181,0.0543,0.0543,0.0181];a=[1.000,-1.76,1.1829,-0.2781];
[x,n]=impseq（0,0,100）;h=filter（b,a,x）;
>>work1（h）;
3、离散卷积的计算
（1）参考实例验证
%矩形序列函数
function[x,n]=rectseq（n0,n1,n2,n3）
n=[n2:
n3];
x=[（n-n0>=0）&（n-n1<0）];
>>[x,n0]=rectseq（0,10,-5,50）;[y,n1]=stepseq（0,-5,50）;h=y.*（0.9.^n1）;
c=conv（x,h）;
>>subplot（3,1,1）;stem（n0,x）;title（'输入序列'）;
xlabel（'）'）;ylabel（'x（n）'）;axis（[-5,50,0,2]）;
subplot（3,1,2）;stem（n1,h）;title（'冲激响应'）;
xlabel（'n'）;ylabel（'h（n）'）;axis（[-5,50,0,2]）;
subplot（3,1,3）;stem（-10:
100,c）;title（'输出'）;
xlabel（'n'）;ylabel（'y（n）'）;axis（[-10,100,0,8]）;
实验二离散傅立叶变换与快速傅立叶变换
一、实验内容
1、用离散傅立叶变换程序处理时间抽样信号，并根据实序列离散傅立叶变换的对称性，初步判定程序的正确性；
function[Xk]=dft（xn,N）
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk;
>>F=50;n=1:
64;T=0.000625;x=cos（2*pi*F*n*T）;X=dft（x,64）;
>>subplot（2,1,1）;stem（n,x）;xlabel（'n'）;ylabel（'x（n）'）;
subplot（2,1,2）;stem（n,X）;xlabel（'k'）;ylabel（'X（k）'）;
2、观察三角波和反三角波序列的时域和幅频特性，用N=8点FFT分析信号序列
和
的幅频特性，观察两者的序列形状和频谱曲线有什么异同？
绘出两序列及其幅频特性曲线。
三角波序列
反三角波序列
functionxc=xc（）
%三角波序列
forn=0:
7
i=n+1;
if（n>=0&&n<=3）
xc（i）=n+1;
elseif（n>=4）
xc（i）=8-n;
else
xc（i）=0;
end
end
end
N=8;
XC=fft（xc,N）;
n=0:
7;
k=0:
7;
subplot（2,2,1）;stem（n,xc）;title（'序列xc'）;
xlabel（'n'）;ylabel（'xc（n）'）;
subplot（2,2,2）;stem（k,XC）;title（'FFT变换'）;
xlabel（'k'）;ylabel（'XC（k）'）;
subplot（2,2,3）;stem（k,abs（XC））;title（'幅度'）;
xlabel（'k'）;ylabel（'|XC（k）|'）;
subplot（2,2,4）;stem（k,angle（XC））;title（'相位'）;
xlabel（'k'）;
>>xc;
functionxd=xd（）
%反三角波序列
forn=0:
7
i=n+1;
if（n>=0&&n<=3）
xd（i）=4-n;
elseif（n>=4）
xd（i）=n-3;
else
xd（i）=0;
end
end
end
N=8;
XD=fft（xd,N）;
n=0:
7;
k=0:
7;
subplot（2,2,1）;stem（n,xd）;title（'序列xd'）;
xlabel（'n'）;ylabel（'xd（n）'）;
subplot（2,2,2）;stem（k,XD）;title（'FFT变换'）;
xlabel（'k'）;ylabel（'XD（k）'）;
subplot（2,2,3）;stem（k,abs（XD））;title（'幅度'）;
xlabel（'k'）;ylabel（'|XD（k）|'）;
subplot（2,2,4）;stem（k,angle（XD））;title（'相位'）;
xlabel（'k'）;
>>xd;
从图中的幅频相位图可以看出，它们的幅度特性是偶对称的，相位特性是奇对称。
实验中的信号序列
和
，在单位圆上的Z变换频谱
和
会相同吗？
如果不同，你能说出哪一个低频分量更多一些吗？
为什么？
functionwork2（xc,xd,N）
XC=0;
XD=0;
w=0:
0.1:
pi;
forn=1:
N
XC=XC+xc（n）.*exp（-j*（n-1）*w）;
XD=XD+xd（n）.*exp（-j*（n-1）*w）;
end
Hmc=abs（XC）;Hmd=abs（XD）;
subplot（2,1,1）;plot（w/pi,Hmc）;title（'Xc（w）幅度响应'）;grid;
subplot（2,1,2）;plot（w/pi,Hmd）;title（'Xd（w）幅度响应'）;grid;
>>xc0=[1,2,3,4,4,3,2,1];xd0=[4,3,2,1,1,2,3,4];
>>work2（xc0,xd0,8）;
3、已知余弦信号如下
当信号频率
，采样间隔
，采样长度
时，对该信号进行傅立叶变换。
用FFT程序分析正弦信号，分别在以下情况下进行，并且分析比较结果
（1）F=50,N=32,T=0.000625；
（2）F=50,N=32,T=0.005；
（3）F=50,N=32,T=0.0046875；
（4）F=50,N=32,T=0.004；
（5）F=50,N=64,T=0.000625
function[x,X]=work3（F,N,T）
n=1:
N;
x=cos（2*pi*F*n*T）;
X=fft（x,N）;
functionwork3_1（x1,x2,x3,x4,x5,X1,X2,X3,X4,X5）
n1=1:
32;
n2=1:
64;
subplot（5,2,1）;stem（n1,x1）;xlabel（'n'）;ylabel（'x1（n）'）;
subplot（5,2,3）;stem（n1,x2）;xlabel（'n'）;ylabel（'x2（n）'）;
subplot（5,2,5）;stem（n1,x3）;xlabel（'n'）;ylabel（'x3（n）'）;
subplot（5,2,7）;stem（n1,x4）;xlabel（'n'）;ylabel（'x4（n）'）;
subplot（5,2,9）;stem（n2,x5）;xlabel（'n'）;ylabel（'x5（n）'）;
subplot（5,2,2）;stem（n1,X1）;xlabel（'k'）;ylabel（'X1（k）'）;
subplot（5,2,4）;stem（n1,X2）;xlabel（'k'）;ylabel（'X2（k）'）;
subplot（5,2,6）;stem（n1,X3）;xlabel（'k'）;ylabel（'X3（k）'）;
subplot（5,2,8）;stem（n1,X4）;xlabel（'k'）;ylabel（'X4（k）'）;
subplot（5,2,10）;stem（n2,X5）;xlabel（'k'）;ylabel（'X5（k）'）;
>>[x1,X1]=work3（50,32,0.000625）;
>>[x2,X2]=work3（50,32,0.005）;
>>[x3,X3]=work3（50,32,0.0046875）;
>>[x4,X4]=work3（50,32,0.004）;
>>[x5,X5]=work3（50,64,0.000625）;
>>work3_1（x1,x2,x3,x4,x5,X1,X2,X3,X4,X5）;
4、选定某一时间信号进行N=64点离散傅立叶变换，并且，对同一信号进行快速傅立叶变换，并比较它们的速度。
>>formatlong;tic;n=1:
64;x=cos（2*pi/32*n）;X1=dft（x,64）;toc;
Elapsedtimeis0.000000seconds.
>>formatlong;tic;n=1:
64;x=cos（2*pi/32*n）;X2=fft（x,64）;toc;
Elapsedtimeis0.000000seconds.
程序太小，程序运算时间约为0.
一、实验要求
1、调试实验程序，并且，给参考程序加注释；
2、完成实验内容2，并对结果进行分析。
实验中的信号序列
和
，在单位圆上的Z变换频谱
和
会相同吗？
如果不同，你能说出哪一个低频分量更多一些吗？
为什么？
3、完成实验内容3，并对结果进行分析；
4、利用编制的计算卷积的计算程序，分别给出一下三组函数的卷积结果。
（1）
（2）
（3）
、
functiony1=work4（）
n1=0:
14;x1=1.^n1;h1=（0.8）.^n1;%
（1）
n21=0:
9;n22=0:
19;x2=1.^n21;h2=0.5*sin（0.5*n22）;%
（2）
n3=0:
9;x3=1-（0.1）*n3;h3=（0.1）*n3;%（3）
%
（1）
l1=length（x1）+length（h1）-1;
X1=fft（x1,l1）;H1=fft（h1,l1）;
y1=ifft（X1.*H1）;
%
（2）
l2=length（x2）+length（h2）-1;
X2=fft（x2,l2）;H2=fft（h2,l2）;
y2=ifft（X2.*H2）;
%（3）
l3=length（x3）+length（h3）-1;
X3=fft（x3,l3）;H3=fft（h3,l3）;
y3=ifft（X3.*H3）;
%showresult
figure
（1）;
subplot（3,1,1）;stem（n1,x1）;xlabel（'n'）;ylabel（'x1（n）'）;
subplot（3,1,2）;stem（n1,h1）;xlabel（'n'）;ylabel（'h1（n）'）;
subplot（3,1,3）;stem（1:
l1,real（y1））;xlabel（'n'）;ylabel（'y1（n）'）;
figure
（2）;
subplot（3,1,1）;stem（n21,x2）;xlabel（'n'）;ylabel（'x2（n）'）;
subplot（3,1,2）;stem（n22,h2）;xlabel（'n'）;ylabel（'h2（n）'）;
subplot（3,1,3）;stem（1:
l2,real（y2））;xlabel（'n'）;ylabel（'y2（n）'）;
figure（3）;
subplot（3,1,1）;stem（n3,x3）;xlabel（'n'）;ylabel（'x3（n）'）;
subplot（3,1,2）;stem（n3,h3）;xlabel（'n'）;ylabel（'h3（n）'）;
subplot（3,1,3）;stem（1:
l3,real（y3））;xlabel（'n'）;ylabel（'y3（n）'）;
（1）
（2）
（3）
、
实验三IIR数字滤波器设计
一、实验内容如下：
1、设计一个巴特沃思数字低通滤波器，设计指标如下：
通带内
幅度衰减不大于1dB；
阻带
幅度衰减不小于15dB；
%非归一化巴特沃斯模拟滤波器原型函数
function[b,a]=buttap_o（N,omega）
[z,p,k]=buttap（N）;
p=p*omega;
k=k*omega^N;
B=real（poly（z））;
b0=k;b=k*B;a=real（poly（p））;
%非归一化巴特沃斯模拟滤波器设计函数
function[b,a]=afd_buttap（wp,ws,Rp,As）
ifwp<=0
error（'通带边缘必须大于0'）;
end
ifws<=0
error（'阻带边缘必须大于通带边缘'）;
end
if（Rp<=0）|（As<0）
error（'通带波动或阻带波动衰减必须大于0'）;
end
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（wp/ws）））;
fprintf（'\n***ButterworthFilterOrder=%2.0f\n',N）;
omega=wp/（（10^（Rp/10）-1）^（1/（2*N）））;
[b,a]=buttap_o（N,omega）;
%计算模拟滤波器频率响应的函数
function[db,mag,pha,w]=freqs_m（b,a,wmax）
w=[0:
1:
500]*wmax/500;
H=freqs（b,a,w）;
mag=abs（H）;
pha=angle（H）;
db=20*log10（（mag+eps）/max（mag））;
%冲激响应不变法设计IIR数字滤波器
function[b,a]=imp_invr（c,d,T）
[R,p,k]=residue（c,d）;
p=exp（p*T）;
[b,a]=residuez（R,p,k）;
b=real（b'）;a=real（a'）;
%计算离散系统频率响应
function[db,mag,pha,grd,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
1:
501））';w=（w（1:
1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
grd=grpdelay（b,a,w）;
2、编制计算设计的数字滤波器幅度特性和相位特性的程序，并进行实验验证。
>>wp=0.2*pi;ws=0.35*pi;Rp=1;As=15;T=1;
op=wp/T;os=ws/T;
[cs,ds]=afd_buttap（op,os,Rp,As）;
[b,a]=imp_invr（cs,ds,T）;
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度响应'）;grid;ylabel（'幅度'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,2）;plot（w/pi,pha）;title（'相位响应'）;grid;ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;grid;ylabel（'对数幅度/dB'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,4）;plot（w/pi,grd）;title（'群延迟'）;grid;ylabel（'样本'）;xlabel（'以\pi为单位的频率'）;
***ButterworthFilterOrder=5
二、实验要求
1、调试实验程序，并且，给实验程序加注释；
2、根据实验结果，给出自己设计的数字滤波器的幅度特性和相位特性；
3、用所设计的滤波器对不同频率的正弦波信号进行滤波，以说明其特性；
4、fp=0.2KHz,Rp=1dB,fs=0.3KHz,As=25dB,T=1ms；分别用脉冲响应不变法及双线性变换法设计一Butterworth数字低通滤波器，观察所设计数字滤波器的幅频特性曲线，记录带宽和衰减量，检查是否满足要求。
比较这两种方法的优缺点。
冲激响应不变法：
>>fp=200;fs=300;Rp=1;As=25;T=0.001;
wp=2*pi*fp*T;ws=2*pi*fs*T;
op=wp/T;os=ws/T;
[cs,ds]=afd_buttap（op,os,Rp,As）;
[b,a]=imp_invr（cs,ds,T）;
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度响应'）;grid;ylabel（'幅度'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,2）;plot（w/pi,pha）;title（'相位响应'）;grid;ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;grid;ylabel（'对数幅度/dB'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,4）;plot（w/pi,grd）;title（'群延迟'）;grid;ylabel（'样本'）;xlabel（'以\pi为单位的频率'）;
***ButterworthFilterOrder=9
双线性变换法：
>>fp=200;fs=300;Rp=1;As=25;T=0.001;
wp=2*pi*fp*T;ws=2*pi*fs*T;
op=（2/T）*tan（wp/2）;os=（2/T）*tan（ws/2）;
[cs,ds]=afd_buttap（op,os,Rp,As）;
[b,a]=bilinear（cs,ds,1/T）;
[db,mag,pha,grd,w]=freqz_m（b,a）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度响应'）;grid;ylabel（'幅度'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,2）;plot（w/pi,pha）;title（'相位响应'）;grid;ylabel（'相位'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,3）;plot（w/pi,db）;title（'幅度响应（dB）'）;grid;ylabel（'对数幅度/dB'）;xlabel（'以\pi为单位的频率'）;
subplot（2,2,4）;plot（w/pi,grd）;title（'群延迟'）;grid;ylabel（'样本'）;xlabel（'以\pi为单位的频率'）;
***ButterworthFilterOrder=6
实验四FIR数字滤波器设计
一、实验内容
1、设计一个FIR数字滤波器，设计指标如下：
通带内
幅度衰减不大于1dB；
阻带
幅度衰减不小于15dB；
%理想低通滤波器的计算函数
functionhd=ideal_lp（wc,M）;
alpha=（M-1）/2;
n=[0:
1:
（M-1）];
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
>>wp=0.2*pi;ws=0.35*pi;N=20;n=[0:
1:
N-1];
wc=（ws+wp）/2;%理想低通截止频率
hd=ideal_lp（wc,N）;%理想低通冲激响应
w_ham=（hamming（N））';%Hamming窗
h
升级会员 