信号分析与处理实验报告程序+报告内容(非常全)

非常全

实验四 离散傅里叶变换及其快速算法 一、实验目的1、掌握计算序列的离散傅里叶变换(DFT)的方法。 2、掌握实现时间抽取快速傅里叶变换(FFT)的编程方法。 3、通过编写程序，深入理解快速傅里叶变换算法(FFT)的含义，比较 DFT 与 FFT 的运算效率。 4、加深对 DFT 与序列的傅里叶变换和 Z 变换之间关系的理解。

二、实验原理 2.1 离散傅里叶变换(DFT)算法有限列长为 N 的序列 x(n) 的 DFT 变换为X (k )

x(n)Wn 0

N 1

nk n

k 0,1,2 N 1

其逆变换为x ( n) 1 N

X (k )Wk 0

N 1

nk N

n 0,1, N 1

2.2 快速傅里叶变换基本思路当 N 很大的时候， 求一个 N 点的 DFT 要完成 N N 次复数乘法和 N ( N 1) 次 复数加法，其计算量相当大。快速傅里叶变换巧妙地利用 WN 因子的周期性和对 称性，构造了一个 DFT 快速算法。其基本思路如下： 1、利用 WN nk 的对称性使 DFT 运算中有些项合并

WN2、利用WNnk

k ( N n )

WN

kn

(WN ) kn

的周期性和对称性使长序列的 DFT 分解为更小点数的 DFT

WN

nk

WN

k ( n N )

(WN

( k N )n

)

为了讨论方便，设 N 2 ,其中 为整数。如果不满足这个条件，可以认为得 加上若干零点来达到。由 DFT 的定义知：

非常全

X (k )

x(n)Wn 0

N 1

nk n

k 0,1,2 N 1

其中 x(n) 是列长为 N (n 0,1 N 1) 的输入序列，把它按 n 的奇偶分成两个序 列：

x(2r ) x1 (r ) x(2r 1) x2 (r )2 又由于 WN e j 2 N j

r 0,1,

N 1 2

e

2 N 2

WN2

,

则 X (k )

n 0 n为偶数

N 1

x(n)Wn

nk

n 0 n为奇数

N 1

nk k x(n)WN X 1 (k ) WN X 2 (k )

上式表明了一个 N 点的 DFT 可以被分解为两个 N/2 点的 DFT。同时，这两个 N/2 点的 DFT 按照上式又可以合成为一个 N 点的 DFT。 为了要用点数为 N/2 点的 X 1 (k)、X 2 (k)来表达 N 点的 X(k)值还必须要用 W 系rk WN 数的周期性，即 W N 2 2 r (k N ) 2

这样可得N N N 2 2 r (k ) N rk X1 ( k ) x1 (r )WN 2 x1 (r )WN 2 r 0 r 0 2 2

即 同理可得

X1 ( X2(k

N k ) X 1 (k ) 2

N k ) X 2 (k ) 2

另外再加上 W N 的对称性WN 2( N k ) N k k WN2 WN WN

就可以将 X(k)的表达式分为前后两个部分： 前半部分 后半部分k X (k ) X 1 (k ) WN X 2 (k )

r 0,1,

N 1 2

N ( k ) N N N X ( k ) X 1 ( k ) WN 2 X 2 ( k ) 2 2 2

非常全

X 1 (k ) W X 2 (k )k N

r 0,1,

N 1 2

N 由以上分析可见，只要求出区间 0, 1 内各个整数 k 值所对应的 X 1 (k)、 2

X 2 (k)的值，即可求出 0, N 1 区间内的全部 X(k)值，这一点恰恰是 FFT 能

大量节省计算的关键所在。 2.3 时间抽取的快速傅里叶变换计算机算法 1、程序输入序列的元素数目必须为 2 的整数次幂，即 N 2M ,整个运算需要 M 级蝶形运算； 2、输入序列应该按二进制的码位倒置排列，输出序列按自然序列排列； 每个蝶形运算的输出数据占用其他输入数据的存储单元，实现“即位运算” ; 每级包括 N/2 个基本蝶形运算，共有 M*N/2 个基本蝶形运算； 3、第 L 级中有 N / (2 L) 个群，群与群的间隔为 2 L ; 4、同一级的各个群的系数 Wn 分布相同，第 L 级的群中有 2L 1 个系数； 5、处于第 L 级的群的系数是 WN N ( p 1)/2L

(p=1,2,3, ., 2L 1 );

6、对于第 L 级的蝶形运算，两个输入数据的间隔为 2L 1 。

三、实验程序及注释 3.1 用 C 语言实现 FFT 与 IFFT#include "stdio.h" #include "math.h" #include <stdlib.h> /*清屏命令 clrscr()需加载的预处理函数*/ #define swap(a,b) temp=(a);(a)=(b);(b)=temp /*宏定义的交换函数*/ /*快速傅利叶变换程序，数组 A、B 分别是带变换序列的实部和虚部 ap 是变换 的类别，ap=1 是做 FFT; ap=-1 是做 IFFT */ void fft(float A[],float B[],unsigned M,int ap) {unsigned long N,I,J,K,L,LE,LE1,P,Q,R; int i; float Wr,Wi,W1r,W1i,WTr,WTi,theta,Tr,Ti,temp; N=pow(2,M); /*N=2M 是序列的总长度*/ J=0; for(I=0;I<N-1;I++) /*码位倒置*/ {if(J>I) {swap(A[I],A[J]); swap(B[I],B[J]);}

非常全

K=N>>1; /*即 k=k/2 */ while(K>=2&&J>=K) {J-=K; K>>=1;} J+=K;} for(L=1;L<=M;L++) /*外层循环由级数 L 控制，执行 M 次*/ {LE=1<<L; /* LE=2L 是群间隔*/ LE1=LE/2; /* LE1=2L-1 是每个群的系数 W 数目*/ Wr=1.0; Wi=0.0; theta=(-1)*ap*3.1415926536/LE1; /*ap 用于在 IFFT 时改变 极性*/ W1r=cos(theta); W1i=sin(theta); for(R=0;R<LE1;R++) /*中层循环由群系数 控制，执行 次*/ {for(P=R;P<N-1;P+=LE) /*R 是群系数的编号，P、Q 是基本蝶形运 算两个输入数据在数组中的编号，循环每 次完成同一个系数 的蝶形运算*/ {Q=P+LE1; Tr=Wr*A[Q]-Wi*B[Q]; Ti=Wr*B[Q]+Wi*A[Q]; /*Tr、Ti 是 的实部和虚部*/ A[Q]=A[P]-Tr; /*即 */ B[Q]=B[P]-Ti; A[P]+=Tr; /*即 */ B[P]+=Ti;} WTr=Wr; /*Wr、Wi 是 的实部和虚部*/ WTi=Wi; /*下面用和差化积公式求 的实虚部 Wr、Wi*/ Wr=WTr*W1r-WTi*W1i; /* */ Wi=WTr*W1i+WTi*W1i; /* */ } } if(ap==-1){ /*ap=-1 是做 IFFT,应使变换后的序列除以 N*/ for(i=0;i<N;i++){ A[i]/=N; B[i]/=N; } } return; } main(){ int i,t,ap; float A[100],B[100]; unsigned M; system("cls"); /*清屏*/ printf("\nplease input bian huan lei bie(FFT:1 or IFFT:-1):");

非常全

scanf("%d",&ap); /*输入序列变换的类别*/ print …… 此处隐藏：13628字，全部文档内容请下载后查看。喜欢就下载吧 ……