实验四离散傅里叶变换.docx
《实验四离散傅里叶变换.docx》由会员分享，可在线阅读，更多相关《实验四离散傅里叶变换.docx（21页珍藏版）》请在冰豆网上搜索。
实验四离散傅里叶变换
实验四离散傅里叶变换
一、实验目的
（1）掌握计算序列的离散傅里叶（DFT）的方法；
（2）掌握实现时间抽取快速傅里叶变换（FFT）的编程方法，比较DFT和FFT的计算效率；
（3）加深对DFT与序列的傅里叶变换和Z变换之间的关系的理解；
（4）复习复数序列的运算方法。
二、实验原理
一个长度为N的序列x（n）的离散傅里叶变换（DFT）及其反变换（IDFT）为
其中，
。
计算DFT的快速算法称为快速傅里叶变换（FFT），主要有时间抽取算法和频率抽取算法。
在编制程序时要注意体会码位倒置、即位运算和蝶形图的实现方法。
单位圆上得Z变换即为序列的傅里叶变换，序列的离散傅里叶变换（DFT）则是序列的傅里叶变换的离散化
通过本实验可以验证这一关系。
三、实验内容与步骤
（1）编制基2时间抽取算法的FFT子程序fft（x,M）,参数x为被分析序列x（n），M为蝶形运算的级数，它与序列长度N的关系是N=
。
调用fft（x，M）的结果是生成一个长度相等的变换序列X（k）。
编制验证子程序功能的主程序。
（2）在子程序fft（x，M）的基础上编制计算IFFT的子程序ifft（X，M），参数X为被分析序列X（k），M为蝶形运算的级数。
编程验证子程序功能的主程序，利用所编制的子程序ifft（X，N）计算步骤
（1）中序列X（k）的离散傅里叶反变换x（n），将计算结果与步骤
（1）对比。
（3）画出步骤
（1）中序列X（k）的幅频特性和相频特性。
四、实验结果
（1）给出本实验内容要求的2个子程序并加以注释
程序一：
fft
#include
#include
#defineswap（a,b）{floatT;T=（a）;（a）=（b）;（b）=T;}
voidfft（floatA[],floatB[],unsignedM）//蝶形运算程序，A存实部，B存虚部，M是级数
{
unsignedlongN,I,J,K,L,LE,LE1,P,Q,R;//LE为群间隔，LE1为群内碟间隔
floatWr,Wi,Wlr,Wli,WTr,WTi,theta,Tr,Ti;
N=1<J=0;
for（I=0;I{
if（J>I）//为了避免交换重复，只有标号J大于I时才执行交换操作
{
swap（A[I],A[J]）;//实部交换
}
if（J>I）
{
swap（B[I],B[J]）;//虚部交换
}
K=N>>1;//K=N>>1表示K=N/2
while（K>=2&&J>=K）//while循环表示须向次高位进一位
{
J-=K;K>>=1;//K>>=1表示K=K/2
}
J+=K;
}
for（L=1;L<=M;L++）//for循环为M级FFT运算
{
LE=1<LE1=LE/2;
Wr=1.0;
Wi=0.0;
theta=（-1）*3.1415926536/LE1;
Wlr=cos（theta）;
Wli=sin（theta）;//将
换为cos（
）-jsin（
）的形式
for（R=0;R{
for（P=R;P{
Q=P+LE1;
Tr=Wr*A[Q]-Wi*B[Q];
Ti=Wr*B[Q]+Wi*A[Q];
A[Q]=A[P]-Tr;
B[Q]=B[P]-Ti;
A[P]+=Tr;
B[P]+=Ti;
}
WTr=Wr;
WTi=Wi;
Wr=WTr*Wlr-WTi*Wli;
Wi=WTr*Wli+WTi*Wlr;
}
}
return;
}
voidmain（）//主程序
{
intM,I,N;
floatA[100],B[100];
printf（"请输入M：
"）;//输入M的值
scanf（"%d",&M）;
N=pow（2,M）;//
printf（"\n请输入序列实部（N个）：
"）;//输入序列实部
for（I=0;Iscanf（"%f",&A[I]）;
for（I=0;Iprintf（"\n请输入序列虚部（N个）：
"）;//输入序列虚部
scanf（"%f",&B[I]）;
printf（"则您所输入的序列为：
\n"）;
for（I=0;Iprintf（"x[%d]=%.2f+j%.2f\n",I,A[I],B[I]）;
fft（A,B,M）;
printf（"则经FFT计算后，结果为：
\n"）;//输出序列
for（I=0;Iprintf（"X[%d]=%.2f+j%.2f\n",I,A[I],B[I]）;
}
运行结果如下：
程序二：
IFFT
#include
#include
voidfft（floatA[],floatB[],unsignedM）;
voidmain（）
{
floatA[1024],B[1024];
unsignedlongN,i,M;
printf（"输入M="）;
scanf（"%d",&M）;
N=1<printf（"\n请输入序列实部（”N”个）;"）;
for（i=0;iscanf（"%f",&A[i]）;
printf（"\n请输入序列虚部（”N”个）;"）;
for（i=0;iscanf（"%f",&B[i]）;
for（i=0;iB[i]=0-B[i];
printf（"则您所输入的序列为：
\n"）;
for（i=0;iprintf（"x[%d]=%.2f+j%.2f\n",i,A[i],B[i]）;
fft（A,B,M）;
printf（"则经IFFT计算后，结果为：
\n"）;//输出序列
for（i=0;iprintf（"X[%d]=%.2f+j%.2f\n",i,A[i]/N,B[i]/N）;
}
#defineswap（a,b）T=（a）;（a）=（b）;（b）=T
voidfft（floatA[],floatB[],unsignedM）//套用fft函数求ifft
{
unsignedlongN,I,J,K,L,LE,LE1,P,Q,R;
floatWr,Wi,W1r,W1i,WTr,WTi,theta,Tr,Ti;
N=1<J=0;
for（I=0;I{
if（J>I）
{
floatT;
swap（A[I],A[J]）;
swap（B[I],B[J]）;
}
K=N>>1;
while（K>=2&&J>=K）
{
J-=K;K>>=1;
}
J+=K;
}
for（L=1;L<=M;L++）
{
LE=1<LE1=LE/2;
Wr=1.0;
Wi=0.0;
theta=（-1）*3.1415926536/LE1;
W1r=cos（theta）;
W1i=sin（theta）;
for（R=0;R{
for（P=R;P{
Q=P+LE1;
Tr=Wr*A[Q]-Wi*B[Q];
Ti=Wr*B[Q]+Wi*A[Q];
A[Q]=A[P]-Tr;
B[Q]=B[P]-Ti;
A[P]+=Tr;
B[P]+=Ti;
}
WTr=Wr;
WTi=Wi;
Wr=WTr*W1r-WTi*W1i;
Wi=WTr*W1i+WTi*W1r;
}
}
return;
}
运行结果为：
实验五、IIR数字滤波器的设计
一、实验目的
（1）掌握设计IIR数字滤波器的冲激响应不变法；
（2）掌握设计IIR数字滤波器的双线变换法；
（3）掌握IIR数字滤波器的现实方法。
二、实验原理
1、设计IIR数字滤波器的冲激响应不变法的步骤
（1）明确设计指标，由于需要设计符合指标的模拟滤波器H（s），因此，首先要将在数字域给出的设计指标转换为模拟域设计的设计指标；
（2）设计符合指标要求的模拟滤波器
，即确定滤波器的阶数N和截止角频率
。
（3）将模拟滤波器
变换为数字滤波器H（z），首先计算
在S平面左半平面N个极点
及留数，然后代入下式得到H（z）
式中
为采样时间；
（4）验证数字滤波器的频率特性，首先计算验证数字滤波器的幅频特性
然后选取包括通带截止角频率和阻带截止角频率在内的若干频率点，计算数字滤波器的幅频特性，并与设计指标比较。
如果不满足要求，就要改进设计，重复上述步骤，直到满足要求。
2、设计IIR数字滤波器的双线性变换法的步骤
（1）明确设计指标，双线性变换法需要对设计指标进行预畸变，才能将数字域给出的设计指标变换为模拟域的设计指标，预畸变公式为
其中，
为采样时间，
和
分别为数字滤波器的通带截止角频率和阻带截止角频率，
和
分别为模拟滤波器的通带截止角频率和阻带截止角频率。
（2）设计符合实际指标要求的模拟滤波器
，即确定滤波器的阶数N和截止角频率
。
（3）将模拟滤波器
变换为数字滤波器H（z），首先计算
在S平面左半平面N个极点
及留数，然后代入下式得到H（z）
（4）验证数字滤波器的频率特性
在设计出数字滤波器的系统函数H（z）后，为了实现数字滤波器算法，需要选择一种实现结够形式。
常用的结构形式有直接型、级联型、并联型等。
三、实验内容与步骤
（1）利用冲激响应不变法设计一个巴特沃斯型IIR数字滤波器，设计指标要求在通带
<=0.2
内，允许幅度误差小于1dB，在阻带
>=0.3
内衰减应大于15dB，通带幅度归一化，使其在w=0处为0dB。
画出数字滤波器的频率特性，选择并联结构实现设计的数字滤波器。
（2）利用双线性变换法设计一个巴特沃斯型IIR数字滤波器，设计指标同步骤
（1）。
画出数字滤波器的频率特性，选择级联型结构实现设计的数字滤波器。
（3）为了验证所设计数字滤波器的实现功能，首先利用实验一的方法将下列分别位于通带、过渡带、阻带的三个模拟信号数字化
将得到的三个序列作为激励信号分别输入到所设计的数字滤波器的输入端，计算输出地响应信号，观察序列幅度的变化。
四、实验结果
（1）给出程序并加以注释
用matlab实现
采样频率10Hz，通带截止频率fp=0.2*piHz，阻带截止频率fs=0.3*piHz
通带衰减小于1dB，阻带衰减大于15dB
%冲激响应不变法
wp=0.2*pi;ws=0.3*pi;
Rp=1;%通带最大衰减Rp=1dB
Rs=15;%阻带最小衰减Rs=20dB
TS=1;FS=1/TS;
Wp=wp/TS;Ws=ws/TS;
[n,wn]=buttord（Wp,Ws,1,15,'s'）%计算巴特沃夫型滤波器的阶数n和截止频率wn
[Bs,As]=butter（n,wn,'s'）%求系统函数，Bs为分子向量，As为分母向量
W=0:
0.01:
pi;
[Bz,Az]=impinvar（Bs,As,FS）%冲激响应不变法
hz=filt（Bz,Az）%求系统函数H（z）
s2=freqz（Bz,Az,W）;
subplot（211）;
plot（W,20*log10（abs（s2）））;%绘制数字滤波器的幅度响应
grid;
title（'冲激响应不变法数字滤波器幅频特性','fontsize',12）;
xlabel（'w（rad）'）;
ylabel（'幅度（dB）'）;
subplot（212）;
plot（W,angle（s2））;%数字低通滤波器的相位响应
grid;
title（'冲激响应不变法数字滤波器相频特性','fontsize',12）;
xlabel（'w（rad）'）;
ylabel（'相位（rad）'）;
hz=sym（filt（Bz,Az））;
h=iztrans（hz）;%求hz的Z反变换
h=simplify（h）%化简后的h
部分运行结果;
4.108e-015+0.0006584z^-1+0.0105z^-2+0.01672z^-3
+0.004232z^-4+0.0001062z^-5
-----------------------------------------------------------
1-3.344z^-1+5.018z^-2-4.219z^-3+2.073z^-4
0.56z^-5+0.0647z^-6
%双线性变换实现ButterWorth低通
%使用双线性变换法由模拟滤波器原型设计数字滤波器
T=0.1;FS=1/T;pi=3.14;
fp=0.2*pi;fs=0.3*pi;
wp=fp/FS*2*pi;
ws=fs/FS*2*pi;
Rp=1;%通带衰减
As=15;%阻带衰减
OmegaP=（2/T）*tan（wp/2）;%PrewarpPrototypePassbandfreq%频率预畸变
OmegaS=（2/T）*tan（ws/2）;%PrewarpPrototypeStopbandfreq
%设计butterworth低通滤波器原型
N=ceil（（log10（（10^（Rp/10）-1）/（10^（As/10）-1）））/（2*log10（OmegaP/OmegaS）））;%求阶数N
OmegaC=OmegaP/（（10^（Rp/10）-1）^（1/（2*N）））;%求
[z,p,k]=buttap（N）;%获取零极点参数
p=p*OmegaC;
k=k*OmegaC^N;
B=real（poly（z））;
b0=k;
cs=k*B;
ds=real（poly（p））;
%双线性变换
[b,a]=bilinear（cs,ds,FS）;%映射为数字的
figure
（1）;%绘制频率特性
freqz（b,a,512,FS）;
%进行滤波验证
figure
（2）;
f1=50;
f2=100;
f3=200;
pi=3.14;
n=0:
（80/1.25）;
x1=4*sin（2*pi*f1*n）;
x2=4*sin（2*pi*f2*n）;
x3=4*sin（2*pi*f3*n）;
subplot（3,2,1）;
stem（x1,'.'）
title（'x1输入信号'）
y1=filter（b,a,x1）;
subplot（3,2,2）;
stem（y1,'.'）
title（'x1滤波之后的信号'）;
subplot（3,2,3）;
stem（x2,'.'）
title（'x2输入信号'）
y2=filter（b,a,x2）;
subplot（3,2,4）;
stem（y2,'.'）
title（'x2滤波之后的信号'）;
subplot（3,2,5）;
stem（x3,'.'）
title（'x3输入信号'）
y3=filter（b,a,x3）;
subplot（3,2,6）;
stem（y3,'.'）
title（'x3滤波之后的信号'）
运行结果：
>>
图像分析：
x1为低频信号，位于通带，几乎没有衰减；x2为过渡带的信号也几乎没有衰减，x3为高频信号，位于阻带，信号大幅衰减。
符合滤波器要求。
选择级联型结构实现数字滤波器
#include
#include
voidIIRDF（floatA[],unsignedlongN）;
voidmain（）
{
floatA[1024];
intN,i,m=100;//改变m值来改变输入信号的频率，取100,200,400//
printf（"函数是：
"）;//提示输入函数//
printf（"x（t）=4*sin（%d*3.1415926*1.25*t）\n",m）;//输入信号
for（i=0;i*1.25<=80;i++）
{
A[i]=4*sin（m*3.1415926*1.25*i/1000）;//对信号采样离散化
}
N=80/1.25;
for（i=0;i<=N;i++）
{
printf（"x（%d）=%.4f",i,A[i]）;//输出离散化的结果
if（（i+1）%3==0）
printf（"\n"）;
}
printf（"\n"）;
IIRDF（A,N）;
printf（"theresultsofIIRDFis\n"）;//提示输出结果//
printf（"\n"）;
for（i=0;i<=N;i++）
{
printf（"y（%d）=%.4f",i,A[i]）;/*输出滤波后的结果*/
if（（i+1）%3==0）
printf（"\n"）;
}
printf（"\n"）;
}
voidIIRDF（floatA[],unsignedlongN）//级联型低通滤波器
{
unsignedlongn;
floatx[3]={0,0,0},y1[3]={0,0,0},y2[3]={0,0,0},y[3]={0,0,0};
for（n=0;n{
x[0]=A[n];
y1[0]=1.31432*y1[1]-0.71489*y1[2]+0.08338*x[0]+0.16676*x[1]+0.08338*x[2];
//第一级
x[2]=x[1];
x[1]=x[0];
y2[0]=1.0541*y2[1]-0.37534*y2[2]+0.08338*y1[0]+0.16676*y1[1]+0.08338*y1[2];//第二级
y1[2]=y1[1];
y1[1]=y1[0];
y[0]=0.94592*y[1]-0.23422*y[2]+0.08338*y2[0]+0.16676*y2[1]+0.08338*y2[2];
//第三级
y2[2]=y2[1];
y2[1]=y2[0];
y[2]=y[1];
y[1]=y[0];
A[n]=y[0];
}
}
结果分析
当输入为100时，结果为：
用matlabfilterdesignedtool仿真如下：
（2）比较应用冲激响应不变法和双线性变换法所设计的IIR数字滤波的幅频特性和相频特性
答：
由理论分析及以上仿真效果知：
a在冲激响应不变法中，数字角频率和模拟角频率之间呈线性关系。
这样，所设计的数字滤波器能很好的保持模拟滤波器原型的频率特性。
但是，由于S平面与Z平面之间的变换不是一一对应的，从模拟滤波器到数字滤波器的变换会发生频率混叠。
b在双线性变换法中，S平面与Z平面之间存在一一映射关系，从而克服了冲激响应不变法产生的频谱混叠现象；但是，由于数字角频率与模拟角频率之间不是线性关系，因此，在变换过程中会发生频率失真。
通过采用预畸变的补偿方法，可以对具有片状常数频率特性的滤波器进行校正。
（3）讨论模拟角频率与数字角频率的关系
答：
概念：
模拟角频率Ω：
每秒经历多少弧度，单位rad/s；
数字频率w：
每个采样点间隔之间的弧度，单位rad。
（w称为正弦序列的数字域频率（也称数字频率），单位是弧度（rad），它表示序列变化的速率，或者说表示相邻两个序列值之间变化的弧度数。
）
关系：
数字频率和模拟角频率的关系：
w=Ω*T。
（这个式子具有普遍意义，它表示凡是由模拟信号采样得到的序列，模拟角频率Ω与序列的数字域频率ω成线性关系。
）
升级会员 