数字信号实验.docx
《数字信号实验.docx》由会员分享，可在线阅读，更多相关《数字信号实验.docx（15页珍藏版）》请在冰豆网上搜索。
数字信号实验
实验步骤1、编写nonrec.m函数文件，实现FIR滤波。
这里给定，，求。
2、编写rec.m函数文件，实现纯递归IIR滤波。
这里给定，，，，求单位取样响应。
3、用nonrec.m和rec.m函数编制dfilter.m函数文件，构成完整的一般IIR滤波器程序，并完成下列信号滤波这里给定系统函数，计算。
4、用help查看内部函数filter.m，了解调用格式，重做3，并和你编写的dfilter.m的结果进行比较。
实验内容1、编程实现FIR滤波2、编写时实现纯递归IIR滤波；差分方程：
系统函数：
3、调用库函数filter.m实现IIR滤波
实验数据1.nonrec.m函数文件functiony=nonrec（h,x）m=length（h）;n=length（x）;fork=1:
（m+n-1）S=0;ifk<=nn1=k;form1=1:
mifn1>=1S=S+h（m1）*x（n1）;y（k）=S;n1=n1-1;endendelsen2=n;form1=（1+k-n）:
mifn2>=1S=S+h（m1）*x（n2）;y（k）=S;n2=n2-1;endendendendend>>h=ones（1,8）;>>x=0:
15;>>y=nonrec（h,x）y=0136101521283644526068768492847565544229152.rec.m函数文件functiony=rec（x,a）La=length（a）;Lx=length（x）;y
（1）=x
（1）;fork=2:
100A=0;ifk<=Lxifk<=La+1form=1:
k-1A=A+a（k-m）*y（m）;y（k）=A+x（k）;endelseforc=1:
LaA=A+a（c）*y（k-c）;y（k）=x（k）+A;endendelseford=1:
LaA=A+a（d）*y（k-d）;y（k）=A;endendendendend
3.dfilter.m函数文件functiony=dfilter（x,a,b）y1=rec（x,a）;y=nonrec（b,y1）;
>>a1=2*0.95*cos（pi/8）;>>a2=-0.95*0.95;>>a=[a1a2];x=[1];>>y=rec（x,a）;>>a=[0.6,-0.5]>>b=[1-31];>>n=0:
63;>>x=sin（3*pi*n/5）;>>y1=dfilter（x,a,b）;>>y2=filter（b,a,x）;>>subplot（2,1,1）,stem（y1）,gridon,title（'dfilter.m计算结果'）>>subplot（2,1,2）,stem（y2）,gridon,title（'filter.m计算结果'）
内部函数filter.m的调用格式是y=filter（b,a,X）编写的dfilter.m调用格式是y=dfilter（x,a,b）实验总结通过本此上机实践熟悉和掌握了数字信号出来的因果性数字系统的时域实现，进一步了解和认识了数字信号的MATLAB实现及用MATLAB实现数字信号处理。
176
本文来自网学（），转载请注明出处：
实验一离散傅里叶变换的性质及应用
一、实验目的
1、了解DFT的性质及应用。
2、熟悉MATLAB编程的特点。
二、实验内容
1、用三种不同的DFT程序计算x（n）=R8（n）的傅里叶变换X（ejw）,并比较三种程序计算机运行时间。
（1）用forloop语句的M函数文件dft1.m，用循环变量逐点计算X（k）；
（2）编写用MATLAB矩阵运算的M函数文件dft2.m,完成下列矩阵运算；
（3）调用FFT库函数，直接计算X（k）；
（4）分别利用上述三种不同方式编写的DFT程序计算序列x（n）的傅立叶变换X（ejw），并画出相应的幅频和相频特性，再比较各个程序的计算机运行时间。
M函数文件如下：
dft1.m:
function[Am,pha]=dft1（x）
N=length（x）;
w=exp（-j*2*pi/N）;
fork=1:
N
sum=0;
forn=1:
N
sum=sum+x（n）*w^（（k-1）*（n-1））;
end
Am（k）=abs（sum）;
pha（k）=angle（sum）;
end
dft2.m:
function[Am,pha]=dft2（x）
N=length（x）;
n=[0:
N-1];
k=[0:
N-1];
w=exp（-j*2*pi/N）;
nk=n'*k;
wnk=w.^（nk）;
Xk=x*wnk;
Am=abs（Xk）;
pha=angle（Xk）;
dft3.m:
function[Am,pha]=dft3（x）
Xk=fft（x）;
Am=abs（Xk）;
pha=angle（Xk）;
源程序及运行结果：
（1）x=[ones（1,8）,zeros（1,248）];
t=cputime;
[Am1,pha1]=dft1（x）;
t1=cputime-t
n=[0:
（length（x）-1）];
w=（2*pi/length（x））*n;
figure
（1）
subplot（2,1,1）,plot（w,Am1,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
subplot（2,1,2）,plot（w,pha1,'r'）;grid;
title（'PhasePart'）;
xlabel（'frequencyinradians'）;
ylabel（'argX[exp（jw）]/radians'）;
运行时间：
t1=2.25s
（2）x=[ones（1,8）,zeros（1,248）];
t=cputime;
[Am2,pha2]=dft2（x）;
t2=cputime-t
n=[0:
（length（x）-1）];
w=（2*pi/length（x））*n;
figure
（2）
subplot（2,1,1）,plot（w,Am2,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
subplot（2,1,2）,plot（w,pha2,'r'）;grid;
title（'PhasePart'）;
xlabel（'frequencyinradians'）;
ylabel（'argX[exp（jw）]/radians'）;
运行时间：
t2=0.77s
（3）x=[ones（1,8）,zeros（1,248）];
t=cputime;
[Am3,pha3]=dft3（x）;
t3=cputime-t;
n=[0:
（length（x）-1）];
w=（2*pi/length（x））*n;
figure（3）
subplot（2,1,1）,plot（w,Am3,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
subplot（2,1,2）,plot（w,pha3,'r'）;grid;
title（'PhasePart'）;
xlabel（'frequencyinradians'）;
ylabel（'argX[exp（jw）]/radians'）;
运行时间：
t3=0
从以上运行结果可以看出，调用FFT库函数直接计算X（k）速度最快，矩阵运算次之，用循环变量逐点计算运行速度最慢。
因此FFT算法大大提高了DFT的实用性。
2、利用DFT实现两序列的卷积运算，并研究DFT点数与混叠的关系。
给定x（n）=nR16（n），h（n）=R8（n），用FFT和IFFT分别求线性卷积和混叠结果输出，并用函数stem（n,y）画出相应图形。
源程序如下：
N=16;
x=[0:
N-1];
h=ones（1,8）;
线性卷积：
Xk1=fft（x,23）;%做23点fft
Hk1=fft（h,23）;
Yk1=Xk1.*Hk1;
y1=ifft（Yk1）;
n=0:
22;
figure
（1）
stem（n,y1）
混叠：
Xk2=fft（x）;
Hk2=fft（h,16）;%做16点fft
Yk2=Xk2.*Hk2;
y2=ifft（Yk2）;
n=0:
15;
figure
（2）
stem（n,y2）
运行结果：
线性卷积：
混叠结果：
将两个信号的DFT相乘，再做反变换，相当于时域做循环卷积。
而循环卷积是线性卷积以N为周期周期性延拓的结果，其中N为离散傅立叶变换的点数。
若x1（n）是N1点序列，x2（n）是N2点序列，则其线性卷积为N1+N2-1点，做N=max（N1,N2）点的傅立叶变换，则结果中将有N2-1点是混叠的。
题目中N1=16,N2=8,由程序运行结果可以看出，第二张图的前7（N2-1）个点为混叠结果，后9个点是线性卷积的结果。
所以，要用DFT来做线性卷积，DFT的点数必须大于或等于线性卷积的结果。
本题中取N1+N2-1=23为DFT的点数，即可得到正确的线性卷积结果。
3、研究高密度频谱与高分辨率频谱
设有连续信号
Xa（t）=cos（2π×6.5×10^3t）+cos（2π×7×10^3t）+cos（2π×9×10^3t）
以采样频率fs=32kHz对信号Xa（t）采样,分析下列三种情况的幅频特性。
（1）采集数据长度N=16点，做N=16点的DFT，并画出幅频特性。
（2）采集数据长度N=16点，补零到256点，做N=256点的DFT，并画出幅频特性。
（3）采集数据长度N=256点，做N=256点的DFT，并画出幅频特性。
观察三幅不同频率特性图，分析和比较它们的特点以及形成的原因。
源程序如下：
fs=32000;%采样频率
N=16;%采集16点
n=0:
N-1;%做16点DFT
xa=cos（2*pi*6.5*10^3*n/fs）+cos（2*pi*7*10^3*n/fs）+cos（2*pi*9*10^3*n/fs）;
[Am1,pha1]=dft3（xa）;
n=[0:
（length（xa）-1）];
w=（2*pi/length（xa））*n;
figure
（1）
plot（w,Am1,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
x=[xa（1:
16）,zeros（1,240）];%补零到256
[Am2,pha2]=dft2（x）;%做256点DFT
n=[0:
（length（x）-1）];
w=（2*pi/length（x））*n;
figure
（2）
plot（w,Am2,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
N0=256;%采集256点
n0=0:
N0-1;%做256点DFT
xa0=cos（2*pi*6.5*10^3*n0/fs）+cos（2*pi*7*10^3*n0/fs）+cos（2*pi*9*10^3*n0/fs）;
[Am3,pha3]=dft3（xa0）;
n=[0:
（length（xa0）-1）];
w=（2*pi/length（xa0））*n;
figure（3）
plot（w,Am3,'b'）;grid;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
运行结果：
N=16点DFT幅频特性
高密度频谱：
高分辨率频谱：
观察运行结果可知，N=16点时的DFT，由于点数太少，很难反映出频谱的细节特征，只能分辨出两个频率分量。
将16点DFT补零到256点，做256点DFT,由结果可以看出，频率间隔缩小（由2?
/16减小为2256），得到了比较平滑的连续曲线，此为高密度频谱。
它是由16点经过补零所得到的，并没有增加任何信息，频率分辨率并未提高，仍然只能看出两个频率分量。
第三次做的是256点的DFT，采集点数大大增加，分辨率提高，可以清楚看到三根谱线的位置，是为高分辨率频谱。
4、实现序列的内插和抽取所对应的傅里叶变换。
给定序列x（n）=[cos（π/36×n）+cos（1.5π/36×n）]R128（n），做Nfft=128点的傅里叶变换，并求
x1（n）=x（4n）和x2（n）=x（n/4）n=4kx2（n）=0n≠4k
对应的傅里叶变换（N=128点）。
比较这三个计算结果得到的幅频特性图，分析其差别产生的原因。
源程序如下：
N=128;%原序列128点变换
t=0:
N-1;
x=cos（pi/36*t）+cos（1.5*pi/36*t）;
Xk=fft（x）;
Am=abs（Xk）;
n=[0:
（length（x）-1）];
w=（2*pi/length（x））*n;
figure
（1）
plot（w,Am,'b'）;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
N=32;%抽取32点
t=0:
N-1;
x2=cos（pi/9*t）+cos（1.5*pi/9*t）;
Xk2=fft（x2,128）;%做128点傅里叶变换
Am2=abs（Xk2）;
n=[0:
127];
w=（2*pi/128）*n;
figure
（2）
plot（w,Am2,'b'）;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
x3=zeros（1,128）;%初始化数组为零
N=32;
fori=0:
N-1%n=4i时赋值
x3（1,4*i+1）=cos（pi/36*i）+cos（1.5*pi/36*i）;
end
Xk3=fft（x3）;%做128点fft
Am3=abs（Xk3）;
n=[0:
127];
w=（2*pi/128）*n;
figure（3）
plot（w,Am3,'b'）;
title（'Magnitudepart'）;
xlabel（'frequencyinradians'）;
ylabel（'|X（exp（jw））|'）;
原序列离散傅里叶变换：
抽取后的离散傅里叶变换：
内插后的离散傅里叶变换：
第一张图为原序列的傅立叶变换。
将原序列经过4倍降采样率后得到第二幅图，该频谱本应为原序列频谱移位后混叠相加的结果，但由于x1（n）取了128点，在32点之后全部为零，信息量减少，而形成了高密度频谱，大致形状与原频谱有所出入。
第三幅图为经过4倍升采样率的结果，相当于频率轴经过4倍的压缩，一个周期内有四个频带经过压缩的原频谱复现。
经过低通滤波可得到原频谱频带压缩后的结果。
三、思考题
1、MATLAB在处理数据运算时与其他高级语言相比有何不同？
答：
MATLAB是以复数矩阵为基本编程单元的一种程序设计语言，我们可以直接以矩阵和数组为单位，进行数组和矩阵的运算。
而其他高级语言没有这个功能。
从第一题的计算中就可以看出，使用矩阵运算比逐点计算的速度要快得多。
另外，MATLAB中有多种数学、工程类的函数可以调用，如这次使用的fft函数，极大的方便了我们进行工程计算。
2、什么是高密度频谱？
什么是高分辨率频谱？
实验内容3中为区分三个频率分量至少应采集几个样本？
答：
高密度频谱是指频谱的频率间隔较小，连续性好的谱。
但频率间隔小通过补零就可以得到，因此它不能反映频率分辨率大小。
高分辨率频谱是指可以很明显的分开各个信号频谱，它与采样点数有直接关系，采样点数越多，分辨率越好。
DFT的频率分辨率为fs/N,实验3中频率最接近的两个频率分量间隔为f1-f2=500Hz,fs=32k，可得出N=32k/500=64，这是最少的采样点数。
3、实验内容4中，分别求X1（ejw）和x2（ejw）与原X（ejw）的关系，并说明为什么抽取后能分辨出两个频率分量？
有无任何条件的限制？
答：
X1（Z）=1/4ΣX（e-j2πk/4×Z1/4）X1（ejw）＝1/4ΣX（ej（w-2πk）/4）
X2（ejw）=X（ej4w）
x（n）经抽取之后，相当于频谱沿频率轴拉伸了4倍，因此原来靠得很近而不能分辨的两根谱线就有可能分开而看清楚。
实验二因果性数字系统的时域实现
一.实验目的
编制实现下列差分方程的程序:
y（n）=∑bkx（n-k）+∑aky（n-k）
二.实验内容
1.编制nonrec.m函数文件,实现FIR滤波y（n）=h（n）*x（n）.这里给定h（n）=R8（n）,x（n）=nR16（n）,求y（n）.
nonrec.m函数文件:
functiony=nonrec（x,h）
x=[x,zeros（1,length（h）-1）];%补零
w=zeros（1,length（h））;
forI=1:
length（x）
forj=length（h）:
-1:
2
w（j）=w（j-1）;%得到每一延时单元输出变量
end
w
（1）=x（i）;
y（i）=w*h’;%Bcofi与wi对应相乘
end
主程序文件:
x=0:
15;
h=ones（1,8）;
y=nonrec（x,h）;
n=0:
22;
stem（n,y）;
分析：
线性卷积y（n）=x（n）*h（n）的长度为16+8-1=23，可利用y（n）=∑h（m）x（n-m）直接计算得
n（n+1）/2,n≤7
y（n）=4（2n-7）,8≤n≤15
（n+8）（23-n）/2,16≤n≤22
即y=[01361015212836445260687684928475
6554422915],与曲线相符。
2.编制rec.m函数文件,实现纯递归IIR滤波y（n）=x（n）+∑aky（n-k）.这里给定a1=2rcosw0,a2=-r2,r=0.95,w0=π/8,求单位抽样响应h（n）.
rec.m函数文件:
functiony=rec（x,a,n）
x=[x,zeros（1,n-length（x））];%补零到所需长度
sum=0;
w=zeros（1,length（a））;
fori=1:
n
y（i）=sum+x（i）;%递归
forj=length（a）:
-1:
2
w（j）=w（j-1）;%延时
end
w
（1）=y（i）;
sum=w*a';
end
主程序文件:
x=[1];
a=[2*0.95*cos（pi/8）,-0.95^2];
h=rec（x,a,75）;%取h（n）的长度为75点
n=0:
74;
stem（n,h）;
分析计算:
由题意,a1=2*0.95*cos（π/8）,a2=-0.952,所以,
得到系统函数H（z）=1/[1-1.9cos（π/8）z-1+0.952z-2],
做逆Z变换得h（n）=0.95ncos（πn/8）+ctg（π/8）*0.95nsin（πn/8）,
利用MATLAB直接画h（n）,即使用下列语句
n=0:
74;
h=0.95.^n.*cos（pi.*n./8）+cot（pi/8）.*（0.95.^n）.*sin（pi.*n./8）;
stem（n,h）;
得到如下曲线
此曲线与用rec函数所画曲线完全一致,可见,用MATLAB编制的FIR滤波程序与理论计算所得结果是相同的.
3.用nonrec.m和rec.m函数编制dfilter.m函数文件,构成完整的一般IIR滤波程序,并完
成下列信号滤波:
x（n）=cos（2πn/5）R64（n）
这里给定系统函数
H（z）=（1-2z-1+z-2）/（1-0.5z-1+0.5z-2）
求y（n）.
dfilter.m函数文件:
functiony=dfilter（x,b,a,n）
y1=nonrec（x,b）;
y=rec（y1,a,n）;
主程序文件:
n=0:
63;
x=cos（2*pi/5*n）;
b=[1,-2,1];%由H（z）得到系数a,b
a=[0.5,-0.5];
y=dfilter（x,b,a,64）;%取y（n）的长度为64点
stem（n,y）;
4.用help查看内部函数filter.m,了解调用格式,重做3,并与你编写的dfilter.m结果进行比较.
用help可以看到内部函数为Y=FILTER（B,A,X）,且有Thefilterisa"DirectFormIITransposed"implementationofthestandarddifferenceequation:
a
（1）*y（n）=b
（1）*x（n）+b
（2）*x（n-1）+...+b（nb+1）*x（n-nb）
-a
（2）*y（n-1）-...-a（na+1）*y（n-na）
因此,调用内部函数filter时,要对原系数a做适当变化:
a1
（1）=1,a1（i）=-a（i）,i>1.
n=0:
63;
x=cos（2*pi/5*n）;
b=[1,-2,1];
a=[0.5,-0.5];
y=dfilter（x,b,a,64）;%调用自己编的dfilter函数
a1=[1,-0.5,0.5];%a变为a1,用于调用内部函数filter
y1=filter（b,a1,x）;
subplot（2,1,1）;
stem（n,y）;
subplot（2,1,2）;
stem（n,y1）;
（曲线见下页）
由上两图可以看到,两函数所得结果完全一致,由此验证了所编dfilter函数的正确
5.重叠保留法做线性卷积
x=0:
134;%x为长序列
y=[];
h=[ones（1,8）,zeros（1,8）];%h为单位抽样响应
H
升级会员 