一个基于C++的FFT实现方法—librow
admin 于 2017年09月21日 发表在 数字信号处理

前几天看到一个不错的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)输出相关的图形,如下表:

注意:本站所有文章除特别说明外,均为原创,转载请务必以超链接方式并注明作者出处。 标签:FFT,信号处理,处理算法