一维FFT变换

140 阅读 0 评论 93 点赞

我是靠谱客的博主笨笨白开水，最近开发中收集的这篇文章主要介绍一维FFT变换，觉得挺不错的，现在分享给大家，希望可以做个参考。

概述

#include "stdafx.h"
#include "D:/CImg/CImg.h"
#include <iostream>

using namespace std;
using namespace cimg_library;

// The lines below are necessary when using a
// non-standard compiler such as visualcpp6.
#ifdef cimg_use_visualcpp6
#define std
#endif
#ifdef min
#undef min
#undef max
#endif

// Define circumference ratio constant
#define PI 3.14159265358

int main(int argc, char* argv[])
{

// Frequency
const int feq = 25;
// Sampling rate
const int sample_rate = 1024;
// Ampitude
const float ampitude = 50.0f;
// Colors
const unsigned char red[] = { 255, 0, 0 }, yellow[] = { 255, 255, 0 },
green[] = { 0, 255, 0 };

// 1-D signal
CImg<float> signal(1, 1024, 1, 1, 0);
// Filling signal with values
cimg_forXY( signal, x, y )
{
signal( x, y ) = ampitude*cos( 2*PI*feq*( y/ 1023.0 ) ) + 125;
}

// Create a black image as drawing panel
CImg<unsigned char> signal_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
signal_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 200, -200, green, 1.0f );
// Draw the signal with cubic spline
signal_visu.draw_graph( signal, red, 1, 2, 0, 255, 0 );
// Display our drawing result
CImgDisplay main_disp( signal_visu, "signal" );

// FFT result container
CImgList<float> fft = signal.get_FFT( 'y' );
// Create a drawing panel
CImg<unsigned char> fft_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
fft_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 300, 0, green, 1.0f );;
// Power of FFT. I scale the power in order to display it properly
CImg<float> power = ( fft(0).get_pow(2.0) + fft(1).get_pow(2.0) ).sqrt() * 0.05;
// Draw the signal with cubic spline
fft_visu.draw_graph( power , red, 1, 2, 0, 255, 0 );

CImgDisplay fft_disp( fft_visu, "FFT: Magitude" );

CImg<float> phase(1, 1024, 1, 1, 0.0f);
cimg_forXY( phase, x, y )
{
phase(x, y) = atan(abs((fft(0))(x, y) / (fft(1))(x, y)));
}
CImg<unsigned char> phase_visu( 800, 300, 1, 3, 0 );
// Draw grid and axis in 'drawing panel'
phase_visu.draw_grid( 50, 50, 0, 0, 0, 0, yellow, 0.5 ).draw_axis( 0, 1000, 300, 0, green, 1.0f );;
// Draw the signal with cubic spline
phase_visu.draw_graph( phase , red, 1, 2, 0, 255, 0 );

CImgDisplay phase_disp( phase_visu, "FFT: Phase");

while ( !main_disp.is_closed )
{
main_disp.wait();
}
return 0;
}

今天给大家介绍的例子是利用CImg进行快速傅立叶变换。傅立叶变换被广泛的应用于信号处理行业，为信号分析提供了一种强有力的分析方法。老傅在他的一篇文章中，提出周期性的函数都可以利用正弦函数来表示。当时，审稿人中有两位当时大名顶顶的数学家：拉格朗日与拉普拉斯，这两个人对于此文章是否可以发表持相反意见。拉格朗日支持方认为该文章不能发表，他指出该方法不能表示带有边角(corner)的函数，而拉普拉斯和其他学者则支持发表。最后，迫于拉格朗日的权威（牛人就是牛人），这篇文章在他西去之后的15年才得以发表。不过这段时间，傅立叶也没有闲着，也就没有太在意这件事情。是的，人有的时候不能太在意一些东西，随遇而安最好。到底上面两位牛人谁的观点对呢？这要一分为二看：站在老拉的角度来看，这个方法是有‘问题’，但我们可以无限逼近元函数，使得它们之间的能量几乎为零，这样一来，傅立叶变换是可行的。此外，'对于离散信号来说，傅立叶信号是准确的，就像7 = 3 + 4一样'.(1)

简单的说，傅立叶变换其实就是统计某一频率的信号在待分析信号中出现的频率（次数），貌似和概率分布差不多(个人观点)。傅立叶变换可以通过多种方法来解决，例如解线性方程组，求得每一幅值大小。这种方法简单直接，但是耗时巨大。我们平时所要分析的数据量往往很大。第二种方法是通过卷积来求。这种方法直接向我们表明了傅立叶变换的实质。让一个频率为f正弦基信号沿着输入信号进行卷积，检测与其频率相同的子信号。我们知道，在一个周期内，正弦曲线有个性质：同频率的正弦函数相乘积分为1，不同频率的则为0。(上面是个不太准确的定义，严格定义请参考高等数学）。这样，‘游走’一圈之后，对应频率位值的值会非常大，而其他地方则为0.（这是理想的结果，现实中的信号多带有噪声，所以结果有些出入）。第三种方法就是顶顶大名的快速傅立叶变换了（FFT）。这个方法可以把计算时间复杂度维持在nlog(n)数量级上，这在数据量非常大的时候，速度能够提高成千倍。FFT的思想早在高斯的时候就诞生了，不过直到1967年的时候才由两位牛人利用计算机实现了。如果当时有计算机的话，估计就轮不到这两哥们了。关于FFT的计算方法，我将在陆续给出详细解释（说实话，我到现在还没有完全弄懂，所以就不误导大家了）。

OK，废话说了太多了，现在来分析一下这个小程序吧。这个程序是利用CImg图像库做的。CImg是一个不错的图像库（当然，它也可以用来处理1D信号了），该库非常的轻巧（就一个头文件，使用的时候include进去就可以了），该库由于利用的模板编程的方法，所以使得该程序具有良好的可重用性。程序的前几行，主要是定义了一些常量，紧接着就是初始化信号，这个信号具50Hz的频率和50.0的幅值，信号的采样频率为1024（这个为了FFT， 2的次方）。接下来，就是把我们的输入信号画出来。然后，利用get_FFT函数进行FFT（参数’y‘表明是在y方向进行1D变换）。最后，求取FFT的Magitude。

在结果图像中，0的位置为直流分量（工程师都这么叫），其实就是输入信号的均值。我们发现在50的位置，值特别大（这里我加了缩放因子的，否则显示出来不太美观）。眼力好的人或许发现这个图像貌似的是对称。恩，是的。所以在显示的时候，我们可以只显示一半。从图中，我们可以知道，输入信号有频率为50的正弦信号组成。实际信号中，可能会有几个峰值，最大的叫first frequency,其次呢？当然是secocnd frequency喽，呵呵。

Magitude用来表示能量，那相位（Phase）呢？它是用来表示正弦信号的初始位置的，比如cos(2×PI×f + 45），那么这个信号的初始相位就为45（这里的单位是度，而函数库中所提够的算法结果一般为弧度）。由于我们的信号初始相位为0，所以你看到第三幅图像中所有的值均为0.