实验指导书 36章.docx
《实验指导书 36章.docx》由会员分享，可在线阅读，更多相关《实验指导书 36章.docx（22页珍藏版）》请在冰豆网上搜索。
实验指导书36章
实验三离散傅里叶变换
一、实验目的
1.学习编制离散傅里叶变换程序.
2.理解DFT在数字信号处理中的核心地位和作用.
二、实验内容
1.编制计算离散傅里叶变换程序.
2.用编制程序处理时间抽样信号.
3.根据实序列离散傅里叶变换的对称性,初步判定程序的正确性.
三、实验说明
1.离散傅里叶变换公式如下
构造离散傅里叶正、反变换函数的MATLAB实现程序如下，其中dft（xn,N）为离散傅里叶正变换，idft（Xk,N）为离散傅里叶反变换:
---------------------------------------------------------------------------------------------------------------
function[Xk]=dft（xn,N）
n=[0:
1:
N-1];
k=n;
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^nk;
Xk=xn*WNnk
function[xn]=idft（Xk,N）
n=[0:
1:
N-1];
k=n;
WN=exp（-j*2*pi/N）;
nk=n'*k;
WNnk=WN.^（-nk）;
xn=（Xk*WNnk）/N;
---------------------------------------------------------------------------------------------------------------
2.利用离散傅里叶变换函数求解序列傅里叶变换的MATLAB实现程序如下（N=8，
a=0.7）。
---------------------------------------------------------------------------------------------------------------
clear
N=8;
a=0.7;
n=[0:
7];
xn=a.^n;
Xk=dft（xn,N）;
subplot（3,1,1）
stem（n,xn,'.k'）;axis（[0,8,0,1.5]）
subplot（3,1,2）
stem（n,abs（Xk）,'.k'）;axis（[0,8,0,5]）
subplot（3,1,3）
stem（n,angle（Xk）,'.k'）;axis（[0,8,-1.5,1.5]）
3.下图所示为
的幅度特性与相位特性。
其中,N为总的抽样点数;T为抽样时间间隔.
在这种条件下分析结果中的
以
点左右对称,说明该程序基本正确,可以进行其它分析.
4.在实验中同学们会发现,抽样信号分析得到的
的幅度与连续傅里叶变换分析该信号幅度不相同,这是值得大家深入讨论的问题.
3.14已知序列
=（3,11,7,0,-1,4,2）,令
为
加入噪声干扰并移位后的序列：
其中,
为具有零均值和单位方差的高斯序列。
用MATLAB语言计算
和
之间的互相关。
3.14解：
程序及运行结果如下。
Clear
n=[1:
7];
xn=zeros（1,length（n））;
xn2=xn;wn=xn;
xn（1:
7）=[3,11,7,0,-1,4,2];%signal
xn2（1:
7）=[4,2,3,11,7,0,-1];%delaysignal
wn（1:
7）=randn（1,7）;%noise
yn=xn+wn;%signal+noise
xyn=conv（xn,yn）;
subplot（2,2,1）;
stem（n,xn,'k'）;title（'x（n）'）;
subplot（2,2,2）;
stem（n,xn2,'k'）;title（'x（n-2）'）;
subplot（2,2,3）;
stem（n,wn,'k'）;title（'w（n）'）;
subplot（2,2,4）;
stem（1:
（2*length（n）-1）,xyn,'k'）;title（'xy（n）'）;
四、实验报告要求
1.整理好经过运行并证明是正确的程序,并且加上详细的注释.
2.用连续傅里叶变换分析（被抽样的）连续信号,将其结果与抽样信号的离散傅里叶变换结果相比较,你能发现什么问题?
如何解释?
3.计算抽样序列的连续傅里叶变换,将其结果与抽样序列的离散傅里叶变换结果相比较,你又能发现什么问题?
如何解释?
五、实验结果分析
实验四快速傅里叶变换
一、实验目的
1.学习时间抽选奇偶分解FFT算法.
2.深入理解和掌握时间抽选奇偶分解FFT计算程序.
3.研究如何利用FFT程序分析确定性时间连续信号.
二、实验内容
1.用MATLAB编程比较DFT和FFT的运算时间
2.用MATLAB编程实现DFT和FFT的运算.
已知有限长序列
长度为N＝4,
且:
1
2
–1
3
用FFT求X（k），再用IFFT求x（n）.
3.用MATLAB程序分析FFT取不同长度时，序列x（n）的频谱变化情况.
三、实验说明
1.可以用以下MATLAB程序比较DFT和FFT的运算时间：
---------------------------------------------------------------------------------------------------------------
N=1024;
M=80;
x=[1:
M,zeros（1,N-M）];
t=cputime;
y1=fft（x,N）;
Time_fft=cputime-t;
t1=cputime;
y2=dft（x,N）;
Time_dft=cputime-t1;
t2=cputime;
---------------------------------------------------------------------------------------------------------------
Time_dft=
6.0290
Time_fft=
0.0100
由此可见FFT算法比直接计算DFT速度快得多。
2.利用快速傅里叶变换函数求解FFT和IFFT运算的MATLAB实现程序如下：
---------------------------------------------------------------------------------------------------------------
clear
xn=[1,2,-1,3];
X=fft（xn）
x=ifft（X）
---------------------------------------------------------------------------------------------------------------
X=[5.0000,2.0000+1.0000i,-5.0000,2.0000-1.0000i]
x=[1,2,-1,3]
3.设x（n）是长度为N=6的矩形序列，用MATLAB分析FFT取不同长度时x（n）的频谱变化。
N=8，32，64时x（n）的FFTMATLAB实现程序如下
---------------------------------------------------------------------------------------------------------------
x=[1,1,1,1,1,1];
N=8;y1=fft（x,N）;
n=0:
N-1;subplot（3,1,1）;stem（n,abs（y1）,'.k'）;axis（[0,9,0,6]）;
N=32;y2=fft（x,N）;
n=0:
N-1;subplot（3,1,2）;stem（n,abs（y2）,'.k'）;axis（[0,40,0,6]）;
N=64;y3=fft（x,N）;
n=0:
N-1;subplot（3,1,3）;stem（n,abs（y3）,'.k'）;axis（[0,80,0,6]）;
---------------------------------------------------------------------------------------------------------------
x（n）的频谱如下图所示,N值取的越大就越接近序列真正的频谱。
4.8设
是长度为ML的长序列，其中
把
分成M段，记为
，
每段长度为L。
设
为L点冲激响应，则
显然，
是
点序列。
在这种方法中，需要保存中间卷积结果，在相加之前进行恰当的重叠，形成
。
（1）利用循环卷积，开发一个MATLAB函数实现重叠相加法；
（2）利用
（1）开发的函数采用基-2FFT，编写一个高速重叠相加分段卷积的MATLAB程序。
4.8解：
程序及运行结果如下。
（1）循环卷积子函数：
functionfn=circonvt（x1，x2，N）；
%circonvt函数实现输入序列x1和x2的循环倦积，fn为输出序列
%N为循环卷积长度
if（length（x1）>N|length（x2>N）%判断输入信号的长度
error（'N的长度必须大于输入数据的长度'）;
end
x1=[x1,zeros（1,N-length（x1））];
x2=[x2,zeros（1,N-length（x2））];
m=0:
N-1;
x=zeros（N,N）;
forn=0:
N-1
x（:
n+1）=x2（mod（（n-m）,N）+1）';
end;
fn=x1*x;
主函数:
function[y]=ovrlpadd（x,h,L）
x=input（'请输入x序列：
'）;
h=input（'请输入y序列：
'）;
L=input（'请输入段长L：
'）;
lenx=length（x）；%x的长度
M=length（h）；%h的长度
N1=L+M-1；%圆周卷积点数，即每一个输出序列Yi的长度
m=rem（lenx,L）;%求余
ifm=0
x=[xzeros（1,L-m）];%末尾补零，使每段长度为N
K=floor（lenx/L）+1;%段数
else
x=x;
K=floor（lenx/L）;
end
ytemp=zeros（1,N1-L）;%N1-N为重叠部分，使其初始化为零
n1=1;n2=L;
fork=1;K
xk=x（n1:
n2）;
Y（k,:
）=circonvt（xk,h,N1）;
fori=1:
N1-L
Y（k,i）=Y（k,i）+ytemp（i）;
ytemp（i）=Y（k,i+L）;
end
y（n1:
n2+M-1）=Y（k,1:
N1）;
n1=n1+L;n2=n2+L;
end
stem（y）;
title（'ovrlpadd'）;
请输入x序列：
[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
请输入y序列：
[1,0,0,-1]
请输入段长L:
4
ans=
1233333333333
333-14-15-16
（2）
function[y]=ovrlpaddfft（x,h,L）
x=input（'请输入x序列：
'）;
h=input（'请输入y序列：
'）;
L=input（'请输入段长L：
'）;
lenx=length（x）；%x的长度
M=length（h）；%h的长度
N1=L+M-1；%满足循环卷积等于线性卷积的长度
N1=2^（ceil（log10（N1）/log10
（2）））;
m=rem（lenx,L）;
ifm=0
x=[xzeros（1,L-m）];%末尾补零，使每段长度为N
K=floor（lenx/L）+1;%段数
else
x=x;
K=floor（lenx/L）;%段数
end
ytemp=zeros（1,N1-L）;%N1-N为重叠部分，使其初始化为零
n1=1;
n2=L;
fork=1;K
z=x（n1:
n2）
xk=fft（z,N1）;
Y（k,:
）=real（ifft（xk.*h））;
fori=1:
N1-L
Y（k,i）=Y（k,i）+ytemp（i）;
ytemp（i）=Y（k,i+L）;
end
y（n1:
n2+N1-L）=Y（k,1:
N1）;%输出结果
n1=n1+L;
n2=n2+L;
end
y=y（1:
（lenx+M-1））;
stem（y）;
title（'ovrlpadd-2fft'）;
请输入x序列：
[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16]
请输入y序列：
[1,0,0,-1]
请输入段长L:
4
ans=
Columns1through12
1.00002.00003.00003.00003.00003.00003.00003.0000
3.00003.00003.00003.0000
Columns13through19
3.00003.00003.00003.0000-14.0000-15.0000-16.0000
四、实验报告要求
1.整理好经过运行并证明是正确的程序,并且加上详细的注释.
2.把N=1024点FFT分析所用的时间与直接计算DFT所用的时间相对比,从感性上理解”快”字.
3.深入讨论如何运用FFT来研究序列的频谱,为第八章的学习做好准备.
五、实验结果分析
实验五IIR滤波器设计
一.实验目的
1.学习模拟滤波器的设计方法.
2.学习模拟---数字变换滤波器设计方法.
3.掌握双线性变换数字滤波器设计法.
4.掌握用频带变换设计数字滤波器的具体方法.
二.实验内容
1.设计一个巴特沃斯模拟低通滤波器，满足以下性能指标：
通带的截止频率
，通带最大衰减
，阻带的截止频率
，阻带最大衰减
。
2.已知模拟滤波器的系统函数为
分别用冲激响应不变变焕法和双线性变换法将
转换为数字滤波器的系统函数
，并画出相应的
和
的频率响应曲线。
采样频率分别为
和
。
3.要求用双线性变换法设计一个数字巴特沃思低通滤波器，其特性曲线如下图所示。
在通带内
，允许幅度误差小于1dB，在阻带
时衰减应大于15dB。
通带幅度归一化，使其在
处为1。
4.用模拟频带变换法，由二阶巴特沃思函数设计截止频率为
抽样频率为
的数字高通滤波器。
三.实验说明
1.求解实验内容1的MATLAB实现程序如下：
---------------------------------------------------------------------------------------------------------------
clear;
closeall
fp=10000;fs=40000;Rp=3;As=35;
[N,fc]=buttord（fp,fs,Rp,As,'s'）
[B,A]=butter（N,fc,'s'）;
[hf,f]=freqs（B,A,1024）;
subplot（3,2,1）;
plot（f,20*log10（abs（hf）/abs（hf
（1））））
grid;xlabel（'f/Hz'）;ylabel（'幅度（dB）'）;
axis（[0,50000,-40,5]）
line（[0,50000],[-3,-3]）;
---------------------------------------------------------------------------------------------------------------
程序运行结果：
N=3
fc=1.0441e+004
其频率特性曲线如下图所示:
频率单位：
2.求解实验内容2的MATLAB实现程序如下：
---------------------------------------------------------------------------------------------------------------
clear;
closeall
b=1000;a=[1,1000];
w=[0:
1000*2*pi];
[hf,w]=freqs（b,a,w）;
subplot（2,3,1）
plot（w/2/pi,abs（hf））;
grid;
xlabel（'f/（Hz）'）;ylabel（'幅度'）;
Fs0=[1000,500];
form=1:
2
Fs=Fs0（m）;
[d,c]=impinvar（b,a,Fs）;
wd=[0:
512]*pi/512;
hw1=freqz（d,c,wd）;
subplot（2,3,2）;
plot（wd/pi,abs（hw1）/abs（hw1
（1）））;holdon;
end
grid;xlabel（'f/（Hz）'）;
text（0.52,0.88,'T=0.002s'）;
text（0.12,0.54,'T=0.001s'）;
form=1:
2
Fs=Fs0（m）;
[f,e]=bilinear（b,a,Fs）;
wd=[0:
512]*pi/512;
hw2=freqz（f,e,wd）;
subplot（2,3,3）;
plot（wd/pi,abs（hw2）/abs（hw2
（1）））;holdon;
end
grid;xlabel（'f/（Hz）'）;
text（0.5,0.74,'T=0.002s'）;
text（0.12,0.34,'T=0.001s'）;
---------------------------------------------------------------------------------------------------------------
运行结果如下图所示:
f/Hz
（a）（b）（c）
模拟滤波器到数字滤波器的转换
由图（b）可见，对冲激响应不变法，采样频率越高（时间
越小），混叠越小；由图（c）可见，对双线性变换法，无频率混叠，但存在非线性失真。
3.求解实验内容3的MATLAB实现程序如下：
---------------------------------------------------------------------------------------------------------------
wp=0.2*pi;
ws=0.3*pi;
Rp=1;
As=15;
T=1;
Fs=1/T;
OmegaP=（2/T）*tan（wp/2）;
OmegaS=（2/T）*tan（ws/2）;
ep=sqrt（10^（Rp/10）-1）;
Ripple=sqrt（1/（1+ep*ep））;
Attn=1/（10^（As/20））;
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;
OmegaC=OmegaP/（（10^（Rp/10）-1）^（1/（2*N）））;
[B,A]=butter（N,OmegaC,'s'）;
W=（0:
500）*pi/500;
[H]=freqs（B,A,W）;
mag=abs（H）;
db=20*log10（（mag+eps）/max（mag））;
[b,a]=bilinear（B,A,T）;
[h,w]=freqz（b,a,1000,'whole'）;
h=（h（1:
501））';
w=（w（1:
501））';
m=abs（H）;
db=20*log10（（m+eps）/max（m））;
figure
（1）;
subplot（2,2,1）;plot（w/pi,mag）;title（'幅度'）
ylabel（’模拟滤波器’）；
axis（[0,1,0,1.1]）
set（gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]）;
set（gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]）;grid
subplot（2,2,2）;plot（w/pi,db1）;title（’幅度（dB）’）
axis（[0,1,-30,5]）
set（gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]）;
set（gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]）;grid
subplot（2,2,3）;plot（w/pi,m）;
xlabel（’频率单位：
pi’）;ylabel（’数字滤波器’）；
axis（[0,1,0,1.1]）
set（gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]）;
set（gca,’YTickmode’,’manual’,’YTick’,[0,Attn,Ripple,1]）;grid
subplot（2,2,4）;plot（w/pi,db2）;
xlabel（’频率单位:
pi’）;
axis（[0,1,-30,5]）
set（gca,’XTickMode’,’manual’,’XTick’,[0,0.2,0.3,1]）;
set（gca,’YTickmode’,’manual’,’YTick’,[-30,-15,-1,0]）;grid
---------------------------------------------------------------------------------------------------------------
运行结果如下图所示。
N=6
OmegaC=0.7273
巴特沃思模拟滤波器及利用双线性变换法设计的数字滤波器
4.求解实验内容4的MATLAB实现程序如下：
---------------------------------------------------------------------------------------------------------------
N=2;
Fs=500;
fch=200;
wch=2*pi*fch/Fs;
[z,p,k]=buttap（N）;
[b,a]=zp2tf（z,p,k）;
[h,w]=freqs（b,a,512）;
Omegach=2*
升级会员 