我自己写的1024点CFFT小工具;
此CFFT算法为第一版，在算法模型上可以继续优化运算速度，生成的程序段长度也可以继续精简
; 但是我保证了在我能轻松愉快地编程的能力范围内，这是兼顾了运算精度与运算速度的最佳版本，如果要继续优化代码的话那我感觉就很难受
了; 下面是工具的测试和使用例，以及头文件：
1. #define _CRT_SECURE_NO_WARNINGS
#include <stdio.h>
#include "cfft.h"
#define POINTS 1024
#if POINTS == 4
#define FFT(a,b,c,d,e) fft_4(a,b,c,d)
#define IFFT(a,b,c,d,e) ifft_4(a,b,c,d)
#elif POINTS == 16
#define FFT(a,b,c,d,e) fft_16(a,b,c,d,e)
#define IFFT(a,b,c,d,e) ifft_16(a,b,c,d,e)
#elif POINTS == 64
#define FFT(a,b,c,d,e) fft_64(a,b,c,d,e)
#define IFFT(a,b,c,d,e) ifft_64(a,b,c,d,e)
#elif POINTS == 256
#define FFT(a,b,c,d,e) fft_256(a,b,c,d,e)
#define IFFT(a,b,c,d,e) ifft_256(a,b,c,d,e)
#elif POINTS == 1024
#define FFT(a,b,c,d,e) fft_1024(a,b,c,d,e)
#define IFFT(a,b,c,d,e) ifft_1024(a,b,c,d,e)
#endif int main() {
float* re1, * im1, * re2, * im2;
int i;
re1 = (float*)malloc(sizeof(float) * POINTS);
im1 = (float*)malloc(sizeof(float) * POINTS);
re2 = (float*)malloc(sizeof(float) * POINTS);
im2 = (float*)malloc(sizeof(float) * POINTS);
for (i = 0; i < POINTS; i++)
{
re1[i] = (float)i * 10;
im1[i] = 0;
}
for (i = 0; i < POINTS; i++)
{
printf("%.6f\t%.6f\n", re1[i], im1[i]);
}
FFT(re1, im1, re2, im2, 1); //if the last parameter is 1,
//then the FFT function will clean up
//the memory that it occupies at the end
printf("fft:\n");
for (i = 0; i < POINTS; i++)
{
printf("%.6f\t%.6f\n", re2[i], im2[i]);
}
IFFT(re2, im2, re1, im1, 1);
printf("ifft(/%d):\n",POINTS);
for (i = 0; i < POINTS; i++)
{
printf("%.6f\t%.6f\n", re1[i]/POINTS, im1[i]/POINTS);
}
free(re1);
free(im1);
free(re2);
free(im2);
return 0; }
头文件：
1. #ifndef _CFFT_TOOL_
#define _CFFT_TOOL_
#include <math.h>
#include <stdlib.h>
#define CFFT_PI 3.1415926535897
#define COEF_NUM_16 (3*3+1)
#define COEF_NUM_64 (15*3+1)
#define COEF_NUM_256 (63*3+1)
#define COEF_NUM_1024 (255*3+1) //declare void fft_4(
float* re_in,
float* im_in,
float* re_out,
float* im_out ); void ifft_4(
float* re_in,
float* im_in,
float* re_out,
float* im_out ); void fft_16(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void addr_sort(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int total ); void fft_coef_gen_16(
float* re_coef,
float* im_coef ); void fft_layer_transform(
float* re_out,
float* im_out,
float* re_coef,
float* im_coef,
int total ); void ifft_16(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void ifft_coef_gen_16(
float* re_coef,
float* im_coef ); void ifft_layer_transform(
float* re_out,
float* im_out,
float* re_coef,
float* im_coef,
int total ); void fft_64(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void fft_coef_gen_64(
float* re_coef,
float* im_coef ); void ifft_coef_gen_64(
float* re_coef,
float* im_coef ); void fft_256(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void ifft_256(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void ifft_coef_gen_256(
float* re_coef,
float* im_coef ); void fft_coef_gen_256(
float* re_coef,
float* im_coef ); void fft_1024(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void ifft_1024(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ); void ifft_coef_gen_1024(
float* re_coef,
float* im_coef ); void fft_coef_gen_1024(
float* re_coef,
float* im_coef ); //define void fft_4(
float* re_in,
float* im_in,
float* re_out,
float* im_out ) {
int i;
for (i = 0; i < 4; i++)
{
if (i == 0)
{
re_out[0] = re_in[0];
im_out[0] = im_in[0];
re_out[1] = re_in[0];
im_out[1] = im_in[0];
re_out[2] = re_in[0];
im_out[2] = im_in[0];
re_out[3] = re_in[0];
im_out[3] = im_in[0];
}
else if (i == 1)
{
re_out[0] += re_in[1];
im_out[0] += im_in[1];
re_out[1] += im_in[1];
im_out[1] -= re_in[1];
re_out[2] -= re_in[1];
im_out[2] -= im_in[1];
re_out[3] -= im_in[1];
im_out[3] += re_in[1];
}
else if (i == 2)
{
re_out[0] += re_in[2];
im_out[0] += im_in[2];
re_out[1] -= re_in[2];
im_out[1] -= im_in[2];
re_out[2] += re_in[2];
im_out[2] += im_in[2];
re_out[3] -= re_in[2];
im_out[3] -= im_in[2];
}
else if (i == 3)
{
re_out[0] += re_in[3];
im_out[0] += im_in[3];
re_out[1] -= im_in[3];
im_out[1] += re_in[3];
re_out[2] -= re_in[3];
im_out[2] -= im_in[3];
re_out[3] += im_in[3];
im_out[3] -= re_in[3];
}
} } void ifft_4(
float* re_in,
float* im_in,
float* re_out,
float* im_out ) {
int i;
for (i = 0; i < 4; i++)
{
if (i == 0)
{
re_out[0] = re_in[0];
im_out[0] = im_in[0];
re_out[1] = re_in[0];
im_out[1] = im_in[0];
re_out[2] = re_in[0];
im_out[2] = im_in[0];
re_out[3] = re_in[0];
im_out[3] = im_in[0];
}
else if (i == 1)
{
re_out[0] += re_in[1];
im_out[0] += im_in[1];
re_out[1] -= im_in[1];
im_out[1] += re_in[1];
re_out[2] -= re_in[1];
im_out[2] -= im_in[1];
re_out[3] += im_in[1];
im_out[3] -= re_in[1];
}
else if (i == 2)
{
re_out[0] += re_in[2];
im_out[0] += im_in[2];
re_out[1] -= re_in[2];
im_out[1] -= im_in[2];
re_out[2] += re_in[2];
im_out[2] += im_in[2];
re_out[3] -= re_in[2];
im_out[3] -= im_in[2];
}
else if (i == 3)
{
re_out[0] += re_in[3];
im_out[0] += im_in[3];
re_out[1] += im_in[3];
im_out[1] -= re_in[3];
re_out[2] -= re_in[3];
im_out[2] -= im_in[3];
re_out[3] -= im_in[3];
im_out[3] += re_in[3];
}
} } void fft_16(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_16);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_16);
re_tmp = (float*)malloc(sizeof(float) * 16);
im_tmp = (float*)malloc(sizeof(float) * 16);
fft_coef_gen_16(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
16
);
for (int i = 0; i < 16; i += 4)
{
fft_4(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i
);
}
fft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
4
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
start_flag = 1;
} } void addr_sort(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int total ) {
int i, j;
int* re_ptr = re_out, * im_ptr = im_out;
for (i = 0; i < 4; i++)
{
for (j = 0; j < total; j++)
{
if (j % 4 == i)
{
*re_out++ = re_in[j];
*im_out++ = im_in[j];
}
}
} } void fft_coef_gen_16(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_16; i++)
{
re_coef[i] = (float)cos(-2 * CFFT_PI * i / 16);
im_coef[i] = (float)sin(-2 * CFFT_PI * i / 16);
} } void fft_layer_transform(
float* re_out,
float* im_out,
float* re_coef,
float* im_coef,
int total ) {
float* re_first_ptr, * im_first_ptr;
float* re_second_ptr, * im_second_ptr;
float* re_coef_ptr, * im_coef_ptr;
float re_4_in[4];
float im_4_in[4];
float re_4_out[4];
float im_4_out[4];
float ac, bd, ad, bc;
int i, j;
re_first_ptr = re_out;
im_first_ptr = im_out;
for (i = 0; i < total; i++)
{
re_second_ptr = re_first_ptr;
im_second_ptr = im_first_ptr;
re_coef_ptr = re_coef;
im_coef_ptr = im_coef;
for (j = 0; j < 4; j++)
{
ac = (*re_second_ptr) * (*re_coef_ptr);
bd = (*im_second_ptr) * (*im_coef_ptr);
ad = (*re_second_ptr) * (*im_coef_ptr);
bc = (*im_second_ptr) * (*re_coef_ptr);
re_4_in[j] = ac - bd;
im_4_in[j] = ad + bc;
re_second_ptr += total;
im_second_ptr += total;
re_coef_ptr += i;
im_coef_ptr += i;
}
fft_4(
re_4_in,
im_4_in,
re_4_out,
im_4_out
);
for (j = 3; j >= 0; j--)
{
re_second_ptr -= total;
im_second_ptr -= total;
*re_second_ptr = re_4_out[j];
*im_second_ptr = im_4_out[j];
}
re_first_ptr++;
im_first_ptr++;
} } void ifft_16(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_16);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_16);
re_tmp = (float*)malloc(sizeof(float) * 16);
im_tmp = (float*)malloc(sizeof(float) * 16);
ifft_coef_gen_16(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation process
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
16
);
for (int i = 0; i < 16; i += 4)
{
ifft_4(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i
);
}
ifft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
4
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
start_flag = 1;
} } void ifft_coef_gen_16(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_16; i++)
{
re_coef[i] = (float)cos(2 * CFFT_PI * i / 16);
im_coef[i] = (float)sin(2 * CFFT_PI * i / 16);
} } void ifft_layer_transform(
float* re_out,
float* im_out,
float* re_coef,
float* im_coef,
int total ) {
float* re_first_ptr, * im_first_ptr;
float* re_second_ptr, * im_second_ptr;
float* re_coef_ptr, * im_coef_ptr;
float re_4_in[4];
float im_4_in[4];
float re_4_out[4];
float im_4_out[4];
float ac, bd, ad, bc;
int i, j;
re_first_ptr = re_out;
im_first_ptr = im_out;
for (i = 0; i < total; i++)
{
re_second_ptr = re_first_ptr;
im_second_ptr = im_first_ptr;
re_coef_ptr = re_coef;
im_coef_ptr = im_coef;
for (j = 0; j < 4; j++)
{
ac = (*re_second_ptr) * (*re_coef_ptr);
bd = (*im_second_ptr) * (*im_coef_ptr);
ad = (*re_second_ptr) * (*im_coef_ptr);
bc = (*im_second_ptr) * (*re_coef_ptr);
re_4_in[j] = ac - bd;
im_4_in[j] = ad + bc;
re_second_ptr += total;
im_second_ptr += total;
re_coef_ptr += i;
im_coef_ptr += i;
}
ifft_4(
re_4_in,
im_4_in,
re_4_out,
im_4_out
);
for (j = 3; j >= 0; j--)
{
re_second_ptr -= total;
im_second_ptr -= total;
*re_second_ptr = re_4_out[j];
*im_second_ptr = im_4_out[j];
}
re_first_ptr++;
im_first_ptr++;
} } void fft_64(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_64);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_64);
re_tmp = (float*)malloc(sizeof(float) * 64);
im_tmp = (float*)malloc(sizeof(float) * 64);
fft_coef_gen_64(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
64
);
for (int i = 0; i < 64; i += 16)
{
if (i == 48 && end_flag == 1)
{
start_flag = 1;
}
fft_16(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
fft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
16
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} } void fft_coef_gen_64(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_64; i++)
{
re_coef[i] = (float)cos(-2 * CFFT_PI * i / 64);
im_coef[i] = (float)sin(-2 * CFFT_PI * i / 64);
} } void ifft_64(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_64);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_64);
re_tmp = (float*)malloc(sizeof(float) * 64);
im_tmp = (float*)malloc(sizeof(float) * 64);
ifft_coef_gen_64(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
64
);
for (int i = 0; i < 64; i += 16)
{
if (i == 48 && end_flag == 1)
{
start_flag = 1;
}
ifft_16(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
ifft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
16
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} } void ifft_coef_gen_64(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_64; i++)
{
re_coef[i] = (float)cos(2 * CFFT_PI * i / 64);
im_coef[i] = (float)sin(2 * CFFT_PI * i / 64);
} } void ifft_coef_gen_256(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_256; i++)
{
re_coef[i] = (float)cos(2 * CFFT_PI * i / 256);
im_coef[i] = (float)sin(2 * CFFT_PI * i / 256);
} } void fft_coef_gen_256(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_256; i++)
{
re_coef[i] = (float)cos(-2 * CFFT_PI * i / 256);
im_coef[i] = (float)sin(-2 * CFFT_PI * i / 256);
} } void ifft_256(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_256);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_256);
re_tmp = (float*)malloc(sizeof(float) * 256);
im_tmp = (float*)malloc(sizeof(float) * 256);
ifft_coef_gen_256(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
256
);
for (int i = 0; i < 256; i += 64)
{
if (i == 192 && end_flag == 1)
{
start_flag = 1;
}
ifft_64(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
ifft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
64
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} } void fft_256(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_256);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_256);
re_tmp = (float*)malloc(sizeof(float) * 256);
im_tmp = (float*)malloc(sizeof(float) * 256);
fft_coef_gen_256(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
256
);
for (int i = 0; i < 256; i += 64)
{
if (i == 192 && end_flag == 1)
{
start_flag = 1;
}
fft_64(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
fft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
64
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} } void ifft_coef_gen_1024(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_1024; i++)
{
re_coef[i] = (float)cos(2 * CFFT_PI * i / 1024);
im_coef[i] = (float)sin(2 * CFFT_PI * i / 1024);
} } void fft_coef_gen_1024(
float* re_coef,
float* im_coef ) {
int i;
for (i = 0; i < COEF_NUM_1024; i++)
{
re_coef[i] = (float)cos(-2 * CFFT_PI * i / 1024);
im_coef[i] = (float)sin(-2 * CFFT_PI * i / 1024);
} } void ifft_1024(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_1024);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_1024);
re_tmp = (float*)malloc(sizeof(float) * 1024);
im_tmp = (float*)malloc(sizeof(float) * 1024);
ifft_coef_gen_1024(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
1024
);
for (int i = 0; i < 1024; i += 256)
{
if (i == 768 && end_flag == 1)
{
start_flag = 1;
}
ifft_256(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
ifft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
256
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} } void fft_1024(
float* re_in,
float* im_in,
float* re_out,
float* im_out,
int end_flag ) {
static int start_flag = 1;
static float* re_coef, * im_coef;
static float* re_tmp, * im_tmp;
if (start_flag)
{
re_coef = (float*)malloc(sizeof(float) * COEF_NUM_1024);
im_coef = (float*)malloc(sizeof(float) * COEF_NUM_1024);
re_tmp = (float*)malloc(sizeof(float) * 1024);
im_tmp = (float*)malloc(sizeof(float) * 1024);
fft_coef_gen_1024(
re_coef,
im_coef
);
start_flag = 0;
}
//transformation proccess
addr_sort(
re_in,
im_in,
re_tmp,
im_tmp,
1024
);
for (int i = 0; i < 1024; i += 256)
{
if (i == 768 && end_flag == 1)
{
start_flag = 1;
}
fft_256(
re_tmp + i,
im_tmp + i,
re_out + i,
im_out + i,
start_flag
);
}
fft_layer_transform(
re_out,
im_out,
re_coef,
im_coef,
256
);
//clean up
if (end_flag)
{
free(re_coef);
free(im_coef);
free(re_tmp);
free(im_tmp);
} }
#endif
333
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
