五邑大学数字信号处理34章实验报告.docx
《五邑大学数字信号处理34章实验报告.docx》由会员分享，可在线阅读，更多相关《五邑大学数字信号处理34章实验报告.docx（17页珍藏版）》请在冰豆网上搜索。
五邑大学数字信号处理34章实验报告
五邑大学
数字信号处理课程设计报告
题目：
FFT和数字巴特沃斯滤波器的设计
院系
专业
姓名学号
指导教师
实验第3章matlab实现FFT应用
一、实验目的
熟悉matlab在数字信号处理中的应用，并掌握基二FFT算法的实现。
二、实验原理
FFT快速算法，matlab fft（x,N）实现应用， 在MATLAB中，可以用函数X=fft（x，N）和x=ifft（X，N）计算N点序列的DFT正、反变换。
N点序列的DFT和IDFT变换定义式如下：
，
利用旋转因子
具有周期性，可以得到快速算法（FFT）。
三、程序代码与实验波形
%计算N点的DFT
function[Xk]=dft（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
End
%计算N点的逆DFT
function[xn]=idft（Xk,N）
n=[0:
1:
N-1];
k=[0:
1:
N-1];
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^（-nk）;
xn=（Xk*WNnk）/N;
end
%循环移位函数
function[y]=cirshift（x,m,N）
iflength（x）>N
error（'N必须大于等于x（n）的长度'）
end
x=[xzeros（1,N-length（x））];
n=0:
N-1;
n=mod（n-m,N）;
y=x（n+1）;
%计算序列x1（n）和x2（n）之间的N点循环卷积函数
function[y]=circonvt（x1,x2,N）
iflength（x1）>N
error（'N必须大于等于x1的长度'）
end
iflength（x2）>N
error（'N必须大于等于x2的长度'）
end
x1=[x1zeros（1,N-length（x1））];
x2=[x2zeros（1,N-length（x2））];
m=[0:
1:
N-1];
x2=x2（mod（-m,N）+1）;
H=zeros（N,N）;
forn=1:
1:
N
H（n,:
）=cirshift（x2,n-1,N）;
end
y=x1*H';
例3.2（P139）
>>x=[1,1,1,1,1,1];n=0:
5;
>>k=-200:
200;w=（pi/100）*k;
>>X=x*（exp（-j*pi/100））.^（n'*k）;
>>magX=abs（X）;angX=angle（X）*180/pi;
>>figure
（1）
>>subplot（2,1,1）;plot（w/pi,magX）;
>>axis（[-1106]）;grid;
>>title（'DTFT的幅度'）;xlabel（'以\pi为单位的频率'）;ylabel（'幅度'）;
>>subplot（2,1,2）;plot（w/pi,angX）;grid;
>>axis（[-11-200200]）;title（'DTFT的相位'）;
>>xlabel（'以\pi为单位的频率'）;ylabel（'相位'）;
>>N=6;
>>X=dft（x,N）;
>>magX=abs（X）;phaX=angle（X）*180/pi;
>>k=0:
5;
>>figure
（2）
>>subplot（2,1,1）;stem（k,magX）;
>>title（'DFT的幅度'）;xlabel（'k'）;
>>subplot（2,1,2）;stem（k,phaX）;
>>title（'DFT的相位'）;xlabel（'k'）;
Xk=
Columns1through5
6.0000-0.0000-0.0000i-0.0000-0.0000i0-0.0000i0.0000-0.0000i
Column6
0.0000-0.0000i
例3.3（P141）
>>n=0:
9;x=0.8*cos（0.47*pi*n）+0.4*cos（0.53*pi*n）;
>>N=10;
>>X=dft（x,N）;
>>magX=abs（X）,angX=angle（X）
>>k=0:
9;
>>subplot（2,1,1）;stem（k,magX）;title（'DFT的幅度'）;
>>xlabel（'k'）;
>>subplot（2,1,2）;stem（k,angX）;title（'DFT的相位'）;
>>xlabel（'k'）;
Xk=
Columns1through5
1.16111.2250+0.7289i2.0522+3.5598i0.9030-2.9037i0.8213-0.6744i
Columns6through10
0.8360-0.0000i0.8213+0.6744i0.9030+2.9037i2.0522-3.5598i1.2250-0.7289i
magX=
1.16111.42544.10903.04091.06270.83601.06273.04094.10901.4254
angX=
00.53671.0478-1.2693-0.6875-0.00000.68751.2693-1.0478-0.5367
例3.4（P143）
>>n=[0:
10];x=8*（0.8）.^n;N=11;
>>y=cirshift（x,6,N）;
>>subplot（2,1,1）;stem（n,x）;title（'序列x（n）'）;
>>xlabel（'n'）;ylabel（'x（n）'）;
>>subplot（2,1,2）;stem（n,y）;title（'x（n）的循环移位'）;
>>xlabel（'n'）;ylabel（'y（n）'）;
例3.5（P144）
>>x1=[1,2,2];x2=[5,4,3,2,1];N=5;
>>y=circonvt（x1,x2,N）
y=
1116211611
例3.6（P144）
>>n=0:
9;x1=（0.8）.^n;
>>n=0:
6;x2=exp（-n）;N=10;
>>y=circonvt（x1,x2,N）
y=
1.09051.20081.08150.90960.74400.60120.48320.38660.30930.2474
例3.7（P146）
>>N=32;fs=100;T=1/fs;
>>n=0:
N-1;r=n*T;x=r.*exp（r）;
>>X=fft（x,N）;
>>magX=abs（X）,phaX=angle（X）
>>subplot（3,1,1）;stem（n,x）;axis（[03201]）;
>>xlabel（'n'）;title（'序列x（n）'）;
>>k=0:
N-1;
>>subplot（3,1,2）;stem（k,magX）;axis（[03208]）;
>>xlabel（'k'）;ylabel（'DFT的幅度'）;
>>subplot（3,1,3）;stem（k,phaX）;axis（[032-44]）;
>>xlabel（'k'）;ylabel（'DFT的相位'）;
magX=
Columns1through10
6.13572.22151.11820.75180.57040.46300.39290.34410.30870.2824
Columns11through20
0.26250.24750.23630.22810.22260.21940.21830.21940.22260.2281
Columns21through30
0.23630.24750.26250.28240.30870.34410.39290.46300.57040.7518
Columns31through32
1.11822.2215
phaX=
Columns1through10
01.57471.72041.83471.94112.04432.14592.24672.34692.4467
Columns11through20
2.54632.64572.74502.84422.94343.04253.1416-3.0425-2.9434-2.8442
Columns21through30
-2.7450-2.6457-2.5463-2.4467-2.3469-2.2467-2.1459-2.0443-1.9411-1.8347
Columns31through32
-1.7204-1.5747
例3.8（P148）
>>N=21;L=256;
>>f1=120;f2=140;fs=400;
>>T=1/fs;ws=2*pi*fs;
>>n=0:
N-1;
>>x=cos（2*pi*f1*n*T）+cos（2*pi*f2*n*T）;
>>X=fftshift（fft（x,L））;
>>w=（-ws/2+（0:
L-1）*ws/L）/（2*pi）;
>>subplot（2,1,1）;plot（w,abs（X））;
>>ylabel（'幅度谱'）;xlabel（'频率（Hz）N=21'）;axis（[-200200015]）;
>>N=11;n=0:
N-1;
>>x=cos（2*pi*f1*n*T）+cos（2*pi*f2*n*T）;
>>X=fftshift（fft（x,L））;
>>subplot（2,1,2）;plot（w,abs（X））;axis（[-200200025]）;
>>ylabel（'幅度谱'）;xlabel（'频率（Hz）N=11'）;
例3.10（P152）
>>n=0:
20;
>>x=cos（0.1*pi*n）;
>>h=（0.8）.^n;
>>L=length（x）+length（h）-1;
>>X=fft（x,L）;
>>H=fft（h,L）;
>>y=ifft（X.*H）
>>subplot（3,1,1）;stem（n,x）;ylabel（'序列x（n）'）;xlabel（'n'）;
>>subplot（3,1,2）;stem（n,h）;ylabel（'序列h（n）'）;xlabel（'n'）;
>>n=0:
L-1;
>>subplot（3,1,3）;stem（n,real（y））;
>>ylabel（'卷积结果y（n）'）;xlabel（'n'）;
y=
Columns1through10
1.00001.75112.20992.35572.19361.75481.09490.2881-0.5785-1.4139
Columns11through20
-2.1311-2.6559-2.9338-2.9348-2.6569-2.1255-1.3914-0.52530.38881.2621
Columns21through30
2.00971.59851.27001.00860.80140.63830.51060.41140.33450.2751
Columns31through40
0.22880.19230.16260.13750.11550.09520.07620.05810.04100.0254
Column41
0.0115
例3.11（P153）
>>b=[0.0181,0.0543,0.0543,0.0181];
>>a=[1,-1.7600,1.1829,-0.2781];
>>ch=impseq（0,0,20）;n=0:
20;
>>h=filter（b,a,ch）;
>>x=cos（0.1*pi*n）+0.32*randn（size（n））;
>>L=length（x）+length（h）-1;
>>X=fft（x,L）;
>>H=fft（h,L）;
>>y=ifft（X.*H）
>>subplot（3,1,1）;plot（n,x）;axis（[020-1.51.5]）;
>>ylabel（'信号X（n）'）;xlabel（'n'）;
>>subplot（3,1,2）;stem（n,h）;axis（[020-0.50.5]）;
>>ylabel（'系统冲激响应h（n）'）;xlabel（'n'）;
>>n=0:
L-1;
>>subplot（3,1,3）;plot（n,y）;axis（[020-1.51.5]）;
>>ylabel（'x（n）滤波的结果y（n）'）;xlabel（'n'）;\
y=
Columns1through10
0.01300.07950.22830.42460.59090.65910.60500.45560.2503-0.0104
Columns11through20
-0.3476-0.7126-0.9934-1.1208-1.0877-0.9072-0.6222-0.3122-0.03860.2081
Columns21through30
0.48250.77170.94820.91290.70010.41540.1572-0.0203-0.1061-0.1199
Columns31through40
-0.0922-0.0496-0.01130.01210.02170.02040.01440.00630.0008-0.0023
Column41
-0.0013
四、实验分析
根据实验绘出的图形，与实际运算出的结果相比较，可知，利用ifft函数与fft函数求出的值与实际求出的IDFT变换与DFT变换值的误差相差较小。
五、实验结论
如果采样数据过少，运用Fourier变换不能分辨其中的频率成分。
添加零后可增加频谱的数据个数，谱的密度增高了，但仍不能分辨，其中的频率成分，即谱的分辨率没有提高。
只有数据点数足够多时才能分辨其中的频率成分。
实验第4章matlab实现数字巴特沃斯滤波器设计
一、实验目的
1、掌握设计IIR滤波器的原理和具体设计方法；
2、会在计算机上用IIR滤波器进行数字滤波；
3、掌握对IIR滤波器特性分析。
二、实验原理
IIR数字滤波器相关知识。
IIR数字滤波器的设计一般是利用目前已经很成熟的模拟滤波器的设计方法来进行设计，通常采用模拟滤波器原型有butterworth函数、chebyshev函数、bessel函数、椭圆滤波器函数等。
三、程序代码与实验波形
%butterworth低通滤波器原型设计函数，要求Ws﹥Wp﹥0，As﹥Rp﹥0。
function[b,a]=afd_butt（Wp,Ws,Rp,As）
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（Wp/Ws）））;
%ceil朝正无穷大方向取整；
fprintf（'\nButterworthFilterOrder=%2.0f\n',N）
OmegaC=Wp/（（10^（Rp/10）-1）^（1/（2*N）））%求对应于N的3db截止频率；
[b,a]=u_buttap（N,OmegaC）;
%非归一化Butterworth模拟低通滤波器原形设计函数
%得到的b,a分别为传输函数分子、分母多项式系数；
function[b,a]=u_buttap（N,Omegac）;
[z,p,k]=buttap（N）;%归一化巴特沃思模拟低通滤波器原形
p=p*Omegac;%将代入上式，相当于分子乘以，极点乘以
k=k*Omegac^N;
B=real（poly（z））;%poly为构造具有指定根的多项式real为求实部
b=k*B;
a=real（poly（p））;
%利用脉冲响应不变法从模拟到数字滤波器变换函数
function[b,a]=imp_invr（c,d,T）
[R,p,k]=residue（c,d）;%部分分式展开
p=exp（p*T）;%从模拟到数字极点对应关系,部分分式系数相同
[b,a]=residuez（R,p,k）;%将部分分式的形式变换成多项式之比的形式
b=real（b'）;%求出数字滤波器系数
a=real（a'）;
%频率响应函数freqz的修正，此函数可获得滤波器的幅值响应、相位响应及群延迟响应
function[db,mag,pha,w]=freqz_m（b,a）
[H,w]=freqz（b,a,1000,'whole'）;%在0-2*pi之间选取N个点计算频率响应
H=（H（1:
501））';%频率响应
w=（w（1:
501））';%频率
mag=abs（H）;%响应幅度
db=20*log10（（mag+eps）/max（mag））;%增益
pha=angle（H）;%相位
%变直接形式为级联形式
function[b0,B,A]=dir2cas（b,a）
b0=b
（1）;b=b/b0;a0=a
（1）;a=a/a0;b0=b0/a0;%以上步骤求出系数
M=length（b）;N=length（a）;
ifN>M
b=[bzeros（1,N-M）];
elseifM>N
a=[azeros（1,M-N）];
else
NM=0;
end
K=floor（N/2）;B=zeros（K,3）;A=zeros（K,3）;
ifK*2==N
b=[b0];a=[a0];
end
broots=cplxpair（roots（b））;%以下程序将每两个极点和两个零点组合成二阶因子
aroots=cplxpair（roots（a））;%roots：
求多项式的根
fori=1:
2:
2*K
Brow=broots（i:
1:
i+1,:
）;
Brow=real（poly（Brow））;
B（fix（i+1）/2,:
）=Brow;
Arow=aroots（i:
1:
i+1,:
）;
Arow=real（poly（Arow））;
A（fix（i+1）/2,:
）=Arow;
end
%此函数产生理想低通滤波器的冲激响应
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;Rp=1;As=15;T=1;
OmegaP=wp/T;OmegaS=ws/T;
[cs,ds]=afd_butt（OmegaP,OmegaS,Rp,As）;
[b,a]=imp_invr（cs,ds,T）
[db,mag,pha,w]=freqz_m（b,a）;
subplot（2,1,1）;plot（w/pi,mag）;
title（'digitalfilterMagnitudeResponse'）
axis（[0,1,0,1.1]）
subplot（2,1,2）;plot（w/pi,db）;
title（'digitalfilterMagnitudeinDB'）
axis（[0,1,-40,5]）;
五、实验分析
本次数字滤波器设计方法是基于MATLAB的数字滤波器的设计，是用学过的数字信号理论为依据，用MATLAB代码来实现，由滤波器的频谱图和滤波前后的信号的频谱图对比可知本设计选用双线性变换法设计的IIR滤波器比较好。
在同样的技术指标的要求下，IIR滤波器所要求的阶数N也比较小，实现起来比较容易。
六、实验结论
在信号与信息的过滤、检测和预测等处理中，都要使用滤波器，数字滤波器是数字信号处理中使用最广泛的一种方法。
IIR数字滤波器的设计过程中，可以借助模拟滤波器的设计成果或直接采用典型的滤波器类型，减少工作量.。
通过这个实验，对设计数字滤波器的整个过程有了很好的掌握。
升级会员 