数字信号作业报告.docx
《数字信号作业报告.docx》由会员分享，可在线阅读，更多相关《数字信号作业报告.docx（25页珍藏版）》请在冰豆网上搜索。
数字信号作业报告
数字信号作业报告
完成时间：
2012年5月1日星期二
3.53利用MATLAB语言，画出下列无限长脉冲响应系统的幅频特性曲线和相频响应特性曲线，并指出系统的类型。
（1）程序：
a=[010];
b=[1-21];
c=[1-0.40];
d=[1-0.880.61];
e=conv（a,b）;
f=conv（c,d）;
[h,w]=freqz（e,f,'whole'）;
subplot（2,1,1）,plot（w/pi,abs（h））;xlabel（'\omega/\pi'）;ylabel（'|h（e^j^\omega）|'）
subplot（2,1,2）,plot（w/pi,angle（h）/pi）;xlabel（'\omega/\pi'）;ylabel（'\phi（\omega/\pi'）
结果：
（2）程序：
a=[0.0534000];
b=[1100];
c=[1-1.0161];
d=[1-0.68300];
e=[1-1.44610.79570];
f=conv（c,c）;
f1=conv（a,b）;
f2=conv（f1,f）;
g=conv（e,d）;
[h,w]=freqz（f2,g,'whole'）;
subplot（2,1,1）,plot（w/pi,abs（h））;xlabel（'\omega/\pi'）;ylabel（'|h（e^j^\omega）|'）
subplot（2,1,2）,plot（w/pi,angle（h）/pi）;xlabel（'\omega/\pi'）;ylabel（'\phi（\omega/\pi'）
结果：
（3）程序：
a=[1-1];
b=[1-1.4990.8482];
c=[1-1.55480.6493];
f=conv（a,a）;
f1=conv（a,a）;
f2=conv（f1,f）;
%f2=conv（f,f）;
g=conv（b,c）;
[h,w]=freqz（f2,g,'whole'）;
subplot（2,1,1）,plot（w/pi,abs（h））;xlabel（'\omega/\pi'）;ylabel（'|h（e^j^\omega）|'）
subplot（2,1,2）,plot（w/pi,angle（h）/pi）;xlabel（'\omega/\pi'）;ylabel（'\phi（\omega/\pi'）
结果：
4.5在变换区间0<=n<=N-1，内计算以下各序列的N点离散傅里叶变换。
声明自己编的两个函数：
单位脉冲序列：
function[x,n]=impseq（np,ns,nf）
%生成x（n）=delta（n-np）;ns<=n<=nf
%-------------------------------
%ns为序列起始位置，nf为序列终止位置，np为脉冲位置
%调用方式[x,n]=impseq（np,ns,nf）
%
if（ns>np）||（ns>nf）||（np>nf）
error（'输入参数位置不满足ns<=np<=nf'）
elsen=[ns:
nf];
x=[（n-np）==0];
end
单位阶跃序列：
function[x,n]=stepseq（np,ns,nf）
n=[ns:
nf];
x=[（n-np）>=0];
（1）X（n）=1程序：
m=2;
n=2^m;
x=1;
y=fft（x,n）
结果：
y=
1111
y=
（2）X（n=δ（n）
程序;
M=4;
n=round（（2^M）*rand
（1））;%如果区间为（a,b）则可用
%round（（b-a）*rand
（1）+a）
m=round（n*rand
（1））;
y=impseq（m,0,2^M）;
fft（y,n）
结果：
ans=
111
（3）X（n）=δ（n-m）
M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
y=impseq（n-m,0,2^M）;
fft（y,n）
结果：
ans=
1.0-0.5000-0.8660i-0.5000+0.8660i
（4）X（n）=Rm（n）
程序：
M=4;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
y=stepseq（0,0,2^M）-stepseq（m,0,2^M）;
fft（y,n）
结果：
ans=
11
（5）X（n）=a^n
程序;M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
a=1/2;
b=0:
n;
y=a.^b;
fft（y,n）
a1=2;
y1=a1.^b;
fft（y1,n）
结果：
ans=
1.75000.6250-0.2165i0.6250+0.2165i
（6）X（n）=e^（j*2*pi/m*n）
程序：
M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
w=（2*pi/（2^M））*m;
y=cos（w*n）+1i*sin（w*n）;
y1=exp（（1i*2*pi/m*n））;
fft（y,n）
fft（y1,n）
结果：
ans=
1.0000-0.0000i1.0000-0.0000i1.0000-0.0000i1.0000-0.0000i
（7）X（n）=e^（j*w0*n）
程序:
M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
w0=10*pi;
y=exp（（1i*w0*n））;
fft（y,n）
结果：
ans=
1.0000-0.0000i1.0000-0.0000i
（8）X（n）=cos（2*pi/N*m*n）
程序：
M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
w=（2*pi/（2^M））*m;
y=cos（w*n）;
fft（y,n）
结果：
ans=
11111
（11）X（n）=nRN（n）;
程序：
M=2;
n=round（（2^M）*rand
（1））;
m=round（n*rand
（1））;
y=n*（stepseq（0,0,2^M）-stepseq（2^M,0,2^M））;
fft（y,5）
结果：
ans=
Columns1through4
8.0000-0.6180-1.9021i1.6180-1.1756i1.6180+1.1756i
Column5
-0.6180+1.9021i
（12）X（n）=ifn=奇数,1
Else0
程序;function[x,n]=jo（ns）
%生成x（n）=delta（n-np）;ns<=n<=nf
%-------------------------------
%ns为序列起始位置，nf为序列终止位置，np为脉冲位置
%调用方式[x,n]=impseq（np,ns,nf）
%
if（mod（ns,2）==1）
x=0;
elsex=ns;
end
M=4;
n=round（（2^M）*rand
（1））;
fft（jo（n））,n
结果：
ans=
22222
4.7已知序列
X（n）=[123321]
（1）求序列x（n）的离散时间傅里叶变换X（e^jw），画出幅频特性曲线和相频特性曲线。
（2）计算序列x（n）的N（N）6）点离散傅里叶变换X（k），画出幅频特性曲线和相频特性曲线。
（3）将X（e^jw）和X（k）的幅频特性曲线和相频特性曲线分别画在同一幅图中，验证X（k）是X（e^jw）[0,2*pi]区间上的N点等间隔采样，采样间隔为2*pi/N.
（4）计算X（k）的N点IDFT,验证DET和IDFT的唯一性。
程序：
xn=[123321];
x1=fft（xn,1024）;
y11=abs（x1）;
y12=angle（x1）;
figure
（1）
subplot（2,1,1）
stem（y11）
subplot（2,1,2）
stem（y12）
%--------------------------------------------------------------%
x2=fft（xn,32）;
y21=abs（x2）;
y22=angle（x2）;
figure
（2）
subplot（2,1,1）
stem（y21）
subplot（2,1,2）
stem（y22）
%--------------------------------------------------------------%
x3=fft（xn,32）;
y31=ifft（x3）;
figure（3）
subplot（2,1,1）
stem（x3）
subplot（2,1,2）
stem（y31）
%--------------------------------------------------------------%
x3=fft（xn,32）;
n=length（x3）;
m=n/（2*pi/32）;
x4=fft（xn,m）;
figure（4）
subplot（2,2,1）
stem（abs（x3））
subplot（2,2,2）
stem（abs（x4））
subplot（2,2,3）
stem（angle（x3））
subplot（2,2,4）
stem（angle（x4））
5.11利用128点的DFT和IDFT计算一个60点序列和一个1200点序列的线性卷积。
试确定利用重叠相加法计算上述线性卷积所需要的最少DFT和IDFT的次数。
测试主程序：
M=8;
%Èç¹ûÇø¼äΪ（a,b）Ôò¿ÉÓÃ
%round（（b-a）*rand
（1）+a
h=[];
x=[];
form=1:
60
n=round（（2^M）*rand
（1））;
h=cat（2,h,n）;
end
form1=1:
1200
n=round（（2^M）*rand
（1））;
x=cat（2,x,n）;
end
my1conv（h,x）
x1=h;
x2=x;
myconv（x1,x2）
子程序：
my1conv（h,x）（利用DFT实现的线性卷积）
functiony=my1conv（h,x）
figure
（1）
%h=input（'ÇëÊäÈë¸ø¶¨ÐòÁÐh='）;
%x=input（'ÇëÊäÈëÐòÁÐx='）;
lenx=length（x）;
lenh=length（h）;
N=lenx+lenh-1;
x=[xzeros（1,N-lenx）];
subplot（221）;
stem（x）;xlabel（'n'）;ylabel（'x（n）'）;gridon;
title（'ÐòÁÐx'）
h=[hzeros（1,N-lenh）];
subplot（222）;
stem（h）;xlabel（'n'）;ylabel（'h（n）'）;gridon;
title（'ÐòÁÐh'）
N=2^（ceil（log2（N）））;
H=fft（h,128）;
X=fft（x,128）;
Y=H.*X;
y=real（ifft（Y,128））;
subplot（223）;
stem（y）;xlabel（'n'）;ylabel（'y（n）'）;gridon;
title（'ÏßÐÔ¾í»ýºó'）
z=conv（x,h）;
subplot（224）;
stem（z）;xlabel（'n'）;ylabel（'z（n）'）;gridon;
title（'´«Í³¾í»ý·½·¨'）
end
myconv（传统的线性卷积）
functiony=myconv（x1,x2）
figure
（2）
%x1=input（'x1='）;
%x2=input（'x2='）;
N1=length（x1）;
M=length（x2）;
L=N1+M-1;
for（n=1:
L）
y（n）=0;
for（m=1:
M）
k=n-m+1;
if（k>+1&&k<=N1）
y（n）=y（n）+x2（m）*x1（k）;
end
end
end
y1=conv（x1,x2）;
nx1=0:
N1-1;
nx2=0:
M-1;
ny=0:
L-1;
subplot（221）;
stem（nx1,x1,'.k'）;xlabel（'n'）;ylabel（'x1（n）'）;gridon;
title（'ÐòÁÐx1'）
subplot（222）;
stem（nx2,x2,'.k'）;xlabel（'n'）;ylabel（'x2（n）'）;gridon;
title（'ÐòÁÐx2'）
subplot（223）;
stem（ny,y,'.k'）;xlabel（'n'）;ylabel（'y（n）'）;gridon;
title（'ÏßÐÔ¾í»ý'）
subplot（224）;
stem（y1）;xlabel（'n'）;ylabel（'y1'）;gridon;
title（'convÖ±½Ó¾í»ý'）
传统卷积方法的结果：
利用DFT实现的结果：
5.14设连续时间信号x（t）=cos（2*pi*f*t+θ）,其中，f=4khz,θ=pi/8.现用FFT对该信号进行频谱分析，要求信号的频率分辨率∆f《=10hz.
1.采样频率fs应取多高？
答：
fs>=2f;
2.采样时间tp应取多长？
答：
tp=1/fs;
3.FFT的变换区间取多少？
答：
N=2f/∆f;
4.仿真结果：
f=4000;
w0=pi/8;
xf=9;
fs=3*f;
t=0:
1/fs:
1;
x=cos（2*pi*f*t+w0）;
tp=1/fs;
n=2*f/xf;
%y=fft（x,n）;
[h,w]=freqz（x）;
subplot（2,1,1）,plot（w/pi,abs（h））;xlabel（'\omega/\pi'）;ylabel（'|x|'）
subplot（2,1,2）,plot（w/pi,angle（h）/pi）;xlabel（'\omega/\pi'）;ylabel（'\phi（x）'）
FIR滤波器的设计
（1）汉明窗:
Fs=2.5*10^6;
fp=10^6;fs=10^6+200000;rs=40;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Bt=ws-wp;
wc=（wp+ws）/2;
N=ceil（10*pi/Bt）
hmn=fir1（N-1,wc/pi,hamming（N））;
rs=100;a=1;
freqz（hmn,1,5000）;
axis（[0.71-10010]）
%N=63;
汉明窗:
（2）汉宁窗：
Fs=2.5*10^6;
fp=10^6;fs=10^6+200000;rs=40;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Bt=ws-wp;
wc=（wp+ws）/2;
N=ceil（10*pi/Bt）
hmn=fir1（N-1,wc/pi,hanning（N））;
rs=100;a=1;
freqz（hmn,1,5000）;
axis（[0.71-8010]）
%N=63;
（3）布莱克曼窗：
Fs=2.5*10^6;
fp=10^6;fs=10^6+150000;rs=40;
wp=2*pi*fp/Fs;
ws=2*pi*fs/Fs;
Bt=ws-wp;
wc=（wp+ws）/2;
N=ceil（10*pi/Bt）
hmn=fir1（N-1,wc/pi,blackman（N））;
rs=100;a=1;
freqz（hmn,2,5000）;
axis（[0.71-15010]）
%N=63;
滤波器滤波：
functionHd=ulbm
%ULBMReturnsadiscrete-timefilterobject.
%
%M-FilegeneratedbyMATLAB（R）7.10andtheSignalProcessingToolbox6.13.
%
%Generatedon:
01-May-201219:
55:
48
%
%FIRWindowLowpassfilterdesignedusingtheFIR1function.
%AllfrequencyvaluesareinMHz.
Fs=2.5;%SamplingFrequency
N=100;%Order
Fc=1;%CutoffFrequency
flag='scale';%SamplingFlag
%Createthewindowvectorforthedesignalgorithm.
win=hamming（N+1）;
%CalculatethecoefficientsusingtheFIR1function.
b=fir1（N,Fc/（Fs/2）,'low',win,flag）;
Hd=dfilt.dffir（b）;
%[EOF]
fs1=2.5*10^6;
%fs1=10^6+150000;
t=0:
1/fs1:
1;
%y=sin（500*10^3*2*pi*t）+sin（2*10^6*2*pi*t）;
%y=sin（500*10^3*2*pi*t）+5*sin（10^6*pi*2*t）;
y=sin（500*10^3*2*pi*t）+5*sin（1.5*10^6*2*pi*t）;
figure
（1）
y1=filter（b,1,y）;
subplot（2,1,1）
plot（y）
subplot（2,1,2）
plot（y1）
%___________________________________
[h1,w1]=freqz（y,'whole'）;
figure
（2）,plot（w1/pi,abs（h1））;xlabel（'\omega/\pi'）;ylabel（'|y|'）
[h2,w2]=freqz（y1,'whole'）;
figure（3）,plot（w2/pi,abs（h2））;xlabel（'\omega/\pi'）;ylabel（'|y1|'）
结果：
时域：
滤波前：
滤波后：
以上是通过FDATOOL工具箱实现。
实验总结：
技术方面：
对于前面的一些题目的练习，我想问题不是很大，就是对MATLAB软件的使用进一步了解。
好多题书上有现成的例题，首先读懂例题再按照题目要求完成。
在第四章的做DFT变换时，遇到了一些问题比如说是单位脉冲序列，单位阶跃信号的表示上，通过查阅资料完成对其的编译。
对于滤波器的设计刚开始通过资料上相应的程序，想自己设计出来，但是发现有些问题，然后使用了FDATOOL这个分析工具，我觉得我自己对滤波器的MATLAB实现还是有些问题，希望通过以后的学习可以得到提高。
个人收获方方面：
我觉得通过对数字信号的这些题目的编程，使我更加了解的一些数字信号过程用软件是如何实现的，起初我认为实现起来是比较简单的，但是在自己的亲身体验下，我才发现了实现起来确实不易，通过网络，一些相关书籍的查询对某些具体问题有了较为深刻的理解。
比如说对低通滤波器的设计，虽然通过自己的努力没有实现但是我想重要的是过程而不是结果。
当然在此次作业中我对MATLAB的使用有了新的认识，这个软件在工程应用上很强大。
所以我想在以后的学习中我要认真的学习MATLAB相关应用，以便在以后的工程实践中，充分的利用解决实际问题。
当然我也很多不足比如说我对滤波器设计的一些具体问题没有很好的理解到位，我想这与我理论知识的掌握不扎实是有关系的。
另一方面我想对一前没有重视对MATLAB的使用是有关系的。
希望我能通过自己的努力可以提高相关能力。
升级会员 