YY分享

前几天看到一个不错的FFT变换类，基于C++语言的，感觉不错，记录在此，万一以后用到也好找。

1. 官网下载librow，并解压：

complex.h 和 complex.cpp 中定义了各种复数运算的基本单元以及基本运算的具体实现。

fft.h 和 fft.c 中定义了相关运算函数，包含傅里叶反变换、FFT具体实现等。

注：FFT相关原理，请查看博文《数字信号处理算法(2)-快速傅里叶变换(FFT)》，此处不贴。

2. 调用接口函数实现FFT，大体需要如下几步：

（1）初始化参数

double x[N]={0.0};
complex *pSignal = new complex[N];
FILE *fp;      //创建文本操作对象

（2）生成原始时域信号：

for(i=0; i<N; i++)
{
  x[i]=sin(2*PI*(5.0*i/(N-1)))+sin(2*PI*(5.0*3.0*i/(N-1)))/3.0+sin(2*PI*(5.0*5.0*i/(N-1)))/5.0;
}

（3）将生成的时域值传递给定义的complex

for(i=0; i<N; i++)
{
  pSignal[i] = x[i];
}

（4）调用librow中对应接口函数

CFFT::Forward(pSignal,N);

（5）保存并打印输出（可选）

for(i=0; i<N; i++)
{
  pSignal[i] = sqrt(pSignal[i].re()*pSignal[i].re()+pSignal[i].im()*pSignal[i].im());
  fprintf(fp,"%lf\t%lf\n",x[i],pSignal[i]);    //存储数值到文本
}

（6）释放内存占用

fclose(fp);   //关闭文本
delete[] pSignal;

3. 整个main.cpp函数如下：

#include <iostream>
#include "fft.h"
#include <math.h>
#include <stdio.h>

//定义PI
#define PI 3.141593
#define N  64

int main()
{
    int i=0;
    double x[N]={0.0};
    complex *pSignal = new complex[N];
    FILE *fp;   //存储最终计算值
    fp=fopen("data.dat","w");   //打开文本

    //1.生成原始信号x[]
    for(i=0; i<N; i++) //产生正弦三次谐波叠加而成的方波
    {
        x[i]=sin(2*PI*(5.0*i/(N-1)))+sin(2*PI*(5.0*3.0*i/(N-1)))/3.0+sin(2*PI*(5.0*5.0*i/(N-1)))/5.0;
    }

    //2.将原始值赋给pSignal[]
    for(i=0; i<N; i++)
    {
        pSignal[i] = x[i];
    }

    //3.计算FFT数据
    CFFT::Forward(pSignal,N);

    //4.输出计算结果
    for(i=0; i<N; i++)
    {
        pSignal[i] = sqrt(pSignal[i].re()*pSignal[i].re()+pSignal[i].im()*pSignal[i].im());
        fprintf(fp,"%lf\t%lf\n",x[i],pSignal[i]);    //存储数值
    }

    //5.释放内存占用
    fclose(fp);    //关闭文本
    delete[] pSignal;
}

4. 通过Qtiplot绘制对应数据的图形，如下：

补充： 通过python自带的FFT函数验证上述算法，如下：

（1）Python代码：

# -*- coding: utf-8 -*-
"""
Created on Tue Sep 19 23:03:53 2017

@author: lu
"""

import matplotlib.pyplot as plt
import numpy as np

N=64
x = np.linspace(0,N-1,N)
y = np.sin(2*np.pi*(5*x/(N-1)))+ np.sin(2*np.pi*(5*3*x/(N-1)))/3+ np.sin(2*np.pi*(5*5*x/(N-1)))/5

plt.figure(figsize=(9,5))
plt.plot(x,y,label="$f(x)$",color="red",linewidth=2)

plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Origin")
plt.ylim(-1.1,1.1)
plt.legend()
plt.show()

y1 = np.fft.fft(y);
y2 = np.sqrt(y1.real**2 + y1.imag**2)
plt.figure(figsize=(9,5))
plt.plot(x,y2,label="$F(s)$",color="red",linewidth=2)
plt.xlabel("sample")
plt.ylabel("F(s)")
plt.title("FFT")
plt.ylim(0,35)
plt.legend()
plt.show()

（2）输出相关的图形，如下表：