7/1-4 快速傅里叶变换FFT的matlab程序以及C++程序

7/1-4 快速傅里叶变换FFT的matlab程序以及C++程序

1.matlab代码

\qquad 用open命令打开fft内置函数后给出的介绍如下：

%FFT Discrete Fourier transform.
%   FFT(X) is the discrete Fourier transform (DFT) of vector X.  For
%   matrices, the FFT operation is applied to each column. For N-D
%   arrays, the FFT operation operates on the first non-singleton
%   dimension.
%
%   FFT(X,N) is the N-point FFT, padded with zeros if X has less
%   than N points and truncated if it has more.
%
%   FFT(X,[],DIM) or FFT(X,N,DIM) applies the FFT operation across the
%   dimension DIM.
%   
%   For length N input vector x, the DFT is a length N vector X,
%   with elements
%                    N
%      X(k) =       sum  x(n)*exp(-j*2*pi*(k-1)*(n-1)/N), 1 <= k <= N.
%                   n=1
%   The inverse DFT (computed by IFFT) is given by
%                    N
%      x(n) = (1/N) sum  X(k)*exp( j*2*pi*(k-1)*(n-1)/N), 1 <= n <= N.
%                   k=1
%
%   See also FFT2, FFTN, FFTSHIFT, FFTW, IFFT, IFFT2, IFFTN.

%   Copyright 1984-2005 The MathWorks, Inc.

%   Built-in function.

现在来对比一下离散傅里叶变换与快速傅里叶变换所需时间：

t=linspace(1,10,100);
 x=sin(t);
 t1=clock;
 y1=fft(x);
 t2=clock;
 etime(t2,t1)

ans =

    0.0090

\qquad 由上述代码得知matlab自带函数fft计算100个信号$(\sin (0:0.01:1))$所需的时间是$0.009$秒。

 t=linspace(1,10,100);
 N=length(t);
 x=sin(t);
 t1=clock;
%  y1=fft(x);
 k=0;
 y2=zeros(1,N);
 for i=1:N
     for j=1:N
         k=k+exp(-i*2*pi*i*j/N)*x(j);
         y2(i)=k;
     end
 end
 t2=clock;
 etime(t2,t1)

ans =

    0.0320

\qquad 由上述代码得知DFT计算100个信号$(\sin (0:0.01:1))$所需的时间是$0.0320$秒。若将计时函数放在第7行，则可省去申请内存的时间，所需时间大约为$0.0150$秒，仍然是比fft慢很多，这仅仅是100个信号，若处理成千上万个信号则用DFT则会慢很多。
\qquad 下面给出FFT的代码：

function f=fftt(x)
N=length(x);
f=zeros(1,N);
for k=1:N/2
    for n=1:N/2
        f(2*k-1)=f(2*k-1)+exp(-1i*2*pi*(2*k-2)*(n-1)/N)*(x(n)+x(n+N/2));
        f(2*k)=f(2*k)+exp(-1i*2*pi*(2*k-1)*(n-1)/N)*(x(n)-x(n+N/2));
    end
end
end

下面对比DFT与FFT计算时间

 t=linspace(0,1,100);
 N=length(t);
 x=sin(t);
 y1=fft(x);
 y2=zeros(1,N);
 t1=clock;
 for k=1:N
     for n=1:N
         y2(k)=y2(k)+exp(-1i*2*pi*(k-1)*(n-1)/N)*x(n);
     end
 end
 t2=clock;
 time1=etime(t2,t1)
 t3=clock;
 y3=fftt(x);
 t4=clock;
 time2=etime(t4,t3)

\qquad 只取100个点会出现如下情况，即用FFT算法不消耗时间，下面改为1000个点，时间如下：

time1 =

    0.2620


time2 =

    0.0940

2.C++代码

#include <iostream>
#include <math.h>
#include <string.h>
using namespace std;
#define pi acos(-1)

int main()
{
    float a[100];
    int N=sizeof(a)/sizeof(a[0]);
    float y[N];
    float f[N];
    for(int l=0;l<N;l++)
    {
        a[l]=sin((float)l/100);
        y[l]=0;
        f[l]=0;
    }
    for(int k=0;k<N/2;k++)
    {
        for(int n=0;n<N/2;n++)
        {
            y[2*k]=y[2*k]+cos(4*k*n*pi/N)*(a[n]+a[n+N/2]);
            y[2*k+1]=y[2*k+1]+(a[n]-a[n+N/2])*cos(2*(2*k+1)*n*pi/N);
        }
    }
    return 0;
}

\qquad 为了省事，我直接将计算中的$e^{i\omega \pi }$定义为$\cos (i\omega \pi )$,舍掉了虚部部分。计算结果与MATLAB相同（实部）。