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

实现(C描述)

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

//#include "complex.h"

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

#define N 8 //64

#defineM 3 //6 //2^m=N

#definePI3.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=0i<Ni++) x_i[i] = 0.0

}

/**

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

* 这个算法有改进的空间

*/

static void bitrev( void )

{

intp=1, q, i

intbit_rev[ N ]//

float xx_r[ N ] //

bit_rev[ 0 ] = 0

while( p <N )

{

for(q=0q<pq++)

{

bit_rev[ q ] = bit_rev[ q ] * 2

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

}

p *= 2

}

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

for(i=0i<Ni++) 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=0i<Ni++)

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

}

/**

*

*/

void fft1( void )

{ fp = fopen("log1.txt", "a+")

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

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

float tw1, tw2

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

for(L=1L<=ML++)

{

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=0j<bj++)

{

/* 求旋转因子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=jk<Nk=k+2*b)

{

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("log2.txt", "a+")

int cur_layer, gr_num, i, k, p

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

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

intstep // 步进

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

/* 对层循环 */

for(cur_layer=1cur_layer<=Mcur_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=0i<gr_numi++)

{

/*

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

*/

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

{

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

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 )

{

inti, 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=0k<=N-1k++)

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

// caculate the DFT

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

{

for(n=0n<=(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=0i<Ni++)

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博客，转载请标明出处：http://blog.csdn.net/sshcx/archive/2007/06/14/1651616.aspx

我给你一个标准遗传算法程序供你参考：

该程序是遗传算法优化BP神经网络函数极值寻优：

%% 该代码为基于神经网络遗传算法的系统极值寻优

%% 清空环境变量

clc

clear

%% 初始化遗传算法参数

%初始化参数

maxgen=100%进化代数，即迭代次数

sizepop=20 %种群规模

pcross=[0.4] %交叉概率选择，0和1之间

pmutation=[0.2] %变异概率选择，0和1之间

lenchrom=[1 1] %每个变量的字串长度，如果是浮点变量，则长度都为1

bound=[-5 5-5 5] %数据范围

individuals=struct('fitness',zeros(1,sizepop), 'chrom',[]) %将种群信息定义为一个结构体

avgfitness=[] %每一代种群的平均适应度

bestfitness=[]%每一代种群的最佳适应度

bestchrom=[] %适应度最好的染色体

%% 初始化种群计算适应度值

% 初始化种群

for i=1:sizepop

%随机产生一个种群

individuals.chrom(i,:)=Code(lenchrom,bound)

x=individuals.chrom(i,:)

%计算适应度

individuals.fitness(i)=fun(x) %染色体的适应度

end

%找最好的染色体

[bestfitness bestindex]=min(individuals.fitness)

bestchrom=individuals.chrom(bestindex,:) %最好的染色体

avgfitness=sum(individuals.fitness)/sizepop%染色体的平均适应度

% 记录每一代进化中最好的适应度和平均适应度

trace=[avgfitness bestfitness]

%% 迭代寻优

% 进化开始

for i=1:maxgen

i

% 选择

individuals=Select(individuals,sizepop)

avgfitness=sum(individuals.fitness)/sizepop

%交叉

individuals.chrom=Cross(pcross,lenchrom,individuals.chrom,sizepop,bound)

% 变异

individuals.chrom=Mutation(pmutation,lenchrom,individuals.chrom,sizepop,[i maxgen],bound)

% 计算适应度

for j=1:sizepop

x=individuals.chrom(j,:)%解码

individuals.fitness(j)=fun(x)

end

%找到最小和最大适应度的染色体及它们在种群中的位置

[newbestfitness,newbestindex]=min(individuals.fitness)

[worestfitness,worestindex]=max(individuals.fitness)

% 代替上一次进化中最好的染色体

if bestfitness>newbestfitness

bestfitness=newbestfitness

bestchrom=individuals.chrom(newbestindex,:)

end

individuals.chrom(worestindex,:)=bestchrom

individuals.fitness(worestindex)=bestfitness

avgfitness=sum(individuals.fitness)/sizepop

trace=[traceavgfitness bestfitness]%记录每一代进化中最好的适应度和平均适应度

end

%进化结束

%% 结果分析

[r c]=size(trace)

plot([1:r]',trace(:,2),'r-')

title('适应度曲线','fontsize',12)

xlabel('进化代数','fontsize',12)ylabel('适应度','fontsize',12)

axis([0,100,0,1])

disp('适应度 变量')

x=bestchrom

% 窗口显示

disp([bestfitness x])

这是一个POA算法（C#）的梯级电站优化计算的程序，来源于百度文库，代码太长，这里只给个连接。

http://wenku.baidu.com/link?url=QsudfgokAexO8mBHaIzz1nz3xfO1Fgn_PApx1ZMIvhq09hGN99vW5AHj9c8NgsmD1Pr6so08l6_hjQbHSgEDNe_t882vty8EFstMni1VzBe

---