WhyCan Forum(哇酷开发者社区)

bunny

会员

注册时间： 2020-05-23

已发帖子： 154

积分： 154

fft又来了，这次更简单，代码才30行

static double buffer[N]; // N = 2^(k+1)
void fft(const double* src_x, const double* src_y, double* dst_x, double* dst_y, int k) {
    int i, n;
    double temp, temp_x, temp_y;
    if (0 == k) {
        *dst_x = *src_x;
        *dst_y = *src_y;
        return;
    }
    n = 1 << (k - 1);
    for (i = 0; i < n; i++) {
        buffer[i] = src_x[i * 2];
        buffer[i + n] = src_y[i * 2];
        buffer[i + 2 * n] = src_x[i * 2 + 1];
        buffer[i + 3 * n] = src_y[i * 2 + 1];
    }
    fft((const double*)buffer + 2 * n, (const double*)buffer + 3 * n, dst_x, dst_y, k - 1);
    fft((const double*)buffer, (const double*)buffer + n, (double*)buffer + 2 * n, (double*)buffer + 3 * n, k - 1);
    for (i = 0; i < n; i++) {
        temp = i * M_PI / n;
        temp_x = cos(temp);
        temp_y = sin(temp);
        buffer[i] = dst_y[i] * temp_y + dst_x[i] * temp_x;
        buffer[i + n] = dst_y[i] * temp_x - dst_x[i] * temp_y;
        dst_x[i] = buffer[i + 2 * n] + buffer[i];
        dst_y[i] = buffer[i + 3 * n] + buffer[i + n];
        dst_x[i + n] = buffer[i + 2 * n] - buffer[i];
        dst_y[i + n] = buffer[i + 3 * n] - buffer[i + n];
    }
}

函数使用说明：
void fft(const double* src_x, const double* src_y, double* dst_x, double* dst_y, int k)
src_x：输入数据实部
src_y：输入数据虚部
dst_x：输出数据实部
dst_y：输出数据虚部
k：表征fft点数为2^k，比如1024点的fft，k=10；4096点的fft，k=12；

bunny

会员

注册时间： 2020-05-23

已发帖子： 154

积分： 154

Re: fft又来了，这次更简单，代码才30行

牛叉, 递归调用, 路过学习

不优雅的代码也不敢往晕哥的坑网里塞啊

bunny

会员

注册时间： 2020-05-23

已发帖子： 154

积分： 154

Re: fft又来了，这次更简单，代码才30行

请问有大佬测试过这个代码吗，我自己写的的计算出来全是0

int main()
{
	double x[16],y[16];
	int i=0;
	for(i=0;i<16;i++)
	{
		y[i]=0;
		x[i]=0;
	}
	x[0]=1.0;
	x[1]=1.0;
	double dx[16]={0},dy[16]={0};

	fft(x,y,dx,dy,4);
	for (i = 0; i < 16; i++) {
        printf((const char*)"dat[%d] = %f + %f * i\n", i, dx[i], dy[i]);
    }
	return 0;
}

你的测试是没问题的，是我代码的疏漏，出现了内存区域的覆盖，第一版本中使用malloc free，没有这个问题，后来改成静态内存的时候逻辑没理顺，出现的内存覆盖的BUG，由于未经严格测试，没有发现这个BUG，经你提醒才抓出来
上传最早的版本（引用了malloc free函数）

void fft(const double* src_x, const double* src_y, double* dst_x, double* dst_y, int k) {
    int i, n;
    double temp, temp_x, temp_y;
    double* buffer;
    if (0 == k) {
        *dst_x = *src_x;
        *dst_y = *src_y;
        return;
    }
    n = 1 << (k - 1);
    buffer = (double*)malloc(4 * n * sizeof(double));
    for (i = 0; i < n; i++) {
        buffer[i] = src_x[i * 2];
        buffer[i + n] = src_y[i * 2];
        buffer[i + 2 * n] = src_x[i * 2 + 1];
        buffer[i + 3 * n] = src_y[i * 2 + 1];
    }
    fft((const double*)buffer + 2 * n, (const double*)buffer + 3 * n, dst_x, dst_y, k - 1);
    fft((const double*)buffer, (const double*)buffer + n, (double*)buffer + 2 * n, (double*)buffer + 3 * n, k - 1);
    for (i = 0; i < n; i++) {
        temp = i * M_PI / n;
        temp_x = cos(temp);
        temp_y = sin(temp);
        buffer[i] = dst_y[i] * temp_y + dst_x[i] * temp_x;
        buffer[i + n] = dst_y[i] * temp_x - dst_x[i] * temp_y;
        dst_x[i] = buffer[i + 2 * n] + buffer[i];
        dst_y[i] = buffer[i + 3 * n] + buffer[i + n];
        dst_x[i + n] = buffer[i + 2 * n] - buffer[i];
        dst_y[i + n] = buffer[i + 3 * n] - buffer[i + n];
    }
    free(buffer);
}

bunny

会员

注册时间： 2020-05-23

已发帖子： 154

积分： 154

Re: fft又来了，这次更简单，代码才30行

更新的fft函数测试截图

bunny

会员

注册时间： 2020-05-23

已发帖子： 154

积分： 154

Re: fft又来了，这次更简单，代码才30行

可以了，多谢大佬:D

也谢谢你帮我抓BUG