利用FFT计算频谱图

给出了完整的FFT的C语言代码，以及MATLAB的仿真效果，主要解决了FFT计算频谱图的问题。注意大部分周期信号的最小非零频率分量就是基波，但是还有很多信号并不是这样：比如y=sin(2πx) +sin(3πx)的信号周期为2，对应的基波角频率为π，而频谱图中角频率为π的那个点的幅度为0。

#include “fft.h“

inline void swap （float &a float &b）
{
    float t;
    t = a;
    a = b;
    b = t;
}

void bitrp （float xreal [] float ximag [] int n）
{
    // 位反转置换 Bit-reversal Permutation
    int i j a b p;

    for （i = 1 p = 0; i < n; i *= 2）
        {
        p ++;
        }
    for （i = 0; i < n; i ++）
        {
        a = i;
        b = 0;
        for （j = 0; j < p; j ++）
            {
            b = （b << 1） + （a & 1）;    // b = b * 2 + a % 2;
            a >>= 1;        // a = a / 2;
            }
        if （ b > i）
            {
            swap （xreal [i] xreal [b]）;
            swap （ximag [i] ximag [b]）;
            }
        }
}

void FFT（float xreal [] float ximag [] int n）
{
    // 快速傅立叶变换，将复数 x 变换后仍保存在 x 中，xreal ximag 分别是 x 的实部和虚部
    float wreal [N / 2] wimag [N / 2] treal timag ureal uimag arg;
    int m k j t index1 index2;

    bitrp （xreal ximag n）;

    // 计算 1 的前 n / 2 个 n 次方根的共轭复数 W‘j = wreal [j] + i * wimag [j]  j = 0 1 ...  n / 2 - 1
    arg = - 2 * PI / n;
    treal = cos （arg）;
    timag = sin （arg）;
    wreal [0] = 1.0;
    wimag [0] = 0.0;
    for （j = 1; j < n / 2; j ++）
        {
        wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
        wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
        }

    for （m = 2; m <= n; m *= 2）
        {
        for （k = 0; k < n; k += m）
            {
            for （j = 0; j < m / 2; j ++）
                {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                treal = wreal [t] * xreal [index2] - wimag [t] * ximag [index2];
                timag = wreal [t] * ximag [index2] + wimag [t] * xreal [index2];
                ureal = xreal [index1];
                uimag = ximag [index1];
                xreal [index1] = ureal + treal;
                ximag [index1] = uimag + timag;
                xreal [index2] = ureal - treal;
                ximag [index2] = uimag - timag;
                }
            }
        }
}

void  IFFT （float xreal [] float ximag [] int n）
{
    // 快速傅立叶逆变换
    float wreal [N / 2] wimag [N / 2] treal timag ureal uimag arg;
    int m k j t index1 index2;

    bitrp （xreal ximag n）;

    // 计算 1 的前 n / 2 个 n 次方根 Wj = wreal [j] + i * wimag [j]  j = 0 1 ...  n / 2 - 1
    arg = 2 * PI / n;
    treal = cos （arg）;
    timag = sin （arg）;
    wreal [0] = 1.0;
    wimag [0] = 0.0;
    for （j = 1; j < n / 2; j ++）
        {
        wreal [j] = wreal [j - 1] * treal - wimag [j - 1] * timag;
        wimag [j] = wreal [j - 1] * timag + wimag [j - 1] * treal;
        }

    for （m = 2; m <= n; m *= 2）
        {
        for （k = 0; k < n; k += m）
            {
            for （j = 0; j < m / 2; j ++）
                {
                index1 = k + j;
                index2 = index1 + m / 2;
                t = n * j / m;    // 旋转因子 w 的实部在 wreal [] 中的下标为 t
                treal = wreal [t] *

 属性            大小     日期    时间   名称
----------- ---------  ---------- -----  ----

     文件        739  2018-07-21 16:28  FFT计算频谱图\C语言实现\Debug\.qmake.stash

     文件      90994  2018-08-20 23:41  FFT计算频谱图\C语言实现\Debug\debug\FFT.exe

     文件      25926  2018-08-20 23:41  FFT计算频谱图\C语言实现\Debug\debug\fft.o

     文件      17165  2018-08-20 23:41  FFT计算频谱图\C语言实现\Debug\debug\main.o

     文件      50744  2018-07-21 19:13  FFT计算频谱图\C语言实现\Debug\FFT

     文件      46808  2018-07-21 19:13  FFT计算频谱图\C语言实现\Debug\fft.o

     文件       7164  2018-07-21 18:53  FFT计算频谱图\C语言实现\Debug\input.txt

     文件      38043  2018-07-21 16:47  FFT计算频谱图\C语言实现\Debug\input1.txt

     文件        111  2018-07-21 16:34  FFT计算频谱图\C语言实现\Debug\input11.txt

     文件      32728  2018-07-21 19:09  FFT计算频谱图\C语言实现\Debug\main.o

     文件      26586  2018-08-20 22:54  FFT计算频谱图\C语言实现\Debug\Makefile

     文件      14483  2018-08-20 22:54  FFT计算频谱图\C语言实现\Debug\Makefile.Debug

     文件      14522  2018-08-20 22:54  FFT计算频谱图\C语言实现\Debug\Makefile.Release

     文件     147558  2018-07-21 18:59  FFT计算频谱图\C语言实现\Debug\output.txt

     文件       6442  2018-08-21 00:53  FFT计算频谱图\C语言实现\fft.cpp

     文件        662  2018-08-21 00:53  FFT计算频谱图\C语言实现\fft.h

     文件        142  2018-07-21 16:24  FFT计算频谱图\C语言实现\FFT.pro

     文件      23607  2018-08-21 00:21  FFT计算频谱图\C语言实现\FFT.pro.user

     文件      23662  2018-07-21 16:23  FFT计算频谱图\C语言实现\FFT.pro.user.6ed20e1

     文件         92  2018-08-20 22:57  FFT计算频谱图\C语言实现\main.cpp

     文件       1164  2018-07-21 16:46  FFT计算频谱图\matlab验证试验\FFT.m

     文件        261  2018-07-21 18:43  FFT计算频谱图\readme.txt

     文件     347069  2018-07-21 19:13  FFT计算频谱图\实验结果\2018-07-21 19-13-22屏幕截图.png

     文件     363558  2018-07-21 19:14  FFT计算频谱图\实验结果\2018-07-21 19-14-07屏幕截图.png

    ..A.SH.       162  2018-08-20 22:47  FFT计算频谱图\实验结果\~$试验结果记录.docx

     文件     792980  2018-07-21 19:25  FFT计算频谱图\实验结果\试验结果记录.docx

     目录          0  2018-08-21 00:45  FFT计算频谱图\C语言实现\Debug\debug

     目录          0  2018-08-20 22:54  FFT计算频谱图\C语言实现\Debug\release

     目录          0  2018-08-21 00:45  FFT计算频谱图\C语言实现\Debug

     目录          0  2018-08-21 00:45  FFT计算频谱图\C语言实现

............此处省略6个文件信息