当前位置:首页 >> 信息与通信 >>

C语言编写FFT程序


DSP 课程作业
用 C 语言编写 FFT 程序
1,快速傅里叶变换 FFT 简介
快速傅氏变换(FFT) ,是离散傅氏变换的快速算法,它是根据离散傅氏变换的奇、 偶、虚、实等特性,对离散傅立叶变换的算法进行改进获得的。它对傅氏变换的理论并 没有新的发现,但是对于在计算机系统或者说数字系统中应用离散傅立叶变换,可以说 是进了一大步。 我们假设 x(n)为 N 项的复数序列,由 DFT 变换,任一 X(m)的计算都需要 N 次复 数乘法和 N-1 次复数加法,而一次复数乘法等于四次实数乘法和两次实数加法,一次复 数加法等于两次实数加法,即使把一次复数乘法和一次复数加法定义成一次“运算” (四 次实数乘法和四次实数加法) ,那么求出 N 项复数序列的 X(m),即 N 点 DFT 变换大约就 需要 N^2 次运算。当 N=1024 点甚至更多的时候,需要 N2=1048576 次运算,在 FFT 中, 利用 WN 的周期性和对称性,把一个 N 项序列(设 N=2k,k 为正整数) ,分为两个 N/2 项的 子序列,每个 N/2 点 DFT 变换需要(N/2)2 次运算,再用 N 次运算把两个 N/2 点的 DFT 变换组合成一个 N 点的 DFT 变换。 这样变换以后, 总的运算次数就变成 N+2 (N/2) 2=N+N2/2。 继续上面的例子,N=1024 时,总的运算次数就变成了 525312 次,节省了大约 50%的运算 量。而如果我们将这种“一分为二”的思想不断进行下去,直到分成两两一组的 DFT 运 算单元,那么 N 点的 DFT 变换就只需要 Nlog2N 次的运算,N 在 1024 点时,运算量仅有 10240 次,是先前的直接算法的 1%,点数越多,运算量的节约就越大,这就是 FFT 的优 越性。 2 , FFT 算法的基本原理 FFT 算法的基本思想:利用 DFT 系数的特性,合并 DFT 运算中的某些项,吧长序列的 DFT—>短序列的 DFT,从而减少其运算量。 FFT 算 法 分 类 : 时 间 抽 选 法 DIT: Decimation-In-Time ; 频 率 抽 选 法 DIF: Decimation-In-Frequency 按时间抽选的基-2FFT 算法 1、算法原理 设序列点数 N = 2L,L 为整数。 若不满足,则补零。N 为 2 的整数幂的 FFT 算法称基-2FFT 算法。将序列 x(n)按 n 的奇偶分成两组: x 2r ? x r N

? ? 1? ? x ? 2r ? 1? ? x2 ? r ?
N ?1 n ?0

r ? 0,1,...,
N ?1 n ?0

2

?1

则 x(n)的 DFT:

nk nk nk X ? k ? ? ? x ? n ?WN ? ? x ? n ?WN ? ? x ? n ?WN n ?0

N ?1

? ? x ? 2r ?W
N ?1 2

N ?1 2

2 ? ? x1 ? r ? ?WN ? ? WNk ? x2 ? r ? ?WN2 ? rk
2

r ?0

2 rk N

? ? ? x ? 2r ? 1?WN
r ?0 N
?1

N ?1 2

2 r ?1? k

rk

r ?0

r ?0

rk k rk ? ? x1 ? r ?WN / 2 ? WN ? x2 ? r ?WN / 2 r ?0 r ?0

N ?1 2

N ?1 2

? X1 ? k ? ? W X 2 ? k ?
k N

(r , k ? 0,1,...

N ? 1) 2

其中

X 1 (k ) ? ? x1 (r )W N ? ? x(2r )W Nrk
rk

N ?1 2

N ?1 2

X 1 (k ) ? ? x1 (r )W Nrk ? ? x(2r )W Nrk
r ?0
2

N ?1 2

r ?0

2

N ?1 2

r ?0

2

r ?0

2

(k ? 0,1,...

N ? 1) 2

再利用周期性求 X(k)的后半部分: N ? X 1 ? k ? , X 2 ? k ? 是以 为周期的 2 N? N? ? ? ? X1 ? k ? ? ? X1 ? k ? X 2 ? k ? ? ? X 2 ?k ? 2? 2? ? ?

