体验入口地址https://qianwen.aliyun.com/chat
问题:帮我用C语言写一个FFT。
通义千问:以下是一个基于C语言的FFT实现,用于计算一个长度为N的复数序列的FFT。
- #include <stdio.h>
- #include <complex.h>
- #define PI 3.14159265358979323846
- // 定义复数结构体
- typedef struct {
- double real;
- double imaginary;
- } Complex;
- // 计算两个复数的加法
- Complex operator+(Complex a, Complex b) {
- Complex c;
- c.real = a.real + b.real;
- c.imaginary = a.imaginary + b.imaginary;
- return c;
- }
- // 计算两个复数的减法
- Complex operator-(Complex a, Complex b) {
- Complex c;
- c.real = a.real - b.real;
- c.imaginary = a.imaginary - b.imaginary;
- return c;
- }
- // 计算两个复数的乘法
- Complex operator*(Complex a, Complex b) {
- Complex c;
- c.real = a.real * b.real - a.imaginary * b.imaginary;
- c.imaginary = a.real * b.imaginary + a.imaginary * b.real;
- return c;
- }
- // 计算复数的模
- double abs(Complex a) {
- return sqrt(a.real * a.real + a.imaginary * a.imaginary);
- }
- // 计算复数的辐角
- double angle(Complex a) {
- return atan2(a.imaginary, a.real);
- }
- // FFT算法
- void fft(Complex *a, int n, int inverse) {
- int i, j, k;
- Complex w, w0, w1, w2, w3;
- // Bit-reverse permutation
- for (i = 1, j = 0; i < n; ++i) {
- if (i > j) {
- Complex temp = a[i];
- a[i] = a[j];
- a[j] = temp;
- }
- k = n >> 1;
- while (k >= 2 && j >= k) {
- j -= k;
- k >>= 1;
- }
- j += k;
- }
- // Power-of-2 FFT
- for (int len = 2; len <= n; len <<= 1) {
- double angle = 2 * PI / len * inverse;
- Complex wlen(cos(angle), sin(angle));
- for (i = 0; i < n; i += len) {
- w = 1;
- for (j = 0; j < len / 2; ++j) {
- w1 = w * a[i + j + len / 2];
- w2 = w * a[i + j];
- a[i + j] = w1 + w2;
- a[i + j + len / 2] = w1 - w2;
- w *= wlen;
- }
- }
- }
- // Bit-reverse permutation
- for (i = 1, j = 0; i < n; ++i) {
- if (i > j) {
- Complex temp = a[i];
- a[i] = a[j];
- a[j] = temp;
- }
- k = n >> 1;
- while (k >= 2 && j >= k) {
- j -= k;
- k >>= 1;
- }
- j += k;
- }
- }
- // 计算复数序列的IFFT
- void ifft(Complex *a, int n) {
- fft(a, n, 1);
- for (int i = 0; i < n; ++i) {
- a[i].real /= n;
- a[i].imaginary /= n;
- }
- }
- int main() {
- int n = 8;
- Complex a[n] = { {0, 1}, {2, 3}, {4, 5}, {6, 7}, {8, 9}, {10, 11}, {12, 13}, {14, 15} };
- fft(a, n, 0);
- for (int i = 0; i < n; ++i) {
- printf("(%f + %fi) ", a[i].real, a[i].imaginary);
- }
- printf("
- ");
- return 0;
- }
在上面的代码中,我们首先定义了一个复数结构体,包含实部和虚部两个成员。然后实现了加法、减法和乘法三个复数运算函数。接下来,我们实现了计算复数模和辐角的函数。接着,我们实现了FFT算法,包括位反转、长度为2的FFT和长度为N的FFT三个部分。最后,我们实现了IFFT算法,并在主函数中演示了如何使用FFT算法计算一个复数序列的FFT。
|