北理工数字信号处理实验报告Word格式.docx
《北理工数字信号处理实验报告Word格式.docx》由会员分享，可在线阅读，更多相关《北理工数字信号处理实验报告Word格式.docx（34页珍藏版）》请在冰豆网上搜索。
WN=exp（-i*2*pi/Nz）;
forj=1:
Nz/2
fork=j:
Nz:
N
kp=k+Nz/2;
t=y（kp）*u;
%蝶形运算乘项
y（kp）=y（k）-t;
%加项
y（k）=y（k）+t;
end
u=u*WN;
end
end
y
y1=fft（x）%对比matlab的fft
y=
Columns1through3
36.0000-4.0000+9.6569i-4.0000+4.0000i
Columns4through6
-4.0000+1.6569i-4.0000-4.0000-1.6569i
Columns7through8
-4.0000-4.0000i-4.0000-9.6569i
y1=
2、编程实现序列长度为N=8的按频率抽取的基2-FFT算法。
L=2048;
M=256;
>
x=ones（1,L）;
n=0:
M-1;
h=cos（0.2*pi*n）;
tic
y1=conv（x,h）;
toc
tic
X=fft（x,L+M-1）;
H=fft（h,L+M-1）;
Y=X.*H;
y2=ifft（Y）;
运行结果如下：
（分别为线形卷积和快速卷积）
Elapsedtimeis0.003352seconds.
Elapsedtimeis0.007167seconds.
L=4096;
Elapsedtimeis0.003666seconds.
Elapsedtimeis0.030134seconds.
3.将上述FFT程序推广大序列长度为N=2v的情况，要求利用原位运算。
l=128;
x=ones（1,l）;
q=0;
H=fft（h,l+M-1）;
X=fft（x,l+M-1）;
y=ifft（Y）;
yp=y;
fori=2:
16
q=q+l;
yq=[zeros（1,q）,y];
yp=[yp,zeros（1,l）];
yp=yq+yp;
Elapsedtimeis0.018875seconds.
4.
n=0:
h=cos（0.2*pi*n）;
x=ones（1,L）;
Lx=2048;
M1=M-1;
N=512;
L=N-M1;
h=[h,zeros（1,N-M）];
x=[zeros（1,M1）,x,zeros（1,N-1）];
a=floor（（Lx+M1-1）/（L））+1;
Y=zeros（1,N）;
fork=0:
a-1
xk=x（k*L+1:
k*L+N）;
b=fft（xk,N）;
C=fft（h,N）;
Z=b.*C;
Y（k+1,:
）=ifft（Z,N）;
Y=Y（:
M:
N）'
;
Y=（Y（:
））'
Elapsedtimeis0.035329seconds.
实验2利用DFT分析信号频谱
一、实验目的
1、加深对DFT原理的理解。
2、应用DFT分析信号的频谱。
3、深刻理解利用DFT分析信号频谱的原理，分析实现过程中出现的现象及解决方法。
二、实验设备与环境
计算机、MATLAB软件环境
1、DFT与DTFT的关系
序列
的N点DFT
实际上就是
序列的DTFT在N个等间隔频率点
上的采样样本。
2、利用DFT求DTFT
方法一：
利用差值运算
其中
为内插函数
方法二：
在实际MATLAB计算中，采用增加数据的长度N，使得到的DFT谱线更加精细，用其包络近似计算DTFT。
如果数据量不足，可采用补零来增加数据长度。
3、利用DFT分析连续时间信号的频谱
1）确定时域采样间隔T，得到离散序列
2）确定截取长度M，得到M点离散序列
，这里
为窗函数
3）确定频域采样点数N，要求
4）利用FFT计算离散序列的N点DFT，得到
5）由
计算
采样点的近似值。
4、可能用到的MATLAB函数与代码
实验中DFT运算采用MATLAB中提供的函数fft来实现。
DTFT可以利用MATLAB矩阵运算的方法进行计算。
四、实验内容
1、
（1）
3;
x=[2-111];
w=-pi:
0.01*pi:
pi;
X=x*exp（-j*n'
*w）;
subplot（211）
plot（w,abs（X））;
axistight
subplot（212）
plot（w,angle（X）/pi）;
（2）
在原基础上增加代码：
x=[2-111];
N=4;
k=0:
Y=fft（x）;
w=0:
2*pi;
m=0:
2*pi/4:
3/2*pi
subplot（211）;
xlabel（'
\Omega/pi'
）;
title（'
DTFT-Magnitude'
holdon
stem（m,abs（Y））;
subplot（212）;
DTFT-Phase'
stem（m,angle（Y））
（3）
在原本基础上设置代码
Z=fft（x,64）;
%n=64DFT
N=0:
63;
stem（abs（Z））;
stem（angle（Z））
2、考察序列
，用DFT估计
的频谱；
将
补零加长到长度为100点序列用DFT估计
的频谱。
要求画出相应波形。
程序代码：
程序输出结果：
补零后：
时，用DFT估计
的频谱，并画出波形。
（3）根据实验结果，分析怎样提高频谱分辨率。
答：
减小截取长度或减小采样间隔的方法可以提高频谱分辨率。
3、已知信号
其中
=1Hz，
=2Hz，
=3Hz。
从
的表达式可以看出，它包含三个频率的正弦波，但是，从其时域波形来看，似乎是一个正弦信号，利用DFT做频谱分析，确定适合的参数，使得到的频谱的频率分辨率符合需要。
程序输出结果：
利用DFT近似分析连续时间信号
分析采用不同的采样间隔和截取长度进行计算的结果，并最终确定适合的参数。
通过改变步长可以达到改变采样间隔的结果，输出不同的频谱。
实验5、脉冲响应不变法设计IIR数字滤波器
1、掌握利用脉冲响应不变法设计IIR数字滤波器的原理及具体方法
2、加深理解数字滤波器和模拟滤波器之间的技术指标转化。
3、掌握脉冲相应不变法设计IIR数字滤波器的优缺点及适用范围。
1、基本原理：
从时域响应出发，使数字滤波器的单位脉冲响应h（n）模仿模拟滤波器的单位冲激响应。
2、设计步骤：
a.确定数字滤波器性能指标。
b.将数字滤波器频率指标转换成相应的模拟滤波器的频率指标
c.设计模拟滤波器
实现方法：
脉冲响应不变法
[r,p,k]=residue（b,a）
[b,a]=residue（r,p,k）
d.模拟滤波器数字化
[bz,az]=impinvar（b,a,fs）
1、设采样频率fs=4khz,采用脉冲响应不变法设计一个三阶巴特沃斯数字低通滤波器，其3db截止频率为fc=1khz.
首先设计巴特沃斯滤波器：
N=3;
Wc=0.25*pi;
[b,a]=butter（N,Wc,'
s'
）
b=
0000.4845
a=
1.00001.57081.23370.4845
因此，
然后脉冲响应不变法计算：
[bz,az]=impinvar（b,a）
bz=
1.1702-1.42120.6078-0.0945
az=
1.0000-1.62860.9286-0.1879
得到：
画出模拟滤波器和数字滤波器的图像进行对比：
[H,w]=freqs（b,a）;
subplot（221）;
plot（w/pi,abs（H））;
gridon;
subplot（222）;
plot（w/pi,20*log10（（abs（H））/max（abs（H））））;
w=[0:
500]*pi/500;
[H,w]=freqz（bz,az）;
subplot（223）;
subplot（224）;
2.设采样频率为fs=10khz，设计数字滤波器满足如下指标
通带截止频率：
1khz，通带波动：
1db，
阻带截止频率：
1.5khz，阻带衰减：
15db
分别采用巴特沃斯、切比雪夫I型、切比雪夫ii型和椭圆模拟原型滤波器及脉冲响应不变法进行设计。
结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求，分析脉冲响应不变法设计IIR数字滤波器的优缺点及适用范围。
（1）巴特沃斯滤波器为：
Wp=0.1*pi;
Ws=0.15*pi;
Rp=1;
As=15;
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（Wp/Ws）））
OC1=Wp/（（10^（Rp/10）-1）^（1/（2*N）））
[b,a]=butter（N,OC1,'
N=
6
OC1=
0.3516
0000000.0019
1.00001.35850.92270.39740.11410.02080.0019
1.0e-003*
0.00000.01250.25730.51940.16370.00510
1.0000-4.65179.1374-9.68295.8300-1.88890.2570
得到bz，az后，画出模拟滤波器和数字滤波器的图像进行对比：
（2）切比雪夫Ⅰ型
Wp=0.2*pi;
Ws=0.3*pi;
Rp=1;
As=15;
e=
e=（10^（Rp/10）-1）^（1/2）
e=0.5088
A=10^（As/20）
A=5.6234
N=ceil（acosh（（A*A-1）^（1/2）/e）/acosh（Ws/Wp））N=4
[b,a]=cheby1（N,Rp,Wp,'
b=00000.0383a=1.00000.59870.57400.18420.0430
（3）切比雪夫Ⅱ型：
N=4;
[b,a]=cheby2（N,As,Ws,'
b=0.1778-0.00001.2637-0.00001.1225a=1.00002.36963.03221.99251.1225
（4）椭圆模拟原型滤波器：
e=0.5088;
A=5.6234;
k1=e/（（A*A-1）^（1/2））
k1=0.0919
k=Wp/Ws
k=0.6667
[K,E]=ellipke（k）
K=2.0290
E=1.2612
[K1,E1]=ellipke（（1-k1*k1）^（1/2））
K1=4.1217
E1=1.0077
[K2,E2]=ellipke（k1）
K2=1.6089
E2=1.5340
[K3,E3]=ellipke（（1-k*k）^（1/2））
K3=2.1483
E3=1.2140
N=（K*K1）/（K2*K3）
N=2.4195
N2=2;
[b,a]=ellip（N2,Rp,As,Ws,'
b=0.177800.9379
a=1.00000.91201.0523
结果分析：
（1）当给定的设计指标相同时，选用椭圆滤波器所要求的结束N最低，切比雪夫滤波器
次之，巴特沃斯滤波器最高；
（2）从通带的相位相应来看，椭圆滤波器虽然提供了最优秀的幅度平方相应，但通带上
的相位相应非线性较大，而巴特沃斯滤波器在通带上具有相当的线性相位，切比雪
夫滤波器的相位特征介于两者之间。
实验七、窗函数法设计FIR数字滤波器
掌握窗函数法设计FIR数字滤波器的原理及具体方法。
二、实验设备及环境
计算机、matlab软件及环境
三、实验基础理论
1.基本原理
窗函数设计法的基本思想为，首先选择一个适当的理想的滤波器Hd（ejw），然后用窗函数截取它的单位脉冲响应hd（n），得到线性相位和因果的FIR滤波器。
这种方法的重点是选择一个合适的窗函数和理想滤波器，使设计的滤波器的单位响应尽力逼向滤波器的单位响应。
2.设计步骤
（1）给定理想滤波器的频率响应Hd（ejw），再通带上具有单位增益和线性相位，在阻带上具有零响应。
（2）确定这个滤波器的单位脉冲响应
为了得到一个h（n）长度为N的因果相位FIR滤波器，我们令
（3）用窗函数截取hd（n）得到所设计FIR数字滤波器h（n）
h（n）=hd（n）w（n）
四、实验内容
1.设计一个数字低通滤波器，其技术指标如下
0.2
分别采用矩形窗、汉宁窗、海明窗、布莱克曼窗、凯瑟窗设计该滤波器。
结合实验结果，分别讨论采用上述方法设计的数字滤波器是否都能满足给定指标要求。
2.设计一个数字低通滤波器，其技术指标如下
下阻带边缘：
下通带边缘：
0.35
上通带边缘：
0.65
上阻带边缘：
五、代码
确定长度：
wp=0.2*pi;
wst=0.3*pi;
width=wst-wp;
N=ceil（1.8*pi/width）+1
19
确定hd（n）
（N-1）;
wc=（wp+wst）/2;
alpha=（N-1）/2;
hd=（wc/pi）*sinc（（wc/pi）*（n-alpha））
wboxcar=boxcar（N）'
h=hd.*wboxcar
stem（n,hd,'
filled'
axistight;
n'
ylabel（'
hd（n）'
[Hr,w1]=zerophase（h）;
subplot（222）;
plot（w1/pi,Hr）;
axistight;
\omega^pi'
\omega'
subplot（223）;
stem（n,h,'
h（n）'
[H,w]=freqz（h,1）;
subplot（224）;
plot（w/pi,20*log10（abs（H）/max（abs（H））））;
xlabel（'
dB'
gridon
由上图可清楚看出：
加矩形窗不满足指标
用海明窗：
wst=0.3*pi;
width=wst-wp;
N=ceil（6.6*pi/width）+1
67
wc=（wp+wst）/2;
alpha=（N-1）/2;
hd=（wc/pi）*sinc（（wc/pi）*（n-alpha））;
whamming=hamming（N）'
h=hd.*whamming;
subplot（221）;
stem（n,hd,'
[Hr,w1]=zerophase（h）;
[H,w]=freqz（h,1）;
plot（w/pi,20*log10（abs（H）/max（abs（H））））;
gridon
加海明窗满足指标
汉宁窗：
代码：
wp=0.2*pi;
wst=0.3*pi;
tr_width=wst-wp;
N=ceil（6.2*pi/tr_width）+1;
n=0:
wc=（wp+wst）/2;
alpha=（N-1）/2;
hd=（wc/pi）*sinc（（wc/pi）*（n-alpha））;
w_hanning=hanning（N）'
h=hd.*w_hanning;
升级会员 