又WN

k?

N 2

k k ? W WN ? ?WN

N 2 N

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

n为偶数

n为奇数

2) 、运算量 当 N = 2L 时,共有 L 级蝶形,每级 N / 2 个蝶形, 每 N N 复数乘法: 个蝶形有 1 次复数乘法 m2 次复数加法。 ? L ? log N
F

复数加法: 比较 DFT

aF ? NL ? N log 2 N

2

2

2

mF ( DFT ) N2 2N ? ? N mF ( FFT ) log 2 N log 2 N 2

3) 、算法特点 原位计算 蝶形运算两节点的第一个节点为 k 值,表示成 L 位二进制数,左移 L – m 位, 把右边空出的位置补零,结果为 r 的二进制数。
r ? X m (k ) ? X m?1 (k ) ? X m?1 ( j )WN ? r ? X m ( j ) ? X m?1 (k ) ? X m?1 ( j )WN

倒位序

x(n) n ? (n2 n1n0 ) 2

蝶形运算 对 N = 2L 点 FFT,输入倒位序,输出自然 序, 第 m 级运算每个蝶形的两节点距离为 2m–1 第 m 级运算:

r ? X m (k ) ? X m?1 (k ) ? X m?1 (k ? 2m?1 )WN ? m ?1 m ?1 r ? X m (k ? 2 ) ? X m?1 (k ) ? X m?1 (k ? 2 )WN

r WN 的确定

蝶形运算两节点的第一个节点为 k 值,表示成 L 位二进制数,左移 L – m 位,把右边空出的位置补零,结果为 r 的二进制数。 存储单元 输入序列 x(n) : N 个存储单元 r 系数 :N / 2 个存储单元 N

W

3,快速傅立叶变换的 C 语言实现方法 我们要衡量一个处理器有没有足够的能力来运行 FFT 算法,根据以上的简单介绍可以得 出以下两点: 1. 处理器要在一个指令周期能完成乘和累加的工作,因为复数运算要多次查表相乘 才能实现。 2. 间接寻址,可以实现增/减 1 个变址量,方便各种查表方法。FFT 要对原始序列进 行反序排列,处理器要有反序间接寻址的能力。 #include <stdio.h> #include <math.h> #include <stdlib.h>

#define N 1000 typedef struct { double real; double img; }complex;

void void void void void void void void void

fft(); ifft(); initW(); /*初始化变化核*/ change(); add(complex ,complex ,complex *); mul(complex ,complex ,complex *); sub(complex ,complex ,complex *); divi(complex ,complex ,complex *); output();

complex x[N], *W;/*输出序列的值*/ int size_x=0;/*输入序列的长度,只限2的N次方*/ double PI; int { int main() i,method;

system("cls"); PI=atan(1)*4;/*pi等于4乘以1.0的正切值*/ printf("Please input the size of x:\n"); /*输入序列的长度*/ scanf("%d",&size_x); printf("Please input the data in x[N]:(such as:5 6)\n"); /*输入序列对应的值*/ for(i=0;i<size_x;i++) scanf("%lf %lf",&x[i].real,&x[i].img); initW(); /*选择FFT或逆FFT运算*/ printf("Use FFT(0) or IFFT(1)?\n"); scanf("%d",&method); if(method==0) fft(); else ifft(); output(); return 0; }

/*进行基-2 FFT运算*/ void fft() { int i=0,j=0,k=0,l=0; complex up,down,product; change(); for(i=0;i< log(size_x)/log(2) ;i++) /*一级蝶形运算*/ { l=1<<i; for(j=0;j<size_x;j+= 2*l ) /*一组蝶形运算*/ { for(k=0;k<l;k++) /*一个蝶形运算*/ { mul(x[j+k+l],W[size_x*k/2/l],&product); add(x[j+k],product,&up); sub(x[j+k],product,&down); x[j+k]=up; x[j+k+l]=down; } } } }

