数字信号处理实验指示书.docx
《数字信号处理实验指示书.docx》由会员分享，可在线阅读，更多相关《数字信号处理实验指示书.docx（10页珍藏版）》请在冰豆网上搜索。
数字信号处理实验指示书
数字信号处理
实验指示书
西南交通大学信息科学与技术学院
张翠芳王小川尹忠科
实验一序列的傅立叶变换
实验目的
进一步加深理解DFS,DFT算法的原理;
1.研究补零问题;
2.快速傅立叶变换（FFT）的应用。
实验步骤
1.复习DFS和DFT的定义，性质和应用；
2.熟悉MATLAB语言的命令窗口、编程窗口和图形窗口的使用；
3.利用提供的程序例子编写实验用程序；
4.按实验内容上机实验，并进行实验结果分析；
5.写出完整的实验报告，并将程序附在后面。
实验内容
1.周期方波序列的频谱
已知
试画出下面四种情况下的的幅度频谱,并分析补零后，对信号频谱的影响。
2.有限长序列x（n）的DFT
已知：
（1）取x（n）（n=0:
10）时，画出x（n）的频谱X（k）的幅度；
（2）将
（1）中的x（n）以补零的方式，使x（n）加长到（n:
0~100）时，画出x（n）的频谱X（k）的幅度；
（3）
（3）取x（n）（n:
0~100）时，画出x（n）的频谱X（k）的幅度。
3.利用FFT进行谱分析
已知：
模拟信号
以t=0.01n（n=0:
N-1）进行采样，求N点DFT的幅值谱。
请分别画出N=45;N=50;N=55;N=60时的幅值曲线。
数字信号处理实验一（MATLAB语言程序示例）
%Example1
L=5;N=20;
n=1:
N;
xn=[ones（1,L）,zeros（1,N-L）];
Xk=dfs（xn,N）;
magXk=abs（[Xk（N/2+1:
N）Xk（1:
N/2+1）]）;
k=[-N/2:
N/2];
figure
（1）
subplot（2,1,1）;
stem（n,xn）;xlabel（'n'）;ylabel（'xtide（n）'）;
title（'DFSofSQ.wave:
L=16,N=64'）;
subplot（2,1,2）;
stem（k,magXk）;
axis（[-N/2,N/2,0,16]）;
xlabel（'k'）;ylabel（'Xtide（k）'）;
%Example2
M=100;
N=100;
n=1:
M;
xn=cos（0.48*pi*n）+cos（0.52*pi*n）;
n1=[0:
1:
N-1];y1=[xn（1:
1:
M）,zeros（1,N-M）];
figure
（1）
subplot（2,1,1）;
stem（n1,y1）;xlabel（'n'）;
title（'signalx（n）,0<=n<=100'）;
axis（[0,N,-2.5,2.5]）;
Y1=fft（y1）;
magY1=abs（Y1（1:
1:
N/2+1））;
k1=0:
1:
N/2;
w1=2*pi/N*k1;
subplot（2,1,2）;
title（'SamplesofDTFTMagnitude'）;
stem（w1/pi,magY1）;
axis（[0,1,0,10]）;
xlabel（'frequencyinpiunits'）;
%example3
figure
（1）
subplot（2,2,1）
N=45;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=fft（x,N）;
plot（q,abs（y））
stem（q,abs（y））
title（'FFTN=45'）
%
subplot（2,2,2）
N=50;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=fft（x,N）;
plot（q,abs（y））
title（'FFTN=50'）
%
subplot（2,2,3）
N=55;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=fft（x,N）;
plot（q,abs（y））
title（'FFTN=55'）
%
subplot（2,2,4）
N=16;n=0:
N-1;t=0.01*n;
q=n*2*pi/N;
x=2*sin（4*pi*t）+5*cos（8*pi*t）;
y=fft（x,N）;
plot（q,abs（y））
title（'FFTN=16'）
function[Xk]=dfs（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;
实验二用双线性变换法设计IIR数字滤波器
一、实验目的
1．熟悉用双线性变换法设计IIR数字滤波器的原理与方法；
2．掌握数字滤波器的计算机仿真方法；
3．通过观察对实际心电图的滤波作用，获得数字滤波器的感性知识。
二、实验内容
1．用双线性变换法设计一个巴特沃斯低通IIR滤波器，设计指标参数为：
在通带内频率低于0.2π时，最大衰减小于1dB；在阻带内[0.3π，π]频率区间上，最小衰减大于15dB。
2．以0.2π为采样间隔，打印出数字滤波器在频率区间[0，0.2π]上的幅值响应曲线。
3．用所设计的滤波器对实际的心电图信号采样序列x（n）（见教科书上254页实验二的第六部分）进行仿真滤波处理，并分别打印出滤波前后的心电图信号波形图，观察总结滤波作用与效果。
三、实验步骤
1．（参考数字信号处理教材实验二）；
2．（参考数字信号处理教材实验二）；
3．用MATLAB语言编写仿真滤波程序，完成实验内容1~3；
4．写出完整的实验报告并回答教科书上的思考题。
四、IIR滤波器设计的常用函数
1.Butter
功能：
Butterwoeth（巴特沃斯）数字滤波器设计。
格式：
[b，a]=butter（n，Wn）
可设计出截止频率为Wn的n阶低通Butterwoeth滤波器，其滤波器为
截止频率是滤波器幅度下降至
处的频率.
2．Buttord
IIR（巴特沃斯）滤波器阶的选择
格式
[n,Wn]=buttord（Wp,Ws,Rp,Rs）----------数字域
[n,Wn]=buttord（Wp,Ws,Rp,Rs,’s’）-----模拟域
说明：
buttord可在给定滤波器性能的情况下，选择模拟或数字Butterword滤波器的最小阶，其中Wp和Ws分别是通带和阻带的截止频率，其值
，当其值为1时，表示0.5
Rp,RS分别是通带和阻带区的波纹系数。
例如：
设计一低通滤波器，通带范围0~100Hz，通带波纹小于3dB，阻带为-30dB，并利用最小的阶来实现，其程序如下
Wp=100/500；Ws=200/500；
[n，Wn]=buttord（Wp,Ws,3，30）；
[b，a]=butter（n，Wn）；
freqz（b，a，512，1000）
图低通滤波器特性（n=8）
实验三用窗函数法设计FIR数字滤波器
1.实验目的及原理方法参考（301~302）页
2.实验内容
●用MATLAB产生各种窗函数
●利用窗函数设计FIR滤波器
参考程序：
（1）用MATLAB实现各种窗函数
m=200;
a=boxcar（m）;矩形窗
b=hanning（m）;汉宁窗
c=hamming（m）;海明窗
d=blackman（m）;布莱克曼窗
e=kaiser（m,7.865）;凯泽窗
m=1:
200;
plot（m,a,'r*',m,b,'g+',m,c,'y*',m,d,'b.'）
（2）利用窗函数设计FIR滤波器
设计具有下列指标p=0.25,Rp=0.25dB，s=0.3,Rp=50dB的低通数字滤波器。
要求：
●选择合适的窗函数;
●画出滤波器的频率特性
几个主要的MATLAB函数
（1）[r,p,k]=residuez（b,a）将有理分式变换成部分分式表示
（2）[b,a]=residuez（r,p,k）将部分分式换成有理分式变表示
式中各参数意义如下：
（3）freqz（b,a）画连续系统的频率特性曲线
（4）impz（b,a，101）画系统的脉冲响应
例如：
N=5;Wn=[100200]/500;
[b,a]=butter（n,Wn）;
impz（b,a,101）
参考程序：
function[db,mag,pha,grd,w]=freqz_m（b,a）;
[H,w]=freqz（b,a,1000,'whole'）;
H=（H（1:
501））';w=（w（1:
501））';
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
pha=angle（H）;
grd=grpdelay（b,a,w）;
（GRPDELAYGroupdelayofadigitalfilter.
[Gd,W]=GRPDELAY（B,A,N）returnslengthNvectorsGdandWcontainingthegroupdelayandthefrequencies（inradians）atwhichitisevaluated.Groupdelayis-d{angle（w）}/dw.ThefrequencyresponseisevaluatedatNpointsequallyspacedaroundtheupperhalfoftheunitcircle.ForanFIRfilterwhereNisapoweroftwo,
thecomputationisdonefasterusingFFTs.Ifyoudon'tspecifyN,itdefaultsto512.）
GRPDELAY（B,A,N,'whole'）usesNpointsaroundthewholeunitcircle.
[Gd,F]=GRPDELAY（B,A,N,Fs）and[Gd,F]=GRPDELAY（B,A,N,'whole',Fs）
givensamplingfrequencyFsinHzreturnavectorFinHz.
Gd=GRPDELAY（B,A,W）andGd=GRPDELAY（B,A,F,Fs）returnthegroupdelayevaluatedatthepointsinW（inradians/sample）orF（inHz）.
GRPDELAY（B,A,...）withnooutputargumentsplotsthegroupdelayinthe
currentfigurewindow.
/求给定的理想低通滤波器的单位冲激响应。
functionhd=ideal_lp（wc,M）;
alpha=（M-1）/2;
n=[0:
（M-1）];
m=n-alpha+eps;
hd=sin（wc*m）./（pi*m）;
wp=0.2*pi;ws=0.3*pi;
tr_width=ws-wp;
M=ceil（6.6*pi/tr_width）+1;
n=[0:
1:
M-1];
wc=（ws+wp）/2;
hd=ideal_lp（wc,M）;
w_ham=（hamming（M））';
h=hd.*w_ham;
[db,mag,pha,grd,w]=freqz_m（h,[1]）;
delta_w=2*pi/1000;
Rp=-（min（db（1:
1:
wp/delta_w+1）））%PassbandRipple
As=-round（max（db（ws/delta_w+1:
1:
501）））%MinStopbandattenuation
%plots
figure
（1）
subplot（2,2,1）;stem（n,hd）;
title（'IdealImpulseresponse'）
axis（[0M-1-0.10.3]）;
ylabel（'hd（n）'）
subplot（2,2,2）;stem（n,w_ham）;
title（'HammingWindow'）
axis（[0M-101.1]）;ylabel（'w（n）'）
subplot（2,2,3）;stem（n,h）;
title（'ActualImpulseResponse'）
axis（[0M-1-0.10.3]）;
ylabel（'h（n）'）
subplot（2,2,4）;plot（w/pi,db）;
title（'MagnitudeResponseindB'）;
grid
axis（[01-10010]）;ylabel（'Decibels'）
升级会员 