fft算法总结

前言

电赛期间总结的FFT算法

一、fft是什么？

FFT（Fast Fourier Transformation）是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。 

二、使用步骤

电赛时总结的fft算法，经过实际测试，可用！

#include <assert.h>
#include <math.h>
#include <stdint.h>
#include <stdio.h>

#define LN 10
#define N (1 << LN)

#define M_PI 3.14159265358979323846

void fft_origin_sort(float *in_real, float *in_imag, float *out_real,
                     float *out_imag, uint32_t ln) {
  uint32_t LEN = 1 << ln;
  for (int i = 0; i < LEN; i++) {
    uint32_t _idx = 0;
    for (int j = 0; j < ln; j++) {
      _idx = (_idx << 1) | ((i >> j) & 0x01);
    }
    if (in_real != NULL && out_real != NULL) {
      out_real[i] = in_real[_idx];
    } else if (out_real != NULL) {
      out_real[i] = 0;
    }
    if (in_imag != NULL && out_imag != NULL) {
      out_imag[i] = in_imag[_idx];
    } else if (out_imag != NULL) {
      out_imag[i] = 0;
    }
  }
}

void fft_origin_Wn(float *Wn_real, float *Wn_imag, uint32_t ln) {
  uint32_t LEN = 1 << ln;
  for (int i = 0; i < LEN / 2; ++i) {
    float k = (float)i / (float)LEN;
    float a = k * (2.0 * M_PI);
    Wn_real[i] = cosf(a);
    Wn_imag[i] = -sinf(a);
  }
}

static inline void cpx_mul(float real1, float imag1, float real2, float imag2,
                           float *out_real, float *out_imag) {
  *out_real = real1 * real2 - imag1 * imag2;
  *out_imag = real1 * imag2 + real2 * imag1;
}

void fft_origin_cal(float *real, float *imag, float *Wn_real, float *Wn_imag,
                    uint32_t ln) {
  for (int i = 0; i < ln; ++i) {
    uint32_t m = 1 << (ln - 1 - i);
    uint32_t n = 1 << i;
    for (int j = 0; j < m; ++j) {
      for (int k = 0; k < n; ++k) {
        uint32_t _idx = j * n * 2 + k;
        uint32_t X2_idx = _idx + n, Wn_idx = k * m;
        float X2Wn_r, X2Wn_i;
        cpx_mul(real[X2_idx], imag[X2_idx], Wn_real[Wn_idx], Wn_imag[Wn_idx],
                &X2Wn_r, &X2Wn_i);
        real[X2_idx] = real[_idx] - X2Wn_r;
        imag[X2_idx] = imag[_idx] - X2Wn_i;
        real[_idx] = real[_idx] + X2Wn_r;
        imag[_idx] = imag[_idx] + X2Wn_i;
      }
    }
  }
}

float vmod(float r, float i) { return sqrtf(r * r + i * i); }

void base_freq(float *real, float *imag, uint32_t ln) {
  const float epsilon = 1e-2;
  float _diff = 0;
  uint32_t LEN = (1 << ln) / 2;
  float A[LEN];
  for (int i = 0; i < LEN; ++i) {
    A[i] = vmod(real[i], imag[i]) / LEN;
  }
  for (int i = 1; i < LEN; ++i) {
    float diff = A[i] - A[i - 1];
    if (_diff > epsilon && diff < -epsilon) {
      printf("%d,%f\n", i - 1, A[i - 1]);
      return;
    }
    _diff = diff;
  }
}

void genx(float *f, float k, uint32_t n) {
  for (int i = 0; i < n; ++i) {
    f[i] = sin(i * (M_PI * 2) / n * k) * 600;
  }
}

int main() {
  float Wn_r[N / 2], Wn_i[N / 2];
  float _r[N] = {0};
  float r[N], i[N];
  fft_origin_Wn(Wn_r, Wn_i, LN);
  for (int w = 1; w < N / 2 - 1; ++w) {
    genx(_r, w, N);
    fft_origin_sort(_r, NULL, r, i, LN);
    fft_origin_cal(r, i, Wn_r, Wn_i, LN);
    base_freq(r, i, LN);
  }
}

---

总结

电赛中的仪器仪表类或者信号处理类题目经常要进行对信号的处理，其中最常见的就是信号的频谱分析，频谱分析最常见的就是FFT，该方法简单易懂，适合电赛时进行学习用来做为信号处理的手段。