void ifft() { int i=0,j=0,k=0,l=size_x; complex up,down; for(i=0;i< (int)( log(size_x)/log(2) );i++) /*一级蝶形运算*/ { l/=2; for(j=0;j<size_x;j+= 2*l ) /*一组蝶形运算*/ { for(k=0;k<l;k++) /*一个蝶形运算*/ { add(x[j+k],x[j+k+l],&up); up.real/=2;up.img/=2;

sub(x[j+k],x[j+k+l],&down); down.real/=2;down.img/=2; divi(down,W[size_x*k/2/l],&down); x[j+k]=up; x[j+k+l]=down; } } } change(); } /*初始化变化核*/ void initW() { int i; W=(complex *)malloc(sizeof(complex) for(i=0;i<size_x;i++) { W[i].real=cos(2*PI/size_x*i); W[i].img=-1*sin(2*PI/size_x*i); } } /*变址计算,将x(n)码位倒置*/ void change() { complex temp; unsigned short i=0,j=0,k=0; double t; for(i=0;i<size_x;i++) { k=i;j=0; t=(log(size_x)/log(2)); while( (t--)>0 ) { j=j<<1; j|=(k & 1); k=k>>1;

*

size_x);

} if(j>i) { temp=x[i]; x[i]=x[j]; x[j]=temp; } } }

void output() /*输出结果*/ { int i; printf("The result are as follows\n"); for(i=0;i<size_x;i++) { printf("%.4f",x[i].real); if(x[i].img>=0.0001) printf("+%.4fj\n",x[i].img); else if(fabs(x[i].img)<0.0001) printf("\n"); else printf("%.4fj\n",x[i].img); } } void add(complex a,complex b,complex *c) { c->real=a.real+b.real; c->img=a.img+b.img; } void mul(complex a,complex b,complex { c->real=a.real*b.real - a.img*b.img; c->img=a.real*b.img + a.img*b.real; } void sub(complex a,complex b,complex *c)

*c)

{ c->real=a.real-b.real; c->img=a.img-b.img; } void divi(complex a,complex b,complex *c) { c->real=( a.real*b.real+a.img*b.img )/( b.real*b.real+b.img*b.img); c->img=( a.img*b.real-a.real*b.img)/(b.real*b.real+b.img*b.img); }


赞助商链接
相关文章:
C语言实现FFT变换
C语言实现FFT变换_计算机软件及应用_IT/计算机_专业资料。#include <iom128.h>...此程序包采用联合体的形式表示一个复数,输入为自然顺序的复 数(输入实数是可令...
C语言实现FFT(快速傅里叶变换)
C语言实现FFT(快速傅里叶变换)_信息与通信_工程科技_专业资料。C语言实现FFT(...此程序包采用联合体的形式表示一个复数,输入为自然顺序的复 数(输入实数是可令...
FFT的C语言算法实现
FFT的C语言算法实现_信息与通信_工程科技_专业资料。一个用C语言写的简易的FFT的算法,FFT 的 C 语言算法实现程序如下: /***FFT***/ #include <stdio.h> ...
用C语言编写FFT算法
以下是 512 点的单边谱 FFT 算法,在 TMS320C6713 的 DSP 开发板上调试通过。 顺便问一句,有人会编 welch 功率谱密度的 C 语言实现程序么? #include "math...
C语言实现FFT
C语言实现FFT_信息与通信_工程科技_专业资料。#include <iom128.h> #include ...FFT_N 128 //定义福利叶变换的点数 struct compx {float real,imag;}; //...
实序列FFT算法的C语言实现
实序列FFT算法的C语言实现_数学_自然科学_专业资料 暂无评价|0人阅读|0次下载|举报文档 实序列FFT算法的C语言实现_数学_自然科学_专业资料。...
C语言FFT程序
C语言FFT程序_工学_高等教育_教育专区。用C语言变得作快速傅立叶变换的程序 /*** 正弦函数快速傅立叶变换 C 程序 函数简介: 输入正弦函数的频率 FC,然后由采...
CCS上FFT的C语言实现
CCS上FFTC语言实现 - 从两个文件分别读入复数的实部和虚部,然后进行FFT运算。... CCS上FFTC语言实现_电脑基础知识_IT/计算机_专业资料。从两个文件分别读入复...
FFT的C语言编程
FFTC 语言编程 1.程序: #include <stdio.h> #include <math.h> #define re 0 /* re=0,用 re 表示实部*/ #define im 1 /* im=1,用 im 表...
C语言实现FFT
快速傅里叶变换FFTC语言... 14页 1财富值 FFT原理及实现 17页 免费如要投诉违规内容,请到百度文库投诉中心;如要提出功能问题或意见建议,请点击此处进行反馈。...
更多相关标签: