admin 管理员组

文章数量: 887021


2023年12月16日发(作者:paycheck)

我自己的一些详细标注,有利于深入了解FFT,后面附加几位网友对FFT的理解及源代码,让广大朋友更迅速的掌握FFT

#include "DSP281x_Device.h" // DSP281x Headerfile Include File,添加所有头文件

#include "DSP281x_Examples.h" // DSP281x Examples Include File,条件编译而已

#include "f2812a.h" //一些变量的宏定义而已

#include"math.h"

#define PI 3.1415926 //前变后常

#define SAMPLENUMBER 128

//#define SAMPLENUMBER 512

void InitForFFT();

void MakeWave();

//void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER]);

int INPUT[SAMPLENUMBER],DATA[SAMPLENUMBER];

float fWaveR[SAMPLENUMBER],fWaveI[SAMPLENUMBER],w[SAMPLENUMBER];

float sin_tab[SAMPLENUMBER],cos_tab[SAMPLENUMBER];

//逐级计算FFT,一级一级递推

void FFT(float dataR[SAMPLENUMBER],float dataI[SAMPLENUMBER])

{

int x0,x1,x2,x3,x4,x5,x6,xx;

int i,j,k,b,p,L;

float TR,TI,temp;

/********** following code invert sequence ************///倒序

for ( i=0;i

{ //128七位二进制表示,/号代表右移嘛

x0=x1=x2=x3=x4=x5=x6=0;

x0=i&0x01; x1=(i/2)&0x01; x2=(i/4)&0x01; x3=(i/8)&0x01;x4=(i/16)&0x01;

x5=(i/32)&0x01; x6=(i/64)&0x01;

xx=x0*64+x1*32+x2*16+x3*8+x4*4+x5*2+x6;//最低位,最高位反过来

dataI[xx]=dataR[i];

}

for ( i=0;i

{

dataR[i]=dataI[i]; dataI[i]=0; //对应过来

}

/************** following code FFT *******************/

for ( L=1;L<=7;L++ )

{ /* for(1) */

b=1; i=L-1;/* b的意义非常重大,b表示当前层不同旋转因子的个数 */

while ( i>0 )

{

b=b*2; i--;

} /* b= 2^(L-1) */

for ( j=0;j<=b-1;j++ ) /* for (2) */

{

p=1; i=7-L;

while ( i>0 ) /* p=pow(2,7-L)*j; */

{

p=p*2; i--;

}

p=p*j;

for ( k=j;k<128;k=k+2*b ) /* for (3) */

{

TR=dataR[k]; TI=dataI[k]; temp=dataR[k+b];

dataR[k]=dataR[k]+dataR[k+b]*cos_tab[p]+dataI[k+b]*sin_tab[p];

dataI[k]=dataI[k]-dataR[k+b]*sin_tab[p]+dataI[k+b]*cos_tab[p];

dataR[k+b]=TR-dataR[k+b]*cos_tab[p]-dataI[k+b]*sin_tab[p];

dataI[k+b]=TI+temp*sin_tab[p]-dataI[k+b]*cos_tab[p]; //递推嘛,防止立马调用结果,

} /* END for (3) */ //引入一个中间变量存原始值,

} /* END for (2) */ //防止上一步对下一步的影响

} /* END for (1) */

for ( i=0;i

{

w[i]=sqrt(dataR[i]*dataR[i]+dataI[i]*dataI[i]);

}

} /* END FFT */

main()

{

int i;

InitForFFT();

MakeWave();

for ( i=0;i

{

fWaveR[i]=INPUT[i];

fWaveI[i]=0.0f;

w[i]=0.0f;

}

}

FFT(fWaveR,fWaveI);//输入波形进行FFT变换,此处引入起始实参即可递推下去

for ( i=0;i

{

DATA[i]=w[i];//变换后的波形转换到输出接口

}

while ( 1 ); // break point

//旋转因子事先初始化好,方便调用

void InitForFFT()

{

int i;

for ( i=0;i

{

sin_tab[i]=sin(PI*2*i/SAMPLENUMBER);

}

}

cos_tab[i]=cos(PI*2*i/SAMPLENUMBER);//旋转因子事先初始化好,方便调用

//利用这个,确实能产生各种各样的谐波

void MakeWave() //利用这个,确实能产生各种各样的谐波

{

int i;

for ( i=0;i

{

INPUT[i]=sin(PI*2*i/SAMPLENUMBER)*1024

+sin(PI*2*i/SAMPLENUMBER*3)*1024/3

+sin(PI*2*i/SAMPLENUMBER*5)*1024/5

+sin(PI*2*i/SAMPLENUMBER*7)*1024/7

+sin(PI*2*i/SAMPLENUMBER*9)*1024/9;

//INPUT[i]=sin(PI*2*i/SAMPLENUMBER*3)*1024;

}

}

FFT原理及实现(Radix-2)

哈! 经过连续几个晚上的奋战, 终于弄懂了FFT推导过程及实现! Happy☺

基2 FFT总的思想是将输入信号对半分割, 再对半分割, 再再对半分割(以下省略10000个再再...☺) 直至分割到2点.

两点DFT简化

假设输入为x[0],x[1]; 输出为X[0],X[1]. 伪代码如下 :

// ------------------------------------------------------------------

#define N 2

#define PI 3.1415926

// ------------------------------------------------------------------

int i, j

for(i=0, X[i]=0.0; i

for(j=0; j

X[i] += x[j] * ( cos(2*PI*i*j/N) - sin(2*PI*i*j/N) );

注意到(我想Audio编解码很多时候都是对cos,sin进行优化!)

i=0

j=0

cos 1

sin 0

tw 1

i=1

cos 1

Sin 0

tw 1

j=1

cos 1

sin 0

tw 1

cos -1

sin 0

tw -1

X[0] = x[0]*(1-0) + x[1]*(1-0) = x[0] + 1*x[1];

X[1] = x[0]*(1-0) + x[1]*(-1-0) = x[0] - 1*x[1];

这就是单个2点蝶形算法.

FFT实现流程图分析(N=8, 以8点信号为例)

FFT implementation of an 8-point DFT as two 4-point DFTs and four 2-point

DFTs

8点FFT流程图(Layer表示层, gr表示当前层的颗粒)

下面以LayerI为例.

LayerI部分, 具有4个颗粒, 每个颗粒2个输入

(注意2个输入的来源, 由时域信号友情提供, 感谢感谢☺)

我们将输入x[k]分为两部分x_r[k], x_i[k]. 具有实部和虚部, 时域信号本没有虚部的, 因此可以让x_i[k]为0. 那么为什么还要画蛇添足分为实部和虚部呢? 这是因为

LayerII, LayerIII的输入是复数, 为了编码统一而强行分的.当然你编码时可以判断当前层是否为1来决定是否分. 但是我想每个人最后都会倾向分的.

旋转因子 tw = cos(2*PI*k/N)-j*sin(2*PI*k/N); 也可以分为实部和虚部, 令其为tw_r, tw_i;

则tw = tw_r - j*tw_i;

X[k] = (x_r[k] + j*x_i[k]) + (tw_r–j*tw_i) * (x_r[k+N/2]+j*x_i[k+N/2])

X_R[k] = x_r[k] + tw_r*x_r[k+N/2] + tw_i*x_i[k+N/2];

X_I[k] = x_i[k] - tw_i*x_r[k+N/2] + tw_r*x_i[k+N/2];

LayerII部分, 具有2个颗粒, 每个颗粒4个输入

(注意4个输入的来源, 由LayerI友情提供, 感谢感谢☺)

LayerIII部分, 具有1个颗粒, 每个颗粒8个输入

(注意8个输入的来源, 由LayerII友情提供, 感谢感谢☺)

LayerI, LayerII, LayerIII从左往右, 蝶形信号运算流非常明显!

假令输入为x[k], x[k+N/2], 输出为X[k], X[k+N/2]. x[k]分解为x_r[k], x_i[k]部分

则该蝶形运算为

X[k]

= (x_r[k]-j*x_i[k]) +

(x_r[k+N/2]-j*x_i[k+N/2])*(cos(2*PI*k/N)-j*sin(2*PI*k/N));

再令cos(2*PI*k/N)为tw1, sin(2*PI*k/N)为tw2 则

X[k] = (x_r[k]-j*x_i[k]) + (x_r[k+N/2]-j*x_i[k+N/2])*(tw1-j*tw2);

X_R[k] = x_r[k] + x_r[k+N/2]*tw1 - x_i[k+N/2]*tw2;

X_I[K] = x_i[k]

x_r[k] = x_r[k] + x_r[k+b]*tw1 + x_i[k+b]*tw2;

x_i[k] = x_i[k] - x_r[k+b]*tw2 + x_i[k+b]*tw1;

譬如8点输入x[8]

1. 先分割成2部分: x[0], x[2], x[4], x[6] 和 x[1], x[3], x[5], x[7]

2. 信号x[0], x[2], x[4], x[6]再分割成x[0], x[4] 和 x[2], x[6]

信号x[1], x[3], x[5], x[7]再分割成x[1], x[5] 和 x[3], x[7]

3. 无法分割了, 已经分割成2点了☺.

如上图:

在LayerI的时候, 我们是对2点进行DFT.( 一共4次DFT )

输入为 x[0]&x[4]; x[2]&x[6]; x[1]&x[5]; x[3]&x[7]

输出为 y[0],y[1]; Y[2],y[3]; Y[4],y[5]; Y[6],y[7];

流程:

I. 希望将输入直接转换为x[0], x[4], x[2], x[6], x[1], x[5], x[3], x[7]的顺序

II. 对转换顺序后的信号进行4次DFT

步骤I代码实现

/**

* 反转算法. 这个算法效率比较低!先用起来在说, 之后需要进行优化.

*/

static void bitrev( void )

{

int p=1, q, i;

int bit_rev[ N ];

float xx_r[ N ];

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q

{

bit_rev[ q ] = bit_rev[ q ] * 2;

bit_rev[ q + p ] = bit_rev[ q ] + 1;

}

p *= 2;

}

for(i=0; i

for(i=0; i

}

// ------------------------ 此刻序列x重排完毕------------------------

步骤II代码实现

int j;

float TR; // 临时变量

float tw1; // 旋转因子

/* 两点DFT */

for(k=0; k

{

// 两点DFT简化告诉我们tw1=1

TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

x_r[k] = TR + tw1*x_r[k+b];

x_r[k+b] = TR - tw1*x_r[k+b];

}

在LayerII的时候, 我们希望得到z, 就需要对y进行DFT.

y[0],y[2]; y[1],y[3]; y[4],y[6]; y[5],y[7];

z[0], z[1]; z[2],z[3]; z[4],z[5]; z[6],z[7];

在LayerIII的时候, 我们希望得到v, 就需要对z进行DFT.

z[0],z[4]; z[1],z[5]; z[2],z[6]; z[3],z[7];

v[0],v[1]; v[2],v[3]; v[4],v[5]; v[6],v[7];

准备

令输入为x[s], x[s+N/2], 输出为y[s], y[s+N/2]

这个N绝对不是上面的8, 这个N是当前颗粒的输入样本总量

对于LayerI而言N是2; 对于LayerII而言N是4; 对于LayerIII而言N是8

复数乘法:(a+j*b) * (c+j*d)

实部 = a*c – bd;

虚部 = ad + bc;

旋转因子:

实现(C描述)

#include

#include

#include

//#include "complex.h"

//

--------------------------------------------------------------------------

#define N 8 //64

#define M 3 //6 //2^m=N

#define PI 3.1415926

//

--------------------------------------------------------------------------

float twiddle[N/2] = {1.0, 0.707, 0.0, -0.707};

float x_r[N] = {1, 1, 1, 1, 0, 0, 0, 0};

float x_i[N]; //N=8

/*

float twiddle[N/2] = {1, 0.9951, 0.9808, 0.9570, 0.9239, 0.8820,

0.8317, 0.7733,

0.7075,

0.6349, 0.5561, 0.4721, 0.3835, 0.2912, 0.1961, 0.0991,

0.0000,-0.0991,-0.1961,-0.2912,-0.3835,-0.4721,-0.5561,-0.6349,

-0.7075,-0.7733,

0.8317,-0.8820,-0.9239,-0.9570,-0.9808,-0.9951}; //N=64

float x_r[N]={1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

1,1,1,1,1,1,1,1,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,

0,0,0,0,0,0,0,0,};

float x_i[N];

*/

FILE *fp;

// ----------------------------------- func

-----------------------------------

/**

* 初始化输出虚部

*/

static void fft_init( void )

{

int i;

for(i=0; i

}

/**

* 反转算法.将时域信号重新排序.

* 这个算法有改进的空间

*/

static void bitrev( void )

{

int p=1, q, i;

int bit_rev[ N ]; //

float xx_r[ N ]; //

bit_rev[ 0 ] = 0;

while( p < N )

{

for(q=0; q

{

bit_rev[ q ] = bit_rev[ q ] * 2;

bit_rev[ q + p ] = bit_rev[ q ] + 1;

}

p *= 2;

}

for(i=0; i

for(i=0; i

/* ------------ add by sshc625 ------------ */

static void bitrev2( void )

{

return ;

}

/* */

void display( void )

{

printf("/n/n");

int i;

for(i=0; i

printf("%f/t%f/n", x_r[i], x_i[i]);

}

/**

*

*/

void fft1( void )

{ fp = fopen("", "a+");

int L, i, b, j, p, k, tx1, tx2;

float TR, TI, temp; // 临时变量

float tw1, tw2;

/* 深M. 对层进行循环. L为当前层, 总层数为M. */

for(L=1; L<=M; L++)

{

fprintf(fp,"----------Layer=%d----------/n", L);

/* b的意义非常重大,b表示当前层的颗粒具有的输入样本点数 */

b = 1;

i = L - 1;

while(i > 0)

{

b *= 2;

i--;

}

// -------------- 是否外层对颗粒循环, 内层对样本点循环逻辑性更强一些呢!

--------------

/*

* outter对参与DFT的样本点进行循环

* L=1, 循环了1次(4个颗粒, 每个颗粒2个样本点)

* L=2, 循环了2次(2个颗粒, 每个颗粒4个样本点)

* L=3, 循环了4次(1个颗粒, 每个颗粒8个样本点)

*/

for(j=0; j

{

/* 求旋转因子tw1 */

p = 1;

i = M - L; // M是为总层数, L为当前层.

while(i > 0)

{

p = p*2;

i--;

}

p = p * j;

tx1 = p % N;

tx2 = tx1 + 3*N/4;

tx2 = tx2 % N;

// tw1是cos部分, 实部; tw2是sin部分, 虚数部分.

tw1 = ( tx1>=N/2)? -twiddle[tx1-N/2] : twiddle[ tx1 ];

tw2 = ( tx2>=N/2)? -twiddle[tx2-(N/2)] : twiddle[tx2];

/*

* inner对颗粒进行循环

* L=1, 循环了4次(4个颗粒, 每个颗粒2个输入)

* L=2, 循环了2次(2个颗粒, 每个颗粒4个输入)

* L=3, 循环了1次(1个颗粒, 每个颗粒8个输入)

*/

for(k=j; k

{

TR = x_r[k]; // TR就是A, x_r[k+b]就是B.

TI = x_i[k];

temp = x_r[k+b];

/*

* 如果复习一下 (a+j*b)(c+j*d)两个复数相乘后的实部虚部分别是什么

* 就能理解为什么会如下运算了, 只有在L=1时候输入才是实数, 之后层的

* 输入都是复数, 为了让所有的层的输入都是复数, 我们只好让L=1时候的

* 输入虚部为0

* x_i[k+b]*tw2是两个虚数相乘

*/

fprintf(fp, "tw1=%f, tw2=%f/n", tw1, tw2);

x_r[k] = TR + x_r[k+b]*tw1 + x_i[k+b]*tw2;

x_i[k] = TI - x_r[k+b]*tw2 + x_i[k+b]*tw1;

x_r[k+b] = TR - x_r[k+b]*tw1 - x_i[k+b]*tw2;

x_i[k+b] = TI + temp*tw2 - x_i[k+b]*tw1;

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f/n", k, x_r[k],

x_i[k]);

fprintf(fp, "k=%d, x_r[k]=%f, x_i[k]=%f/n", k+b, x_r[k+b],

x_i[k+b]);

} //

} //

} //

}

/**

* ------------ add by sshc625 ------------

* 该实现的流程为

* for( Layer )

* for( Granule )

* for( Sample )

*

*

*

*

*/

void fft2( void )

{ fp = fopen("", "a+");

int cur_layer, gr_num, i, k, p;

float tmp_real, tmp_imag, temp; // 临时变量, 记录实部

float tw1, tw2;// 旋转因子,tw1为旋转因子的实部cos部分, tw2为旋转因子的虚部sin部分.

int step; // 步进

int sample_num; // 颗粒的样本总数(各层不同, 因为各层颗粒的输入不同)

/* 对层循环 */

for(cur_layer=1; cur_layer<=M; cur_layer++)

{

/* 求当前层拥有多少个颗粒(gr_num) */

gr_num = 1;

i = M - cur_layer;

while(i > 0)

{

i--;

gr_num *= 2;

}

/* 每个颗粒的输入样本数N' */

sample_num = (int)pow(2, cur_layer);

/* 步进. 步进是N'/2 */

step = sample_num/2;

/* */

k = 0;

/* 对颗粒进行循环 */

for(i=0; i

{

/*

* 对样本点进行循环, 注意上限和步进

*/

for(p=0; p

{

// 旋转因子, 需要优化...

tw1 = cos(2*PI*p/pow(2, cur_layer));

tw2 = -sin(2*PI*p/pow(2, cur_layer));

tmp_real = x_r[k+p];

tmp_imag = x_i[k+p];

temp = x_r[k+p+step];

/*(tw1+jtw2)(x_r[k]+jx_i[k])

*

* real : tw1*x_r[k] - tw2*x_i[k]

* imag : tw1*x_i[k] + tw2*x_r[k]

* 我想不抽象出一个

* typedef struct {

* double real; // 实部

* double imag; // 虚部

* } complex; 以及针对complex的操作

* 来简化复数运算是否是因为效率上的考虑!

*/

/* 蝶形算法 */

x_r[k+p] = tmp_real + ( tw1*x_r[k+p+step] -

tw2*x_i[k+p+step] );

x_i[k+p] = tmp_imag + ( tw2*x_r[k+p+step] +

tw1*x_i[k+p+step] );

/* X[k] = A(k)+WB(k)

* X[k+N/2] = A(k)-WB(k) 的性质可以优化这里*/

// 旋转因子, 需要优化...

tw1 = cos(2*PI*(p+step)/pow(2, cur_layer));

tw2 = -sin(2*PI*(p+step)/pow(2, cur_layer));

x_r[k+p+step] = tmp_real + ( tw1*temp -

tw2*x_i[k+p+step] );

x_i[k+p+step] = tmp_imag + ( tw2*temp +

tw1*x_i[k+p+step] );

printf("k=%d, x_r[k]=%f, x_i[k]=%f/n", k+p, x_r[k+p],

x_i[k+p]);

printf("k=%d, x_r[k]=%f, x_i[k]=%f/n", k+p+step,

x_r[k+p+step], x_i[k+p+step]);

}

/* 开跳!:) */

k += 2*step;

}

}

}

/*

* 后记:

* 究竟是颗粒在外层循环还是样本输入在外层, 好象也差不多, 复杂度完全一样.

* 但以我资质愚钝花费了不少时间才弄明白这数十行代码.

* 从中我发现一个于我非常有帮助的教训, 很久以前我写过一部分算法, 其中绝大多数都是递归.

* 将数据量减少, 减少再减少, 用归纳的方式来找出数据量加大代码的规律

* 比如FFT

* 1. 先写死LayerI的代码; 然后再把LayerI的输出作为LayerII的输入, 又写死代码; ......

* 大约3层就可以统计出规律来. 这和递归也是一样, 先写死一两层, 自然就出来了!

* 2. 有的功能可以写伪代码, 不急于求出结果, 降低复杂性, 把逻辑结果定出来后再添加.

* 比如旋转因子就可以写死, 就写1.0. 流程出来后再写旋转因子.

* 寥寥数语, 我可真是流了不少汗! Happy!

*/

void dft( void )

{

int i, n, k, tx1, tx2;

float tw1,tw2;

float xx_r[N],xx_i[N];

/*

* clear any data in Real and Imaginary result arrays prior to DFT

*/

for(k=0; k<=N-1; k++)

xx_r[k] = xx_i[k] = x_i[k] = 0.0;

// caculate the DFT

for(k=0; k<=(N-1); k++)

{

for(n=0; n<=(N-1); n++)

{

tx1 = (n*k);

tx2 = tx1+(3*N)/4;

tx1 = tx1%(N);

tx2 = tx2%(N);

if(tx1 >= (N/2))

tw1 = -twiddle[tx1-(N/2)];

else

tw1 = twiddle[tx1];

if(tx2 >= (N/2))

tw2 = -twiddle[tx2-(N/2)];

else

tw2 = twiddle[tx2];

xx_r[k] = xx_r[k]+x_r[n]*tw1;

xx_i[k] = xx_i[k]+x_r[n]*tw2;

}

xx_i[k] = -xx_i[k];

}

// display

for(i=0; i

printf("%f/t%f/n", xx_r[i], xx_i[i]);

}

//

---------------------------------------------------------------------------

int main( void )

{

fft_init( );

bitrev( );

// bitrev2( );

//fft1( );

fft2( );

display( );

system( "pause" );

// dft();

return 1;

}

#include

#include

/*********************************************************************

快速福利叶变换C函数

函数简介:此函数是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依

赖硬件。此函数采用联合体的形式表示一个复数,输入为自然顺序的复

数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的

复数

使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的

应该为2的N次方,不满足此条件时应在后面补0

函数调用:FFT(s);

时 间:2010-2-20

版 本:Ver1.0

参考文献:

**********************************************************************/

#include

#define PI 3.97932384626433832795028841971 //定义圆周率值

#define FFT_N 128 //定义福利叶变换的点数

struct compx {float real,imag;}; //定义一个复数结构

struct compx s[FFT_N]; //FFT输入和输出:从S[1]开始存放,根据大小自己定义

/*******************************************************************

函数原型:struct compx EE(struct compx b1,struct compx b2)

函数功能:对两个复数进行乘法运算

输入参数:两个以联合体定义的复数a,b

输出参数:a和b的乘积,以联合体的形式输出

*******************************************************************/

struct compx EE(struct compx a,struct compx b)

{

struct compx c;

=**;

=*+*;

return(c);

}

/*****************************************************************

函数原型:void FFT(struct compx *xin,int N)

函数功能:对输入的复数组进行快速傅里叶变换(FFT)

输入参数:*xin复数结构体组的首地址指针,struct型

*****************************************************************/

void FFT(struct compx *xin)

{

int f,m,nv2,nm1,i,k,l,j=0;

struct compx u,w,t;

nv2=FFT_N/2; //变址运算,即把自然顺序变成倒位序,采用雷德算法

nm1=FFT_N-1;

for(i=0;i

{

if(i

{

t=xin[j];

xin[j]=xin[i];

xin[i]=t;

}

k=nv2; //求j的下一个倒位序

while(k<=j) //如果k<=j,表示j的最高位为1

{

j=j-k; //把最高位变成0

k=k/2; //k/2,比较次高位,依次类推,逐个比较,直到某个位为0

}

j=j+k; //把0改为1

}

{

int le,lei,ip; //FFT运算核,使用蝶形运算完成FFT运算

f=FFT_N;

for(l=1;(f=f/2)!=1;l++) //计算l的值,即计算蝶形级数

;

for(m=1;m<=l;m++) // 控制蝶形结级数

{ //m表示第m级蝶形,l为蝶形级总数l=log(2)N

le=2<<(m-1); //le蝶形结距离,即第m级蝶形的蝶形结相距le点

lei=le/2; //同一蝶形结中参加运算的两点的距离

=1.0; //u为蝶形结运算系数,初始值为1

=0.0;

=cos(PI/lei); //w为系数商,即当前系数与前一个系数的商

=-sin(PI/lei);

for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结,即计算系数不同的蝶形结

{

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算,即计算系数相同蝶形结

{

ip=i+lei; //i,ip分别表示参加蝶形运算的两个节点

t=EE(xin[ip],u); //蝶形运算,详见公式

xin[ip].real=xin[i].;

xin[ip].imag=xin[i].;

xin[i].real=xin[i].real+;

xin[i].imag=xin[i].imag+;

}

u=EE(u,w); //改变系数,进行下一个蝶形运算

}

}

}

}

/************************************************************

函数原型:void main()

函数功能:测试FFT变换,演示函数使用方法

输入参数:无

输出参数:无

************************************************************/

void main()

{

int i;

for(i=0;i

{

s[i].real=sin(2*3.9793*i/FFT_N); //实部为正弦波FFT_N点采样,赋值为1

s[i].imag=0; //虚部为0

}

FFT(s); //进行快速福利叶变换

for(i=0;i

s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);

while(1);

}

#include

#include

/*********************************************************************

快速福利叶变换C程序包

函数简介:此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依

赖硬件。此程序包采用联合体的形式表示一个复数,输入为自然顺序的复

数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的

复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,

以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度

使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的

应该为2的N次方,不满足此条件时应在后面补0。若使用查表法计算sin值和

cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表

函数调用:FFT(s);

时 间:2010-2-20

版 本:Ver1.1

参考文献:

**********************************************************************/

#include

#define FFT_N 128 //定义福利叶变换的点数

#define PI 3.97932384626433832795028841971 //定义圆周率值

struct compx {float real,imag;}; //定义一个复数结构

struct compx s[FFT_N]; //FFT输入和输出:从S[0]开始存放,根据大小自己定义

float SIN_TAB[FFT_N/2]; //定义正弦表的存放空间

/*******************************************************************

函数原型:struct compx EE(struct compx b1,struct compx b2)

函数功能:对两个复数进行乘法运算

输入参数:两个以联合体定义的复数a,b

输出参数:a和b的乘积,以联合体的形式输出

*******************************************************************/

struct compx EE(struct compx a,struct compx b)

{

struct compx c;

=**;

=*+*;

return(c);

}

/******************************************************************

函数原型:void create_sin_tab(float *sin_t)

函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同

输入参数:*sin_t存放正弦表的数组指针

输出参数:无

******************************************************************/

void create_sin_tab(float *sin_t)

{

int i;

for(i=0;i

sin_t[i]=sin(2*PI*i/FFT_N);

}

/******************************************************************

函数原型:void sin_tab(float pi)

函数功能:采用查表的方法计算一个数的正弦值

输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换

输出参数:输入值pi的正弦值

******************************************************************/

float sin_tab(float pi)

{

int n;

float a;

n=(int)(pi*FFT_N/2/PI);

if(n>=0&&n

a=SIN_TAB[n];

else if(n>=FFT_N/2&&n

a=-SIN_TAB[n-FFT_N/2];

return a;

}

/******************************************************************

函数原型:void cos_tab(float pi)

函数功能:采用查表的方法计算一个数的余弦值

输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换

输出参数:输入值pi的余弦值

******************************************************************/

float cos_tab(float pi)

{

float a,pi2;

pi2=pi+PI/2;

if(pi2>2*PI)

pi2-=2*PI;

a=sin_tab(pi2);

return a;

}

/*****************************************************************

函数原型:void FFT(struct compx *xin,int N)

函数功能:对输入的复数组进行快速傅里叶变换(FFT)

输入参数:*xin复数结构体组的首地址指针,struct型

输出参数:无

*****************************************************************/

void FFT(struct compx *xin)

{

int f,m,nv2,nm1,i,k,l,j=0;

struct compx u,w,t;

nv2=FFT_N/2;

nm1=FFT_N-1;

for(i=0;i

{

if(i

{

t=xin[j];

xin[j]=xin[i];

xin[i]=t;

}

k=nv2;

while(k<=j)

{

j=j-k;

k=k/2;

}

j=j+k;

}

{

int le,lei,ip;

f=FFT_N;

for(l=1;(f=f/2)!=1;l++)

;

for(m=1;m<=l;m++)

{

le=2<<(m-1);

lei=le/2;

=1.0;

=0.0;

//=cos(PI/lei);

// =-sin(PI/lei);

=cos_tab(PI/lei);

=-sin_tab(PI/lei); //变址运算,即把自然顺序变成倒位序,采用雷德算法

//如果i

//求j的下一个倒位序

//如果k<=j,表示j的最高位为1

//把最高位变成0

//k/2,比较次高位,依次类推,逐个比较,直到某个位为0//把0改为1

//FFT运算核,使用蝶形运算完成FFT运算

//计算l的值,即计算蝶形级数

// 控制蝶形结级数

//m表示第m级蝶形,l为蝶形级总数l=log(2)N

//le蝶形结距离,即第m级蝶形的蝶形结相距le点

//同一蝶形结中参加运算的两点的距离

//u为蝶形结运算系数,初始值为1

//不适用查表法计算sin值和cos值

//w为系数商,即当前系数与前一个系数的商

for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结,即计算系数不同的蝶形结

{

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算,即计算系数相同蝶形结

{

ip=i+lei; //i,ip分别表示参加蝶形运算的两个节点

t=EE(xin[ip],u); //蝶形运算,详见公式

xin[ip].real=xin[i].;

xin[ip].imag=xin[i].;

xin[i].real=xin[i].real+;

xin[i].imag=xin[i].imag+;

}

u=EE(u,w); //改变系数,进行下一个蝶形运算

}

}

}

}

/************************************************************

函数原型:void main()

函数功能:测试FFT变换,演示函数使用方法

输入参数:无

输出参数:无

************************************************************/

void main()

{

int i;

create_sin_tab(SIN_TAB);

for(i=0;i

{

s[i].real=sin(2*3.9793*i/FFT_N); //实部为正弦波FFT_N点采样,赋值为1

s[i].imag=0; //虚部为0

}

FFT(s); //进行快速福利叶变换

for(i=0;i

s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);

while(1);

}

#include

#include

/*********************************************************************

快速福利叶变换C程序包

函数简介:此程序包是通用的快速傅里叶变换C语言函数,移植性强,以下部分不依

赖硬件。此程序包采用联合体的形式表示一个复数,输入为自然顺序的复

数(输入实数是可令复数虚部为0),输出为经过FFT变换的自然顺序的

复数.此程序包可在初始化时调用create_sin_tab()函数创建正弦函数表,

以后的可采用查表法计算耗时较多的sin和cos运算,加快可计算速度.与

Ver1.1版相比较,Ver1.2版在创建正弦表时只建立了1/4个正弦波的采样值,

相比之下节省了FFT_N/4个存储空间

使用说明:使用此函数只需更改宏定义FFT_N的值即可实现点数的改变,FFT_N的

应该为2的N次方,不满足此条件时应在后面补0。若使用查表法计算sin值和

cos值,应在调用FFT函数前调用create_sin_tab()函数创建正弦表

函数调用:FFT(s);

时 间:2010-2-20

版 本:Ver1.2

参考文献:

**********************************************************************/

#include

#define FFT_N 128 //定义福利叶变换的点数

#define PI 3.97932384626433832795028841971 //定义圆周率值

struct compx {float real,imag;}; //定义一个复数结构

struct compx s[FFT_N]; //FFT输入和输出:从S[0]开始存放,根据大小自己定义

float SIN_TAB[FFT_N/4+1]; //定义正弦表的存放空间

/*******************************************************************

函数原型:struct compx EE(struct compx b1,struct compx b2)

函数功能:对两个复数进行乘法运算

输入参数:两个以联合体定义的复数a,b

输出参数:a和b的乘积,以联合体的形式输出

*******************************************************************/

struct compx EE(struct compx a,struct compx b)

{

struct compx c;

=**;

=*+*;

return(c);

}

/******************************************************************

函数原型:void create_sin_tab(float *sin_t)

函数功能:创建一个正弦采样表,采样点数与福利叶变换点数相同

输入参数:*sin_t存放正弦表的数组指针

输出参数:无

******************************************************************/

void create_sin_tab(float *sin_t)

{

int i;

for(i=0;i<=FFT_N/4;i++)

sin_t[i]=sin(2*PI*i/FFT_N);

}

/******************************************************************

函数原型:void sin_tab(float pi)

函数功能:采用查表的方法计算一个数的正弦值

输入参数:pi 所要计算正弦值弧度值,范围0--2*PI,不满足时需要转换

输出参数:输入值pi的正弦值

******************************************************************/

float sin_tab(float pi)

{

int n;

float a;

n=(int)(pi*FFT_N/2/PI);

if(n>=0&&n<=FFT_N/4)

a=SIN_TAB[n];

else if(n>FFT_N/4&&n

{

n-=FFT_N/4;

a=SIN_TAB[FFT_N/4-n];

}

else if(n>=FFT_N/2&&n<3*FFT_N/4)

{

n-=FFT_N/2;

a=-SIN_TAB[n];

}

else if(n>=3*FFT_N/4&&n<3*FFT_N)

{

n=FFT_N-n;

a=-SIN_TAB[n];

}

return a;

}

/******************************************************************

函数原型:void cos_tab(float pi)

函数功能:采用查表的方法计算一个数的余弦值

输入参数:pi 所要计算余弦值弧度值,范围0--2*PI,不满足时需要转换

输出参数:输入值pi的余弦值

******************************************************************/

float cos_tab(float pi)

{

float a,pi2;

pi2=pi+PI/2;

if(pi2>2*PI)

pi2-=2*PI;

a=sin_tab(pi2);

return a;

}

/*****************************************************************

函数原型:void FFT(struct compx *xin,int N)

函数功能:对输入的复数组进行快速傅里叶变换(FFT)

输入参数:*xin复数结构体组的首地址指针,struct型

输出参数:无

*****************************************************************/

void FFT(struct compx *xin)

{

int f,m,nv2,nm1,i,k,l,j=0;

struct compx u,w,t;

nv2=FFT_N/2; //变址运算,即把自然顺序变成倒位序,采用雷德算法

nm1=FFT_N-1;

for(i=0;i

{

if(i

{

t=xin[j];

xin[j]=xin[i];

xin[i]=t;

}

k=nv2; //求j的下一个倒位序

while(k<=j) //如果k<=j,表示j的最高位为1

{

j=j-k; //把最高位变成0

k=k/2; //k/2,比较次高位,依次类推,逐个比较,直到某个位为0

}

j=j+k; //把0改为1

}

{

int le,lei,ip; //FFT运算核,使用蝶形运算完成FFT运算

f=FFT_N;

for(l=1;(f=f/2)!=1;l++) //计算l的值,即计算蝶形级数

;

for(m=1;m<=l;m++) // 控制蝶形结级数

{ //m表示第m级蝶形,l为蝶形级总数l=log(2)N

le=2<<(m-1); //le蝶形结距离,即第m级蝶形的蝶形结相距le点

lei=le/2; //同一蝶形结中参加运算的两点的距离

=1.0; //u为蝶形结运算系数,初始值为1

=0.0;

//=cos(PI/lei); //不适用查表法计算sin值和cos值

// =-sin(PI/lei);

=cos_tab(PI/lei); //w为系数商,即当前系数与前一个系数的商

=-sin_tab(PI/lei);

for(j=0;j<=lei-1;j++) //控制计算不同种蝶形结,即计算系数不同的蝶形结

{

for(i=j;i<=FFT_N-1;i=i+le) //控制同一蝶形结运算,即计算系数相同蝶形结

{

ip=i+lei; //i,ip分别表示参加蝶形运算的两个节点

t=EE(xin[ip],u); //蝶形运算,详见公式

xin[ip].real=xin[i].;

xin[ip].imag=xin[i].;

xin[i].real=xin[i].real+;

xin[i].imag=xin[i].imag+;

}

u=EE(u,w); //改变系数,进行下一个蝶形运算

}

}

}

}

/************************************************************

函数原型:void main()

函数功能:测试FFT变换,演示函数使用方法

输入参数:无

输出参数:无

************************************************************/

void main()

{

int i;

create_sin_tab(SIN_TAB);

for(i=0;i

{

s[i].real=sin(2*3.9793*i/FFT_N); //实部为正弦波FFT_N点采样,赋值为1

s[i].imag=0; //虚部为0

}

FFT(s);

for(i=0;i

部部分

s[i].real=sqrt(s[i].real*s[i].real+s[i].imag*s[i].imag);

while(1);

}

//进行快速福利叶变换

//求变换后结果的模值,存入复数的实

看到的跟大家分享一下。。。。

FFT是离散傅立叶变换的快速算法,可以将一个信号变换 //本身就离散有限,满足DFT条件

到频域。有些信号在时域上是很难看出什么特征的,但是如

果变换到频域之后,就很容易看出特征了。这就是很多信号

分析采用FFT变换的原因。另外,FFT可以将一个信号的频谱//不满足的通过抽样截断,提取一部分咯

提取出来,这在频谱分析方面也是经常用的。

虽然很多人都知道FFT是什么,可以用来做什么,怎么去

做,但是却不知道FFT之后的结果是什意思、如何决定要使用

多少点来做FFT。

现在圈圈就根据实际经验来说说FFT结果的具体物理意义。

一个模拟信号,经过ADC采样之后,就变成了数字信号。采样

定理告诉我们,采样频率要大于信号频率的两倍,这些我就

不在此罗嗦了。

采样得到的数字信号,就可以做FFT变换了。N个采样点,

经过FFT之后,就可以得到N个点的FFT结果。为了方便进行FFT

运算,通常N取2的整数次方。

假设采样频率为Fs,信号频率F,采样点数为N。那么FFT//截断时间可以以一个周期去理解

之后结果就是一个为N点的复数。每一个点就对应着一个频率

点。这个点的模值,就是该频率值下的幅度特性。具体跟原始

信号的幅度有什么关系呢?假设原始信号的峰值为A,那么FFT

的结果的每个点(除了第一个点直流分量之外)的模值就是A

的N/2倍。而第一个点就是直流分量,它的模值就是直流分量

的N倍。而每个点的相位呢,就是在该频率下的信号的相位。

第一个点表示直流分量(即0Hz),而最后一个点N的再下一个

点(实际上这个点是不存在的,这里是假设的第N+1个点,也

可以看做是将第一个点分做两半分,另一半移到最后)则表示

采样频率Fs,这中间被N-1个点平均分成N等份,每个点的频率

依次增加。例如某点n所表示的频率为:Fn=(n-1)*Fs/N。 //0到N-1(此处1到N)刚好N个点嘛第一个,第二个

由上面的公式可以看出,Fn所能分辨到频率为为Fs/N,如果

采样频率Fs为1024Hz,采样点数为1024点,则可以分辨到1Hz。

1024Hz的采样率采样1024点,刚好是1秒,也就是说,采样1秒//频域的分辨率1HZ转换过去时域的一个周期,或采样的时间

时间的信号并做FFT,则结果可以分析到1Hz,如果采样2秒时

间的信号并做FFT,则结果可以分析到0.5Hz。如果要提高频率

分辨力,则必须增加采样点数,也即采样时间。频率分辨率和//当采样间距一定时

采样时间是倒数关系。

假设FFT之后某点n用复数a+bi表示,那么这个复数的模就是

An=根号a*a+b*b,相位就是Pn=atan2(b,a)。根据以上的结果,

就可以计算出n点(n≠1,且n<=N/2)对应的信号的表达式为:

An/(N/2)*cos(2*pi*Fn*t+Pn),即2*An/N*cos(2*pi*Fn*t+Pn)。

对于n=1点的信号,是直流分量,幅度即为A1/N。 //抽样后FFT分析的转换公式而已,位置,整体波形不变整体大小有专门的转换公式

由于FFT结果的对称性,通常我们只使用前半部分的结果,//频域里的幅值跟时域正弦里的幅值是有联系的

即小于采样频率一半的结果。

好了,说了半天,看着公式也晕,下面圈圈以一个实际的

信号来做说明。

假设我们有一个信号,它含有2V的直流分量,频率为50Hz、//信号里往往有不同的频率成分

相位为-30度、幅度为3V的交流信号,以及一个频率为75Hz、

相位为90度、幅度为1.5V的交流信号。用数学表达式就是如下:

S=2+3*cos(2*pi*50*t-pi*30/180)+1.5*cos(2*pi*75*t+pi*90/180)

式中cos参数为弧度,所以-30度和90度要分别换算成弧度。

我们以256Hz的采样率对这个信号进行采样,总共采样256点。

按照我们上面的分析,Fn=(n-1)*Fs/N,我们可以知道,每两个

点之间的间距就是1Hz,第n个点的频率就是n-1。我们的信号 //采样频率,就是得到频域后分析的周期,一个周期的最大末端值0HZ到(n-1)HZ

有3个频率:0Hz、50Hz、75Hz,应该分别在第1个点、第51个点、//把一段时域信号,在频域只是点化了,各种各样的不同频率在时域混合在频域给分开了嘛。

第76个点上出现峰值,其它各点应该接近0。实际情况如何呢?//信号里只有这几个频率成分,没其他的嘛

我们来看看FFT的结果的模值如图所示。

图1 FFT结果

从图中我们可以看到,在第1点、第51点、和第76点附近有

比较大的值。我们分别将这三个点附近的数据拿上来细看:

1点: 512+0i

2点: -2.6195E-14 - 1.4162E-13i

3点: -2.8586E-14 - 1.1898E-13i

50点:-6.2076E-13 - 2.1713E-12i

51点:332.55 - 192i

52点:-1.6707E-12 - 1.5241E-12i

75点:-2.2199E-13 -1.0076E-12i

76点:3.4315E-12 + 192i

77点:-3.0263E-14 +7.5609E-13i

很明显,1点、51点、76点的值都比较大,它附近的点值

都很小,可以认为是0,即在那些频率点上的信号幅度为0。

接着,我们来计算各点的幅度值。分别计算这三个点的模值,

结果如下:

1点: 512

51点:384

76点:192

按照公式,可以计算出直流分量为:512/N=512/256=2;

50Hz信号的幅度为:384/(N/2)=384/(256/2)=3;75Hz信号的

幅度为192/(N/2)=192/(256/2)=1.5。可见,从频谱分析出来

的幅度是正确的。

然后再来计算相位信息。直流信号没有相位可言,不用管

它。先计算50Hz信号的相位,atan2(-192, 332.55)=-0.5236,

结果是弧度,换算为角度就是180*(-0.5236)/pi=-30.0001。再

计算75Hz信号的相位,atan2(192, 3.4315E-12)=1.5708弧度,

换算成角度就是180*1.5708/pi=90.0002。可见,相位也是对的。

根据FFT结果以及上面的分析计算,我们就可以写出信号的表达//通过FFT得出未知信号的信息

式了,它就是我们开始提供的信号。

总结:假设采样频率为Fs,采样点数为N,做FFT之后,某

一点n(n从1开始)表示的频率为:Fn=(n-1)*Fs/N;该点的模值

除以N/2就是对应该频率下的信号的幅度(对于直流信号是除以

N);该点的相位即是对应该频率下的信号的相位。相位的计算

可用函数atan2(b,a)计算。atan2(b,a)是求坐标为(a,b)点的角

度值,范围从-pi到pi。要精确到xHz,则需要采样长度为1/x秒

的信号,并做FFT。要提高频率分辨率,就需要增加采样点数,

这在一些实际的应用中是不现实的,需要在较短的时间内完成

分析。解决这个问题的方法有频率细分法,比较简单的方法是

采样比较短时间的信号,然后在后面补充一定数量的0,使其长度

达到需要的点数,再做FFT,这在一定程度上能够提高频率分辨力。//这就是为何扩充成2的幂级数

具体的频率细分法可参考相关文献。

[附录:本测试数据使用的matlab程序]

close all; %先关闭所有图片

Adc=2; %直流分量幅度

A1=3; %频率F1信号的幅度

A2=1.5; %频率F2信号的幅度

F1=50; %信号1频率(Hz)

F2=75; %信号2频率(Hz)

Fs=256; %采样频率(Hz)

P1=-30; %信号1相位(度)

P2=90; %信号相位(度)

N=256; %采样点数

t=[0:1/Fs:N/Fs]; %采样时刻

%信号

S=Adc+A1*cos(2*pi*F1*t+pi*P1/180)+A2*cos(2*pi*F2*t+pi*P2/180);

%显示原始信号

plot(S);

title('原始信号');

figure;

Y = fft(S,N); %做FFT变换

Ayy = (abs(Y)); %取模

plot(Ayy(1:N)); %显示原始的FFT模值结果

title('FFT 模值');

figure;

Ayy=Ayy/(N/2); %换算成实际的幅度

Ayy(1)=Ayy(1)/2;

F=([1:N]-1)*Fs/N; %换算成实际的频率值

plot(F(1:N/2),Ayy(1:N/2)); %显示换算后的FFT模值结果

title('幅度-频率曲线图');

figure;

Pyy=[1:N/2];

for i="1:N/2"

Pyy(i)=phase(Y(i)); %计算相位

Pyy(i)=Pyy(i)*180/pi; %换算为角度

end;

plot(F(1:N/2),Pyy(1:N/2)); %显示相位图

title('相位-频率曲线图');

看完这个你就明白谐波分析了。


本文标签: 信号 频率 输入 蝶形 采样