基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教

实现(C描述)

#include <stdioh>

#include <mathh>

#include <stdlibh>

//#include "complexh"

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

#define N 8 //64

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

#define PI 31415926

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

float twiddle[N/2] = {10, 0707, 00, -0707};

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

float x_i[N]; //N=8

/

float twiddle[N/2] = {1, 09951, 09808, 09570, 09239, 08820, 08317, 07733,

07075, 06349, 05561, 04721, 03835, 02912, 01961, 00991,

00000,-00991,-01961,-02912,-03835,-04721,-05561,-06349,

-07075,-07733, 08317,-08820,-09239,-09570,-09808,-09951}; //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<N; i++) x_i[i] = 00;

}

/

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

这个算法有改进的空间

/

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<p; q++)

{

bit_rev[ q ] = bit_rev[ q ] 2;

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

}

p = 2;

}

for(i=0; i<N; i++) xx_r[ i ] = x_r[ i ];

for(i=0; i<N; i++) x_r[i] = xx_r[ bit_rev[i] ];

}

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

static void bitrev2( void )

{

return ;

}

/ /

void display( void )

{

printf("\n\n");

int i;

for(i=0; i<N; i++)

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

}

/

/

void fft1( void )

{ fp = fopen("log1txt", "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<b; j++)

{

/ 求旋转因子tw1 /

p = 1;

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

while(i > 0)

{

p = p2;

i--;

}

p = p j;

tx1 = p % N;

tx2 = tx1 + 3N/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<N; k=k+2b)

{

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

TI = x_i[k];

temp = x_r[k+b];

/

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

就能理解为什么会如下运算了, 只有在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 + temptw2 - 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("log2txt", "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<gr_num; i++)

{

/

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

/

for(p=0; p<sample_num/2; p++)

{

// 旋转因子, 需要优化

tw1 = cos(2PIp/pow(2, cur_layer));

tw2 = -sin(2PIp/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 : tw1x_r[k] - tw2x_i[k]

imag : tw1x_i[k] + tw2x_r[k]

我想不抽象出一个

typedef struct {

double real; // 实部

double imag; // 虚部

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

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

/

/ 蝶形算法 /

x_r[k+p] = tmp_real + ( tw1x_r[k+p+step] - tw2x_i[k+p+step] );

x_i[k+p] = tmp_imag + ( tw2x_r[k+p+step] + tw1x_i[k+p+step] );

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

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

// 旋转因子, 需要优化

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

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

x_r[k+p+step] = tmp_real + ( tw1temp - tw2x_i[k+p+step] );

x_i[k+p+step] = tmp_imag + ( tw2temp + tw1x_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 += 2step;

}

}

}

/

后记:

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

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

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

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

比如FFT

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

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

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

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

寥寥数语, 我可真是流了不少汗! 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] = 00;

// caculate the DFT

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

{

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

{

tx1 = (nk);

tx2 = tx1+(3N)/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<N; 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;

}

本文来自CSDN博客，转载请标明出处：>

两个问题：

（1）函数名fft与matlab自带的fft函数冲突，建议改名后重新存盘

（2）最好缺两end

改后的程序：

1),stem((0:7),abs(hc(1:8)));

grid;

title('8点正三角波的幅频特性');

hd(1:8)=fft(xd(1:8));

subplot(212),stem((0:7),abs(hd(1:8)));

grid;

title('8点反三角波的幅频特性');

elseif

aa==2

subplot(211),stem((0:15),xc(1:16));

title('8点正三角波补0到16点后的时间波形');

grid;

subplot(212),stem((0:15),xd(1:16));

title('8点反三角波补0到16点后的时间波形');

grid;

elseif

aa==3

hc(1:16)=fft(xc(1:16));

subplot(211),stem((0:15),abs(hc(1:16)));

title('8点正三角波补0到16点后的幅频特性');

end

end

function X=myfft(x)

%myfft函数 用递归实现

N=length(x);

t=log2(N);

t1=floor(t);

t2=ceil(t);

if t1~=t2; %若x的长度N不为2的整数次幂，则补0至最接近的2的整数次幂

x=[x zeros(1,2^t2-N)];

N=2^t2;

end

w0=exp(-j2pi/N);

X=zeros(1,N);

if N==2

X(1)=x(1)+x(2);

X(2)=x(1)-x(2);

else

n=1:N/2;

xe(n)=x(2n-1);

xo(n)=x(2n);

XE=myfft(xe); %递归调用

XO=myfft(xo);

for n=1:N/2

X(n)=XE(n)+XO(n)(w0^(n-1));

X(n+N/2)=XE(n)-XO(n)(w0^(n-1));

end

end

先输入x的值，如x=[1，2，3，4，5，6，7，8]

将文件名改为myfft，再调用myfft(x)，比较与fft（x）的结果。

Y(1:halfLength+1)中1:halfLength+1是索引，而索引必须是正整数，所以，不能从0开始，要都加1，但是指却是从0值开始到最后一个值。

至于f的算法是不一样的，在f

=((0:halfLength)+1)Fs1/n这里，0:halfLength)+1是数值，不是索引，结果是一个矢量。

以上就是关于基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教全部的内容，包括:基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教、西储大学轴承数据fft程序问题、请各位看一下下面的matlab 7.0下关于fft算法程序有什么问题，查不出错误，请高手帮帮忙非常感谢！等相关内容解答，如果想了解更多相关内容，可以关注我们，你们的支持是我们更新的动力！