本文目录一览：
- 用c语言实现FFT
- matlab fft2的c代码
- 基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教)
- 怎样用C语言实现FFT算法啊？
- 求FFT的C语言程序……最好是1024点的……希望大家帮帮我！
- C语言 1024点快速傅里叶变换（FFT）程序，最好经过优化，执行速度快
用c语言实现FFT
float ar[1024], ai[1024]; /* 原始数据实部，虚部 */
float a[2050];
void fft(int nn) /* nn数据长度 */
{
int n1, n2, i, j, k, l, m, s, l1;
float t1, t2, x, y;
float w1, w2, u1, u2, z;
float fsin[10] = {0.000000, 1.000000, 0.707107, 0.3826834, 0.1950903, 0.09801713, 0.04906767, 0.02454123, 0.01227154, 0.00613588};
float fcos[10] = {-1.000000, 0.000000, 0.7071068, 0.9238796, 0.9807853, 0.99518472, 0.99879545, 0.9996988, 0.9999247, 0.9999812};
switch (nn) {
case 1024: s = 10; break;
case 512: s = 9; break;
case 256: s = 8; break;
}
n1 = nn / 2;
n2 = nn - 1;
j = 1;
for (i = 1; i <= nn; i++) {
a[2 * i] = ar[i - 1];
a[2 * i + 1] = ai[i - 1];
}
for (l = 1; l <= n2; l++) {
if (l < j) {
t1 = a[2 * j];
t2 = a[2 * j + 1];
a[2 * j] = a[2 * l];
a[2 * j + 1] = a[2 * l + 1];
a[2 * l] = t1;
a[2 * l + 1] = t2;
}
k = n1;
while (k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
}
for (i = 1; i <= s; i++) {
u1 = 1;
u2 = 0;
m = (1 << i);
k = m >> 1;
w1 = fcos[i - 1];
w2 = -fsin[i - 1];
for (j = 1; j <= k; j++) {
for (l = j; l <= nn; l = l + m) {
l1 = l + k;
t1 = a[2 * l1] * u1 - a[2 * l1 + 1] * u2;
t2 = a[2 * l1] * u2 + a[2 * l1 + 1] * u1;
a[2 * l1] = a[2 * l] - t1;
a[2 * l1 + 1] = a[2 * l + 1] - t2;
a[2 * l] = a[2 * l] + t1;
a[2 * l + 1] = a[2 * l + 1] + t2;
}
z = u1 * w1 - u2 * w2;
u2 = u1 * w2 + u2 * w1;
u1 = z;
}
}
for (i = 1; i <= nn / 2; i++) {
ar[i] = 4 * a[2 * i + 2] / nn; /* 实部 */
ai[i] = -4 * a[2 * i + 3] / nn; /* 虚部 */
a[i] = 4 * sqrt(ar[i] * ar[i] + ai[i] * ai[i]); /* 幅值 */
}
}
打字不易，如满意，望采纳。
matlab fft2的c代码
傅立叶变换的c语言源代码 128点DIT FFT函数：
/* 采样来的数据放在dataR[]数组中，运算前dataI[]数组初始化为0 */
void FFT(float dataR[], float dataI[]) {
int x0, x1, x2, x3, x4, x5, x6;
int L, j, k, b, p;
float TR, TI, temp;
/********** following code invert sequence ************/
for (i = 0; i < 128; i++) {
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 < 128; i++) {
dataR[i] = dataI[i];
dataI[i] = 0;
}
/************** following code FFT *******************/
for (L = 1; L <= 7; L++) { /* for(1) */
b = 1;
i = L - 1;
while (i > 0) {
b = b * 2;
i--;
} /* b= 2^(L-1) */
for (j = 0; j < b; 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 < 32; i++) { /* 只需要32次以下的谐波进行分析 */
w[i] = sqrt(dataR[i] * dataR[i] + dataI[i] * dataI[i]);
w[i] = w[i] / 64;
}
w[0] = w[0] / 2;
} /* END FFT */
基于FFT的算法优化 要C语言完整程序（利用旋转因子的性质），有的请留言，答谢！！！(有核心代码，望指教)
实现(C描述)
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
//#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 < N; i++) x_i[i] = 0.0;
}
/**
* 反转算法.将时域信号重新排序.
* 这个算法有改进的空间
*/
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("log1.txt", "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 = 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] : -twiddle[tx1 - N / 2];
tw2 = (tx2 < N / 2) ? twiddle[tx2] : -twiddle[tx2 - N / 2];
/*
* inner对颗粒进行循环
* L=1, 循环了4次(4个颗粒, 每个颗粒2个输入)
* L=2, 循环了2次(2个颗粒, 每个颗粒4个输入)
* L=3, 循环了1次(1个颗粒, 每个颗粒8个输入)
*/
for (k = j; k < N; k = 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部分.
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(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 < 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博客，转载请标明出处：
怎样用C语言实现FFT算法啊？
- 二维FFT相当于对行和列分别进行一维FFT运算。具体的实现办法如下： 先对各行逐一进行一维FFT，然后再对变换后的新矩阵的各列逐一进行一维FFT。相应的伪代码如下所示：
for (int i = 0; i < M; i++)
FFT_1D(ROW[i], N);
for (int j = 0; j < N; j++)
FFT_1D(COL[j], M);
其中，ROW[i]表示矩阵的第i行。注意这只是一个简单的记法，并不能完全照抄。还需要通过一些语句来生成各行的数据。同理，COL[i]是对矩阵的第i列的一种简单表示方法。 所以，关键是一维FFT算法的实现。 2. 例程：
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#define N 1000
/*定义复数类型*/
typedef struct {
double real;
double img;
} complex;
complex x[N], *W; /*输入序列,变换核*/
int size_x = 0; /*输入序列的大小，在本程序中仅限2的次幂*/
double PI; /*圆周率*/
void fft(); /*快速傅里叶变换*/
void initW(); /*初始化变换核*/
void change(); /*变址*/
void add(complex, complex, complex*); /*复数加法*/
void mul(complex, complex, complex*); /*复数乘法*/
void sub(complex, complex, complex*); /*复数减法*/
void output();
int main() {
int i;
system("cls");
PI = atan(1) * 4;
printf("Please input the size of x:\n");
scanf("%d", &size_x);
printf("Please input the data in x[N]:\n");
for (i = 0; i < size_x; i++)
scanf("%lf%lf", &x[i].real, &x[i].img);
initW();
fft();
output();
return 0;
}
/*快速傅里叶变换*/
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 initW() {
int i;
W = (complex*)malloc(sizeof(complex) * size_x);
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;
}
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) {
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->real = a.real - b.real;
c->img = a.img - b.img;
}
求FFT的C语言程序……最好是1024点的……希望大家帮帮我！
float ar[1024], ai[1024]; /* 实部，虚部 */
float a[2050]; /* 实际值 */
void fft() {
int n1, n2, i, j, k, l, m, s = 10, nn = 1024, l1;
float t1, t2, x, y;
float w1, w2, u1, u2, z;
float fsin[10] = {0.000000, 1.000000, 0.707107, 0.3826834, 0.1950903, 0.09801713, 0.04906767, 0.02454123, 0.01227154, 0.00613588};
float fcos[10] = {-1.000000, 0.000000, 0.7071068, 0.9238796, 0.9807853, 0.99518472, 0.99879545, 0.9996988, 0.9999247, 0.9999812};
n1 = nn / 2;
n2 = nn - 1;
j = 1;
for (i = 1; i <= nn; i++) {
a[2 * i] = ar[i - 1];
a[2 * i + 1] = ai[i - 1];
}
for (l = 1; l <= n2; l++) {
if (l < j) {
t1 = a[2 * j];
t2 = a[2 * j + 1];
a[2 * j] = a[2 * l];
a[2 * j + 1] = a[2 * l + 1];
a[2 * l] = t1;
a[2 * l + 1] = t2;
}
k = n1;
while (k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
}
for (i = 1; i <= s; i++) {
u1 = 1;
u2 = 0;
m = (1 << i);
k = m >> 1;
w1 = fcos[i - 1];
w2 = -fsin[i - 1];
for (j = 1; j <= k; j++) {
for (l = j; l < nn; l = l + m) {
l1 = l + k;
t1 = a[2 * l1] * u1 - a[2 * l1 + 1] * u2;
t2 = a[2 * l1] * u2 + a[2 * l1 + 1] * u1;
a[2 * l1] = a[2 * l] - t1;
a[2 * l1 + 1] = a[2 * l + 1] - t2;
a[2 * l] = a[2 * l] + t1;
a[2 * l + 1] = a[2 * l + 1] + t2;
}
z = u1 * w1 - u2 * w2;
u2 = u1 * w2 + u2 * w1;
u1 = z;
}
}
for (i = 1; i <= nn / 2; i++) {
ar[i] = a[2 * i + 2] / nn;
ai[i] = -a[2 * i + 3] / nn;
a[i] = 4 * sqrt(ar[i] * ar[i] + ai[i] * ai[i]);
}
}
C语言 1024点快速傅里叶变换（FFT）程序，最好经过优化，执行速度快
void fft() {
int nn, n1, n2, i, j, k, l, m, s, l1;
float ar[1024], ai[1024]; // 实部 虚部
float a[2050];
float t1, t2, x, y;
float w1, w2, u1, u2, z;
float fsin[10] = {0.000000, 1.000000, 0.707107, 0.3826834, 0.1950903, 0.09801713, 0.04906767, 0.02454123, 0.01227154, 0.00613588}; // 优化
float fcos[10] = {-1.000000, 0.000000, 0.7071068, 0.9238796, 0.9807853, 0.99518472, 0.99879545, 0.9996988, 0.9999247, 0.9999812};
nn = 1024;
s = 10;
n1 = nn / 2;
n2 = nn - 1;
j = 1;
for (i = 1; i <= nn; i++) {
a[2 * i] = ar[i - 1];
a[2 * i + 1] = ai[i - 1];
}
for (l = 1; l <= n2; l++) {
if (l < j) {
t1 = a[2 * j];
t2 = a[2 * j + 1];
a[2 * j] = a[2 * l];
a[2 * j + 1] = a[2 * l + 1];
a[2 * l] = t1;
a[2 * l + 1] = t2;
}
k = n1;
while (k < j) {
j = j - k;
k = k / 2;
}
j = j + k;
}
for (i = 1; i <= s; i++) {
u1 = 1;
u2 = 0;
m = (1 << i);
k = m >> 1;
w1 = fcos[i - 1];
w2 = -fsin[i - 1];
for (j = 1; j <= k; j++) {
for (l = j; l < nn; l = l + m) {
l1 = l + k;
t1 = a[2 * l1] * u1 - a[2 * l1 + 1] * u2;
t2 = a[2 * l1] * u2 + a[2 * l1 + 1] * u1;
a[2 * l1] = a[2 * l] - t1;
a[2 * l1 + 1] = a[2 * l + 1] - t2;
a[2 * l] = a[2 * l] + t1;
a[2 * l + 1] = a[2 * l + 1] + t2;
}
z = u1 * w1 - u2 * w2;
u2 = u1 * w2 + u2 * w1;
u1 = z;
}
}
for (i = 1; i <= nn / 2; i++) {
ar[i] = a[2 * i + 2] / nn;
ai[i] = -a[2 * i + 3] / nn;
a[i] = 4 * sqrt(ar[i] * ar[i] + ai[i] * ai[i]); // 幅值
}
}
