数字信号处理西电上机实验.docx
《数字信号处理西电上机实验.docx》由会员分享，可在线阅读，更多相关《数字信号处理西电上机实验.docx（31页珍藏版）》请在冰豆网上搜索。
数字信号处理西电上机实验
数字信号处理实验报告
实验一：
信号、系统及系统响应
一、实验目的：
（1）熟悉连续信号经理想采样前后的频谱变化关系，加深对时域采样定理的理解。
（2）熟悉时域离散系统的时域特性。
（3）利用卷积方法观察分析系统的时域特性。
（4）掌握序列傅里叶变换的计算机实现方法，利用序列的傅里叶变换对连续信号、离散信号及系统响应进行频域分析。
二、实验原理与方法：
（1）时域采样。
（2）LTI系统的输入输出关系。
三、实验内容、步骤
（1）认真复习采样理论、离散信号与系统、线性卷积、序列的傅里叶变换及性质等有关内容，阅读本实验原理与方法。
（2）编制实验用主程序及相应子程序。
①信号产生子程序，用于产生实验中要用到的下列信号序列：
a.xa（t）=A*e^-at*sin（Ω0t）u（t）
A=444.128;a=50*sqrt
（2）*pi;
b.单位脉冲序列：
xb（n）=δ（n）
c.矩形序列：
xc（n）=RN（n）,N=10
②系统单位脉冲响应序列产生子程序。
本实验要用到两种FIR系统。
a.ha（n）=R10（n）;
b.hb（n）=δ（n）+2.5δ（n-1）+2.5δ（n-2）+δ（n-3）
③有限长序列线性卷积子程序
用于完成两个给定长度的序列的卷积。
可以直接调用MATLAB语言中的卷积函数conv。
conv用于两个有限长度序列的卷积，它假定两个序列都从n=0开始。
调用格式如下：
y=conv（x,h）
四、实验内容
调通并运行实验程序，完成下述实验内容：
①分析采样序列的特性。
a.取采样频率fs=1kHz,即T=1ms。
b.改变采样频率，fs=300Hz，观察|X（ejω）|的变化，并做记录（打印曲线）；进一步降低采样频率，fs=200Hz，观察频谱混叠是否明显存在，说明原因，并记录（打印）这时的|X（ejω）|曲线。
②时域离散信号、系统和系统响应分析。
a.观察信号xb（n）和系统hb（n）的时域和频域特性；利用线性卷积求信号xb（n）通过系统hb（n）的响应y（n），比较所求响应y（n）和hb（n）的时域及频域特性，注意它们之间有无差别，绘图说明，并用所学理论解释所得结果。
b.观察系统ha（n）对信号xc（n）的响应特性。
③卷积定理的验证。
五、思考题
（1）在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同?
它们所对应的模拟频率是否相同?
为什么?
（2）在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？
为什么？
六、参考资料
数字信号处理教程——Matlab释义与实现第二版，陈怀琛著，电子工业出版社，2008年
七、实验程序及结果
本实验使用自定义函数的方法产生信号：
门函数：
function[y,n]=gate（np,ns,nf）
ifns>np|ns>nf|np>nf
error（'ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf'）;
else
n=[ns:
nf];
y=（n>=0&nend
冲击函数：
function[y,n]=impseq（np,ns,nf）
ifns>np|ns>nf|np>nf
error（'ÊäÈëλÖòÎÊý²»Âú×ãns<=np<=nf'）;
else
n=[ns:
nf];
y=[（n-np）==0];
end
时域函数：
function[y,ts]=sig（t0,tp,t1）;
A=444.128;a=50.*sqrt
（2）.*pi;w0=50.*sqrt
（2）.*pi;
ts=t0:
tp:
t1;
y=A.*exp（-a.*ts）.*sin（w0.*ts）.*（ts>=0）;
;hb=impseq（0,0,10）+2.5.*impseq（1,0,10）+2.5.*impseq（2,0,10）+impseq（3,0,10）;
;[ha,N]=gate（10,0,10）;
第一部分主程序：
[y1,t1]=sig（0,0.001,1）;
%t1=0:
0.001:
1;
%y1=0.5.*sin（2.*pi.*15.*t1）;
figure;
plot（t1,y1）;
axis（[00.15-10150]）;
title（'f=1000ʱÓòÐźÅ'）;
figure;
f=-500:
500;
yy1=fftshift（fft（y1））;
plot（f,abs（yy1））;
title（'f=1000ƵÆ×'）;
[y2,t2]=sig（0,1/300,1）;
figure;
plot（t2,y2）;
axis（[00.15-10150]）;
title（'f=300ʱÓòÐźÅ'）;
figure;
yy2=fftshift（fft（y2））;
plot（-150:
150,abs（yy2））;
title（'f=300ƵÆ×'）;
[y3,t3]=sig（0,0.02,1）;
figure;
plot（t3,y3）;
title（'f=200ʱÓòÐźÅ'）;
axis（[00.15-85]）;
figure;
yy3=fftshift（fft（y3））;
plot（-25:
25,abs（yy3））;
title（'f=200ƵÆ×'）;
第二部分主程序：
xb=impseq（0,0,100）;
figure;
stem（0:
length（xb）-1,xb）;
axis（[-1,3,0,2]）;
figure;
y1=fft（xb）;
plot（0:
length（y1）-1,y1）;
hb=impseq（0,0,100）+2.5.*impseq（1,0,100）+2.5.*impseq（2,0,100）+impseq（3,0,100）;
figure;
stem（0:
length（hb）-1,hb）;
axis（[-1,5,0,3]）
figure
y2=fftshift（fft（hb））;
plot（-（length（y2）-1）/2:
（length（y2）-1）/2,abs（y2））;
figure;
y3=conv（double（xb）,double（hb））;
plot（0:
length（y3）-1,y3）;
y4=fftshift（fft（y3））;
figure;
plot（-（length（y4）-1）/2:
（length（y4）-1）/2,abs（y4））;
xc=gate（10,0,100）;
ha=xc;
y5=fftshift（fft（conv（double（xc）,double（ha））））;
figure;
plot（-（length（y5）-1）/2:
（length（y5）-1）/2,abs（y5））;
实验结果：
第一部分：
采样频率为200Hz时时域恢复200Hz时频谱
采样频率为300Hz时时域恢复300Hz时频谱
采样频率为1000Hz是时域恢复1000Hz时频谱
第二部分：
Xb频域分析xb时域分析
hb频域分析hb时域分析
信号通过xb频域信号通过xc频域
卷积定理的验证：
%%%%%卷积定理的验证
y1=[0,1,2,3,4,5,6,7,8];
y2=[8,7,6,5,4,3,2,1,0];
y3=conv（y1,y2）;
f1=fft（y1,20）;
f2=fft（y2,20）;
f3=fft（y3,20）;
f=f1.*f2;
subplot（2,1,1）;
stem（0:
length（f）-1,f,'.'）;
subplot（2,1,2）;
stem（0:
length（f3）-1,f3,'.'）;
时域卷积后求频谱和频域相乘
可见，f=conv（y1,y2）的频谱和y1,y2的频谱相乘后结果相同。
即满足卷积定理
八、回答问题
（1）在分析理想采样序列特性的实验中，采样频率不同时，相应理想采样序列的傅里叶变换频谱的数字频率度量是否都相同？
它们所对应的模拟频率是否相同？
为什么？
答：
数字频率度量不相同，但他们所对应的模拟频率相同。
由w=Ω*Ts公式得，采样间隔变化时模拟频率对应的数字频率会有相应的变化，故其度量会有所变化。
而且采样频率的大小直接关系到能否将能否将原始信号恢复出来。
（2）在卷积定理验证的实验中，如果选用不同的频域采样点数M值，例如，选M=10和M=20，分别做序列的傅里叶变换，求得的结果有无差异？
答：
有差异，所到的结果点数不同。
九、总结
本实验主要是后续实验的基础。
涉及内容也比较浅显。
单位冲击序列与hb卷积后得到的频谱与hb原频谱相同。
原因很简单，是因为单位冲击序列卷积任何函数仍然是原函数。
实验二：
用FFT作谱分析
一、实验目的
（1）进一步加深DFT算法原理和基本性质的理解（因为FFT只是DFT的一种快速算法，所以FFT的运算结果必然满足DFT的基本性质）。
（2）熟悉FFT算法原理和FFT子程序的应用。
（3）学习用FFT对连续信号和时域离散信号进行谱分析的方法，了解可能出现的分析误差及其原因，以便在实际中正确应用FFT。
二、实验步骤
（1）复习DFT的定义、性质和用DFT作谱分析的有关内容。
（2）复习FFT算法原理与编程思想，并对照DIT-FFT运算流图和程序框图，读懂本实验提供的FFT子程序。
（3）编制信号产生子程序，产生以下典型信号供谱分析用：
（4）编写主程序。
下图给出了主程序框图，供参考。
本实验提供FFT子程序和通用绘图子程序。
（5）按实验内容要求，上机实验，并写出实验报告。
三、实验内容
（1）对2中所给出的信号逐个进行谱分析。
（2）令x（n）=x4（n）+x5（n），用FFT计算8点和16点离散傅里叶变换，
X（k）=DFT［x（n）］
（3）令x（n）=x4（n）+jx5（n），重复
（2）。
四、思考题
（1）在N=8时，x2（n）和x3（n）的幅频特性会相同吗?
为什么?
N=16呢?
（2）如果周期信号的周期预先不知道，如何用FFT进行谱分析?
五、参考资料
数字信号处理教程——Matlab释义与实现第二版，陈怀琛著，电子工业出版社，2008年
数字信号处理第三版，高西泉等著，西安电子科技大学出版社，第三版
六、实验程序及结果
x1=[1,1,1,1,0,0,0,0,0,0];
x2=[1,2,3,4,4,3,2,1,0,0];
x3=[4,3,2,1,1,2,3,4,0,0];
n=0:
1000;
x4=cos（n）;
x5=sin（n）;
t=0:
0.001:
1;
x6=cos（8*pi*t）+cos（16*pi*t）+cos（20*pi*t）;
%figure;
%title（'Öð¸öÆ×·ÖÎö'）
subplot（3,2,1）;
stem（0:
length（x1）-1,x1）;
subplot（3,2,2）;
stem（0:
length（x2）-1,x2）;
subplot（3,2,3）;
stem（0:
length（x3）-1,x3）;
subplot（3,2,4）;
stem（0:
20,x4（1:
21））;
subplot（3,2,5）;
stem（0:
20,x5（1:
21））;
subplot（3,2,6）;
stem（0:
20,x6（1:
21））;
figure;
stem（0:
200,x6（1:
201））;
y1=abs（fft（x1,8））;
y2=abs（fft（x2,8））;
y3=abs（fft（x3,8））;
y4=abs（fft（x4,8））;
y5=abs（fft（x5,8））;
y6=abs（fft（x6,8））;
figure;
subplot（3,2,1）;
stem（0:
length（y1）-1,y1）;
subplot（3,2,2）;
stem（0:
length（y2）-1,y2）;
subplot（3,2,3）;
stem（0:
length（y3）-1,y3）;
subplot（3,2,4）;
stem（0:
length（y1）-1,y1）;
subplot（3,2,5）;
stem（0:
length（y5）-1,y5）;
subplot（3,2,6）;
stem（0:
length（y6）-1,y6）;
%%%figure;
%%%stem（0:
15,y6（1:
16））;
%%%%%%%%%%%%%%%
y1=abs（fft（x1,16））;
y2=abs（fft（x2,16））;
y3=abs（fft（x3,16））;
y4=abs（fft（x4,16））;
y5=abs（fft（x5,16））;
y6=abs（fft（x6,16））;
figure;
subplot（3,2,1）;
stem（0:
length（y1）-1,y1）;
subplot（3,2,2）;
stem（0:
length（y2）-1,y2）;
subplot（3,2,3）;
stem（0:
length（y3）-1,y3）;
subplot（3,2,4）;
stem（0:
length（y6）-1,y4）;
subplot（3,2,5）;
stem（0:
length（y6）-1,y5）;
subplot（3,2,6）;
stem（0:
length（y6）-1,y6）;
%figure;
%stem（0:
15,y6（1:
16））;
x7=x4+j*x5;
y71=abs（fft（x7,8））;
y72=abs（fft（x7,16））;
figure;
subplot（1,2,1）;
stem（0:
length（y71）-1,y71）;
subplot（1,2,2）;
stem（0:
length（y72）-1,y72）;
程序结果如下所示：
每个信号时域分析
因为第六个信号不易在上图中展示（包含信息不到一个周期）所以单独画出图形如下：
x6的时域波形
逐个进行8点离散傅里叶变换
逐个进行16点离散傅里叶变换
x（n）8点和16点离散傅里叶变换比较
七、回答问题
（1）在N=8时，x2（n）和x3（n）的幅频特性会相同吗?
为什么?
N=16呢?
答：
N=8时幅频特性一样，N=16时幅频特性不一样。
Fn所能分辨到频率为为Fs/N，如果采样频率Fs为1024Hz，采样点数为1024点，则可以分辨到1Hz。
1024Hz的采样率采样1024点，刚好是1秒，也就是说，采样1秒时间的信号并做FFT，则结果可以分析到1Hz，如果采样2秒时间的信号并做FFT，则结果可以分析到0.5Hz。
如果要提高频率分辨力，则必须增加采样点数，也即采样时间。
频率分辨率和采样时间是倒数关系。
总结：
采样点N的不同，分辨率不同，在N=8，俩函数不同的部分没有分辨出来，因此特性相同，而N=16时已经足以分辨出不同的频率，因此频谱不同。
（2）如果周期信号的周期预先不知道，如何用FFT进行谱分析?
答：
设一个定长的m值,先取2m,看2m/m的误差是否大，如大的话再取4m,看4m/2m的误差是否大，如不大，4m（4倍的m值）则可近似原来点的谱分析。
八、总结
本实验主要掌握fft函数的用法，包括采样点的意义。
这个意义已经在七、回答问题中描述，不再赘述。
实验三：
用窗函数法设计FIR数字滤波器
一、实验目的
（1）掌握用窗函数法设计FIR数字滤波器的原理和方法。
（2）熟悉线性相位FIR数字滤波器特性。
（3）了解各种窗函数对滤波特性的影响。
二、实验内容及步骤
（1）复习用窗函数法设计FIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。
（2）编写程序。
①编写能产生矩型窗、升余弦窗、改进升余弦窗和二阶升余弦窗的窗函数子程序。
②编写主程序。
其中幅度特性要求用dB表示。
窗函数法设计滤波器主程序框图
三、实验内容
四、思考题
（1）如果给定通带截止频率和阻带截止频率以及阻带最小衰减，如何用窗函数法设计线性相位低通滤波器?
写出设计步骤。
（2）如果要求用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，试求理想带通的单位脉冲响应hd（n）。
五、参考资料
数字信号处理教程——Matlab释义与实现第二版，陈怀琛著，电子工业出版社，2008年
数字信号处理第三版，高西泉等著，西安电子科技大学出版社，第三版
XX知道
六、实验程序及结果
myfreqz（b,a）函数，是自编的改进freqz（）的函数。
function[db,mag,pha,grd,w]=myfreqz（b,a）;
%advancedfreqz
%[db,mag,pha,grd,w]=myfreqz（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）;
以下是主程序部分
wp=0.2*pi;ws=0.3*pi;deltaw=ws-wp;
N0=ceil（6.6*pi/deltaw）;
N=N0+mod（N0+1,2）;
wdham=（hamming（N））';
wc=（ws+wp）/2;
%hd=ideallp（wc,N）;
tao=（N-1）/2;
n=[0:
（N-1）];
m=n-tao+eps;
hd=sin（wc*m）./（pi*m）;
h=hd.*wdham;
[db,mag,pha,grd,w]=myfreqz（h,[1]）;
dw=2*pi/1000;
Rp=-（min（db（1:
wp/dw+1）））;
As=-round（max（db（ws/dw+1:
501）））;
figure;
subplot（2,2,1）;
stem（0:
length（hd）-1,hd,'.'）;
subplot（2,2,2）;
stem（0:
length（wdham）-1,wdham,'.'）;
subplot（2,2,3）;
stem（0:
length（h）-1,h,'.'）;
subplot（2,2,4）;
plot（1/length（db）:
1/length（db）:
1,db）;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wn1=boxcar（N）;
figure;
subplot（2,2,1）;
stem（0:
N-1,wn1,'.'）;
title（'boxcar'）;
wn2=bartlett（N）;
subplot（2,2,2）;
stem（0:
N-1,wn2,'.'）;title（'bartlett'）;
wn3=hanning（N）;
subplot（2,2,3）;
stem（0:
N-1,wn3,'.'）;title（'hanning'）;
%wn4=hamming（N）;
wn4=blackman（N）;
subplot（2,2,4）;
stem（0:
N-1,wn4,'.'）;title（'blackman'）;
%wn=（N,beta）;
得到的结果如下图：
其中，图一是理想脉冲响应，第一行图二是汉明窗；图三是实际脉冲响应，图四是幅度响应。
一下是几个典型窗函数的图像。
，汉明窗上例已画出，不再画。
每个窗的名字已在每个子图上方标出。
典型窗函数
七、回答问题
（1）给定通带截止频率和阻带截止频率以及阻带最小衰减，用窗函数法设计线性相位低通滤波器的设计步骤：
答：
技术指标Wp=0.2*pi,Ws=0.4*pi,Ap=0.25dB,As=50dB
方法一
选择海明窗
clearall;
Wp=0.2*pi;
Ws=0.4*pi;
tr_wide=Ws-Wp;%过渡带宽度
N=ceil（6.6*pi/tr_wide）+1;%滤波器长度
n=0:
1:
N-1;
Wc=（Wp+Ws）/2;%理想低通滤波器的截止频率
hd=ideal_lp1（Wc,N）;%理想滤波器的单位冲击响应
w_ham=（hamming（N））';%海明窗
h=hd.*w_ham;%实际海明窗的响应
[db,mag,pha,w]=freqz_m2（h,[1]）;%计算实际滤波器的幅度响应
delta_w=2*pi/1000;
Ap=-（min（db（1:
1:
Wp/delta_w+1）））%实际通带纹波
As=-round（max（db（Ws/delta_w+1:
1:
501）））%实际阻带纹波
subplot（221）
stem（n,hd）
title（'理想单位脉冲响应hd（n）'）
subplot（222）
stem（n,w_ham）
title（'海明窗'）
subplot（223）
stem（n,h）
title（'实际单位脉冲响应hd（n）'）
subplot（224）
plot（wi/pi,db）
title（'幅度响应（dB）'）
axis（[0,1,-100,10]）
方法二
Window=blackman（16）;
b=fir1（15,0.3*pi,'low',Window）;
freqz（b,128）
（2）用窗函数法设计带通滤波器，且给定上、下边带截止频率为ω1和ω2，理想带通的单位脉冲响应hd（n）的求解过程：
答：
由hd（n）=1/2π∫w1woe-jwaejnwdw+1/2π∫w2w3e-jwaejnwdw
（其中w0=-w0-wc,w1=-w0+wc,w2=w0-wc,w3=w0+wc）
计算整理后可得：
hd（n）=2/（（n-a）*π）*sin[（n-a）wc]*cos[（n-a）w0]
=2wc/π*sa[（n-a）wc]*cos[（n-a）w0]
实验四：
用双线性变换法设计IIR数字滤波器
一、实验目的
（1）掌握用双线性变换法设计IIR数字滤波器的原理和方法。
（2）熟悉IIR数字滤波器特性。
二、实验内容及步骤
（1）复习用双线性变换法设计IIR数字滤波器一节内容，阅读本实验原理，掌握设计步骤。
（2）编写程序。
①设计产生巴特沃斯滤波器原型，采用双线性变换法设计巴特沃斯数字低通滤波器。
②编写主程序。
其中幅度特性要求用dB表示。
三、上机实验内容
（1）实验所需设计的IIR数字低通滤波器的指标为wp=0.2*pi，ws=0.3*pi，Rp=1dB，As=15dB，滤波器
升级会员 