视频1 视频21 视频41 视频61 视频文章1 视频文章21 视频文章41 视频文章61 推荐1 推荐3 推荐5 推荐7 推荐9 推荐11 推荐13 推荐15 推荐17 推荐19 推荐21 推荐23 推荐25 推荐27 推荐29 推荐31 推荐33 推荐35 推荐37 推荐39 推荐41 推荐43 推荐45 推荐47 推荐49 关键词1 关键词101 关键词201 关键词301 关键词401 关键词501 关键词601 关键词701 关键词801 关键词901 关键词1001 关键词1101 关键词1201 关键词1301 关键词1401 关键词1501 关键词1601 关键词1701 关键词1801 关键词1901 视频扩展1 视频扩展6 视频扩展11 视频扩展16 文章1 文章201 文章401 文章601 文章801 文章1001 资讯1 资讯501 资讯1001 资讯1501 标签1 标签501 标签1001 关键词1 关键词501 关键词1001 关键词1501 专题2001
快速傅里叶变换FFT算法-精简版
2025-09-30 23:15:06 责编:小OO
文档
1  一维DFT的快速算法—FFT

当序列的点数不超过时,它的点定义为

                                   (1)

反变换定义为

                                 (2)

二者形式相似,快速算法的原理一样,这里先就其正变换进行讨论。令,当依次取为时,可表示为如下的方程组:

              (3)

由上式可见,直接按照定义计算点序列的点时,每行含个复乘和个加,从而直接按定义计算点的总计算量为个复乘和个加。当较大时,很大,计算量过大不仅耗时长,还会因字长有限而产生较大的误差,甚至造成计算结果不收敛。所谓快速傅里叶变换就是能大大减少计算量而完成全部点计算的算法。下面介绍两种经典的的快速算法:频域抽取的算法和时域抽取的算法。

1.1  频域抽取的基2算法

1.1.1 正变换的计算

这里仅介绍基2算法,即是2的整次幂的情况。由定义

                                     (4)

把分成两半,即和,代入(4)式得

              (5)

由于

(5)式两项又可合并为

                   (6)

当为偶数时,注意到, ,(6)式变为

                                 (7)

当为奇数时,,(6)式变为

                      (8)

这样就把一个点序列()的点()计算化成了两个点序列(和)的点(和)计算。由划分成

和的计算量为个加,即

和个乘,即

由于算出的点,是的点()中为偶数的那一半,由算出的则是为偶数的那一半,故需要把偶数的抽出来放在一起作为的()输出,同时把奇数的抽出来放在一起作为的()输出。由于是频域变量,故这种算法称为频域抽取的算法。

接着,两个点仍可用上述方法各经个乘个加划分成两个点(同时还要做相应的频域抽取),从而共划分成4个点,总划分计算量仍是个加和个乘。这样的划分可一步步做下去,不难看出,每步的总划分计算量都是个加和个乘。

经过步的划分后就划成了个如下两点的计算问题

                                      (9)

上式所需计算量是2个加和1个乘,于是完成个两点的总计算量仍是个加和个乘。从而完成全部点的总计算量个加和个乘,这比直接按定义计算所需的个乘和加要少得多。例如,,,用算法计算所需的乘法个数为,而直接按定义计算所需的乘法个数为,二者相差倍。若直接计算需半小时,而用计算只需9s即可完成,可见其效率之高,而且越大,的效率提高越明显。

图1  频域抽取的8点FFT计算流图

一般情况下,由于做了次分奇偶的抽取,此算法最后的个两点计算出的不是顺序抽取的。次序的变化可用二进码来说明:第一次抽取所分的奇偶是由二进码第1位是1或0来区分的,该位为0时为偶数,该位为1时为奇数,第二次抽奇偶是由二进码第2位是1或0来区分的……,每次抽取都是把偶数项放在前(左)边,把奇数项放在后(右)边,从而抽取以后数的二进码是按照二进制位从左向右依次排列的,和普通二进制数从右向左依次排列的的规律正好相反,所以称为倒位序。在计算出之后要把倒位序变成顺序。

1.1.2 逆变换的计算

所谓逆变换是指由求的计算,若直接按定义

做计算,则除了求和号和正变换相同的计算量外,每算一个都还需再多做一个乘的乘法运算。故按定义完成全部点的总计算量是个加和个乘。下面从图导出它的快速算法,先讨论第3列的2点的逆运算如何完成。由式(7)得,

由上式不难解出

                                                 (10)

图2  频域抽取的8点IFFT计算流图

此计算过程如图2所示,可以看出:左边各列的划分计算也都是由个碟形运算来完成的,只是各碟形运算所乘的相移因子不同。把每个碟形运算都用图的办法变成对应的逆运算,并把它们按输入在左、输出在右重新排列,就得到了全部点的计算流图。给出了的示例,图中先对顺序输入的做次的频域抽取,并把个乘的运算合成了一个乘的运算放在了最前边,然后就开始做求逆的碟形运算。

1.2 时域抽取的基2算法

比较正变换和反变换的定义式

可见,去掉乘的运算,把换成,交换、和、,反变换定义式就变成了正变换的定义式。对图2做这些变换,则得到图3的流程图。对图1做这些变换,则得到图4的流程图。这就是时域抽取的算法流图。进行碟形运算之前,先要对顺序的时域输入序列进行次的奇偶抽取,故称为时域抽取的算法。

图3  时域抽取的8点FFT计算流图

比较图2和图3不难看出,两种算法的计算量是完全一样的。这里先算出个两点的

                    (11)

图4  时域抽取的8点IFFT计算流图

然后把个两点的组合成个4点的,再把个4点的组合成个8点的,经过次的组合之后,就得到了顺序点计算结果。算法原理参考文献【4】。

1.3  一维基2 FFT算法编程

以频域抽取的基2 FFT正变换为例,对FFT的信号流图进行讨论,以找到FFT算法的规律。

1)分级

在进行变换的过程中,从点到两点共分了M级,如图1所示,从左到右依次为级,级,…,级。

2)倒位序

在频域抽取的基2 FFT算法中,输出数据不是按照序列的先后顺序排列的,这是由于变换过程中,输出按奇、偶抽取的缘故。如果将序列中标号用二进制值表示,那么在FFT信号流图输入端,位于处,称为倒序。以8点FFT为例,顺序和倒序的关系如表1所示。

表1  顺序和倒序对照表

顺序倒序
十进制数二进制数二进制数十进制数
0

1

2

3

4

5

6

7

0  0  0

0  0  1

0  1  0

0  1  1

1  0  0

1  0  1

1  1  0

1  1  1

0  0  0

1  0  0

0  1  0

1  1  0

0  0  1

1  0  1

0  1  1

1  1  1

0

4

2

6

1

5

3

7

从表1可以看出,一个自然顺序二进制数,是在最低位加1,逢2向左移位;而倒序数的顺序是在最高位加1,逢2向右移位。用表示顺序数,表示倒序数,表示位权重。对于一个倒序数来说,下一个倒序数可以按下面的方法求:先对最高位加1,相当于十进制运算。如果,说明二进制最高位为0,则直接由得到下一个倒序值;如果,说明二进制最高

位为1,则将的最高位变为0,通过实现,同时令,接着判断次高位是否为0,直到位为0时,令。归结起来算法流程图如图5所示。

输出

     

     图5  倒位序程序流程图              图6  频域抽取FFT程序流程图

3)蝶形运算单元和同址计算

频域抽取的信号流图中,基本的运算结构如图7所示,该运算结构的形状像蝴蝶,故称为“蝶形运算单元”。

在这一结构中,DFT和IDFT运算关系分别为

         (12)

(a) 两点DFT的计算 

(b) 两点IDFT的计算

图7  频域抽取FFT算法流图中第m级碟形单元

而在时域抽取的信号流图中,基本的运算结构如图8所示。在这一结构中,DFT和IDFT运算关系分别为

           (13)

(a) 两点DFT的计算

(b) 两点IDFT的计算 

图8  时域抽取FFT算法流图中第m级碟形单元

其中,、分别表示该蝶形运算单元的上下节点的序号。可以看出参与运算的输入序号,输出序号仍为,并且该运算不涉及到其它的点,因此我们可以将输出的结果仍然放在数组中,称这样的操作为同址计算。也就是说,共同占有同一个存储单元。

4)寻址和相移因子的计算

时域抽取基2 FFT信号流图中,每一级有个蝶形单元。每一级的个蝶形单元又可以分为若干组,每一组具有相同的结构和因子的分布。

如图1所示,第1级分为1组,第2级分为2组,…,第级分为组。在第级中,相邻组之间的间距(也即每个分组所含节点数)为,每个蝶形单元的上下节点之间的距离(也即每个分组所含碟形单元数)为。每组的相移因子

,其中

综合以上各步骤,得到频域抽取FFT程序流程图如图6所示。采用类似的步骤可得到频域抽取IFFT流程图、时域抽取FFT与IFFT流程图(略)。

频域抽取FFT算法、时域抽取FFT算法的Visual C++源程序分别见附录1.(1),1.(2)。在Matlab中,傅里叶变换及其逆变换分别用dwt()和idwt()函数实现。

2  一维FFT算法的应用

2.1 利用FFT计算连续时间信号的傅里叶变换

设是连续时间信号,并假设时,则其傅里叶变换由下式给出

令是一个固定的正实数,是一个固定的正整数。当时,利用FFT算法可计算。

已知一个固定的时间间隔,选择足够小,使得每一个秒的间隔内,的变化很小,则式中积分可近似为

     

                                            (27)

假设足够大,对于所有的整数,幅值很小,则式(27)变为

                                       (28)

当时,式(28)两边的值为

                  (29)

其中代表抽样信号的点。最后令,则上式变为

                           (30)

首先用FFT算法求出,然后可用上式求出时的。

应该强调的是,式(28)只是一个近似表示,计算得到的只是一个近似值。通过取更小的抽样间隔,或者增加点数,可以得到更精确的值。如果时,幅度谱很小,对应于奈奎斯特抽样频率,抽样间隔选择比较合适。如果已知信号只在时间区间内存在,可以通过对时的抽样信号补零,使足够大。

例1  利用FFT计算傅里叶变换

如图12所示的信号

其傅里叶变换为:

利用下面的命令,可得到的近似值和准确值。                             图12  连续时间信号x(t)

              

N=input('Input N:');

T=input('Input T:');

%计算X(w)近似值

t=0:T:2;

x=[t-1 zeros(1,N-length(t))];

X=fft(x);

gamma=2*pi/(N*T);

k=0:10/gamma;

Xapp=(1-exp(-i*k*gamma*T))/(i*k*gamma)*X;

%计算真实值X(w)

w=0.05:0.05:10;

Xact=exp(-i*w)*2*i.*(w.*cos(w)-sin(w))./(w.*w);

plot(k*gamma,abs(Xapp(1:length(k))),'o',w,abs(Xact));

legend('近似值','真实值');

xlabel('频率(rad/s)');ylabel('|X|')

运行程序后输入N=128,T=0.1,此时,得到实际的和近似的傅里叶变换的幅度谱如图13所示,此时近似值已经相当准确。通过增加NT可以增加更多的细节,减少T使得到的值更精确。再次运行程序后输入N=512,T=0.05,此时,得到实际的和近似的傅里叶变换的幅度谱如图14所示。

图13  N=128,T=0.1时的幅度谱

图14  N=512,T=0.05时的幅度谱

2.2 利用FFT对离散信号进行滤波

利用FFT算法对信号进行滤波的步骤如下:1)通过采样将信号离散化;2)对离散化信号进行傅里叶变换;3)对变换后的系数进行处理,将绝对值大于某一阈值的系数置为0,保留剩余的系数;4)利用IFFT算法对处理后的信号进行逆傅里叶变换。

例4  股票价格分析

首先进入网址http://finance.yahoo.com/q/hp?s=YHOO,点击网页底部位置的Download To Spreadsheet按钮,即可把以Excel表格格式存储的价格数据下载到本地计算机。表格从1列至第6列分别给出了从1996年4月12日至2007年5月30日的交易期里每天的开盘价、最高价、最低价、收盘价、成交量以及趋势。数据下载完成后,需要颠倒顺序,使得最早时间的数据首先显示。然后另存到Matlab所在的目录中,并重新命名为“yhoodata.csv”。 本次分析选择开盘价,时间是从2007年1月1日至2007年5月30日的=102个交易日期。

令代表一支股票的开盘价。为了便于分析,可以先从中减去跃变,得到,即

                   (32)

输入以下命令,得到的频谱如图19所示。

o=csvread('yhoodata.csv',2700,1,[2700 1 2801 1])';

N=102;

for n=1:N

    x(n)=o(n)-o(1)-((o(N)-o(1))/(N-1))*(n-1);

end

X=fft(x);

k=0:N-1;

stem(k,abs(X));

axis([0 101 0 300]);

xlabel('k');ylabel('|X[k]|')

图19  x[n]的幅度谱

可以看出上图中有5个较强谱分量,频率()分别对应和。保留这5个频率分量的系数,将其他频率分量的系数置为0,然后再进行逆傅里叶变换,得到滤波后的近似值。输入如下Matlab程序,得到真实值与滤波后的近似值,如图20所示。

plot(x);hold on;

fliterX=[X(1:2) 0 0 X(5) zeros(1,102-9) X(99) 0 0 X(102)];

fliterx=ifft(fliterX);

plot(real(fliterx),'r');

axis([1 102 0 7]);

xlabel('n');ylabel('x[n]的值和它的近似值');

legend('x[n]真实值','x[n]近似值')

图20  x[n]的真实值与滤波后的近似值

从上图可以看出,滤波后的近似值既大致上保留了真实值的变化趋势,而且与其十分接近。与滤波前比较,滤波后的图形要比滤波前平滑得多。再由式(32)即可求得

                     (33)

输入如下Matlab程序,画出真实开盘价与近似开盘价的图形。如图21所示,可以看出是近似基础上的平滑值。

plot(o);hold on;

for n=1:N

    oapp(n)=fliterx(n)+o(1)+((o(N)-o(1))/(N-1))*(n-1);

end

plot(oapp,'r');

axis([1 102 25 34]);

xlabel('n');ylabel('o[n]的值和它的近似值');

legend('o[n]真实值','o[n]近似值')

图21  o[n]的真实值与滤波后的近似值

附 录:一维DFT的基2 FFT 算法Visual C++程序

设原始信号由一个模拟频率为3265Hz的正弦信号,以8000Hz的频率采样获得

(1) 频域抽取的FFT和IFFT算法

#include

#include

#include

#define pi 3.14159265

#define N 

typedef std::complex complex;

/*----频域抽取的FFT算法----*/

void fft(complex f[])

{

    int i,j,k,m,n,l,r,M;

    int la,lb,lc;

    complex t;

    /*----计算分解的级数M=log2(N)----*/

    for(i=N,M=1;(i=i/2)!=1;M++); 

    /*----FFT算法----*/

for(m=1;m<=M;m++)

    {

        la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数    

        lb=la/2;    //lb代表第m级每个分组所含碟形单元数

                   //同时它也表示每个碟形单元上下节点之间的距离

        /*----碟形运算----*/

        for(l=1;l<=lb;l++)

        {

            r=(l-1)*pow(2,m-1);

            for(n=l-1;n            {

                lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     

                t=f[n]+f[lc];

                f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),-sin(2*pi*r/N));

                f[n]=t;

            }

        }

    }

    /*----按照倒位序重新排列变换后信号----*/

for(i=1,j=N/2;i<=N-2;i++)

    {

        if(i        {

            t=f[j];

            f[j]=f[i];

            f[i]=t;

        }

        k=N/2;

        while(k<=j)

        {

            j=j-k;

            k=k/2;

        }

        j=j+k;

    }

}

/*----频域抽取的IFFT算法----*/

void ifft(complex f[])

{

    int i,j,k,m,n,l,r,M;

    int la,lb,lc;

    complex t;

    /*----计算分解的级数M=log2(N)----*/

    for(i=N,M=1;(i=i/2)!=1;M++); 

    

    /*----按照倒位序重新排列原信号----*/

for(i=1,j=N/2;i<=N-2;i++)

    {

        if(i        {

            t=f[j];

            f[j]=f[i];

            f[i]=t;

        }

        k=N/2;

        while(k<=j)

        {

            j=j-k;

            k=k/2;

        }

        j=j+k;

    }

    /*----将信号乘以1/N----*/

for(i=0;i    /*----IFFT算法----*/

for(m=1;m<=M;m++)

    {

        la=pow(2,m); //la=2^m代表第m级每个分组所含节点数        

        lb=la/2;    //lb代表第m级每个分组所含碟形单元数

                   //同时它也表示每个碟形单元上下节点之间的距离

        /*----碟形运算----*/

        for(l=1;l<=lb;l++)

        {

            r=(l-1)*pow(2,M-m);

            for(n=l-1;n            {

                lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     

                t=f[lc]*complex(cos(2*pi*r/N),sin(2*pi*r/N));

                f[lc]=f[n]-t;

                f[n]=f[n]+t;

            }

        }

    }

}

/*----显示信号数据----*/

void display(complex f[])

{

    int i;

for(i=0;i    {

        cout.width(9);

        cout.setf(ios::fixed);

        cout.precision(4);

        cout<        cout.flags(ios::showpos);

        cout.width(9);

        cout.setf(ios::fixed);

        cout.precision(4);

        cout<        cout.unsetf(ios::showpos);

        if((i+1)%3==0) cout<    }

}

/*----主函数----*/

void main()

{

    int i;

    complex f[N];

for(i=0;i        f[i]=complex(0.8*sin(2*pi*3265*i/8000),0.0);

    cout<    display(f);

    fft(f);

    cout<    display(f);

    ifft(f);

    cout<    display(f);

}

运行结果为:

原信号:

   0.0000  +0.0000i        0.4366  +0.0000i       -0.7317  +0.0000i

   0.77  +0.0000i       -0.5917  +0.0000i        0.2020  +0.0000i

   0.2532  +0.0000i       -0.6263  +0.0000i        0.79  +0.0000i

  -0.7085  +0.0000i        0.3909  +0.0000i        0.0534  +0.0000i

  -0.4803  +0.0000i        0.7516  +0.0000i       -0.7793  +0.0000i

   0.5545  +0.0000i       -0.1499  +0.0000i       -0.3032  +0.0000i

   0.6581  +0.0000i       -0.7997  +0.0000i        0.6821  +0.0000i

  -0.3435  +0.0000i       -0.1065  +0.0000i        0.5219  +0.0000i

  -0.7682  +0.0000i        0.7656  +0.0000i       -0.5148  +0.0000i

   0.0971  +0.0000i        0.3520  +0.0000i       -0.6870  +0.0000i

   0.7994  +0.0000i       -0.6527  +0.0000i        0.2945  +0.0000i

   0.1592  +0.0000i       -0.5612  +0.0000i        0.7814  +0.0000i

  -0.7484  +0.0000i        0.4728  +0.0000i       -0.0440  +0.0000i

  -0.3991  +0.0000i        0.7128  +0.0000i       -0.7955  +0.0000i

   0.6204  +0.0000i       -0.2442  +0.0000i       -0.2111  +0.0000i

   0.5980  +0.0000i       -0.7911  +0.0000i        0.7278  +0.0000i

  -0.4287  +0.0000i       -0.0094  +0.0000i        0.4445  +0.0000i

  -0.7354  +0.0000i        0.7881  +0.0000i       -0.5853  +0.0000i

   0.1929  +0.0000i        0.2621  +0.0000i       -0.6321  +0.0000i

   0.7973  +0.0000i       -0.7041  +0.0000i        0.3826  +0.0000i

   0.0628  +0.0000i       -0.4878  +0.0000i        0.7548  +0.0000i

  -0.7772  +0.0000i

FFT变换后的信号:

  -0.2416  +0.0000i       -0.2415  -0.0146i       -0.2413  -0.0294i

  -0.2409  -0.0443i       -0.2402  -0.0595i       -0.2394  -0.0751i

  -0.2384  -0.0911i       -0.2371  -0.1078i       -0.2355  -0.1253i

  -0.2336  -0.1438i       -0.2314  -0.1634i       -0.2286  -0.1844i

  -0.2253  -0.2072i       -0.2214  -0.2322i       -0.2165  -0.2600i

  -0.2106  -0.2911i       -0.2032  -0.3268i       -0.1939  -0.3683i

  -0.1818  -0.4177i       -0.1658  -0.4784i       -0.1439  -0.5557i

  -0.1124  -0.6588i       -0.03  -0.8062i        0.0168  -1.0398i

   0.1783  -1.4797i        0.6372  -2.6746i        8.8462 -23.4497i

  -1.6196  +2.9360i       -0.9624  +1.2195i       -0.7711  +0.6680i

  -0.6881  +0.3740i       -0.6501  +0.1707i       -0.63  +0.0000i

  -0.6501  -0.1707i       -0.6881  -0.3740i       -0.7711  -0.6680i

  -0.9624  -1.2195i       -1.6196  -2.9360i        8.8462 +23.4497i

   0.6372  +2.6746i        0.1783  +1.4797i        0.0168  +1.0398i

  -0.03  +0.8062i       -0.1124  +0.6588i       -0.1439  +0.5557i

  -0.1658  +0.4784i       -0.1818  +0.4177i       -0.1939  +0.3683i

  -0.2032  +0.3268i       -0.2106  +0.2911i       -0.2165  +0.2600i

  -0.2214  +0.2322i       -0.2253  +0.2072i       -0.2286  +0.1844i

  -0.2314  +0.1634i       -0.2336  +0.1438i       -0.2355  +0.1253i

  -0.2371  +0.1078i       -0.2384  +0.0911i       -0.2394  +0.0751i

  -0.2402  +0.0595i       -0.2409  +0.0443i       -0.2413  +0.0294i

  -0.2415  +0.0146i

IFFT变换后的信号:

   0.0000  -0.0000i        0.4366  +0.0000i       -0.7317  +0.0000i

   0.77  -0.0000i       -0.5917  -0.0000i        0.2020  +0.0000i

   0.2532  -0.0000i       -0.6263  +0.0000i        0.79  -0.0000i

  -0.7085  +0.0000i        0.3909  -0.0000i        0.0534  -0.0000i

  -0.4803  +0.0000i        0.7516  -0.0000i       -0.7793  +0.0000i

   0.5545  +0.0000i       -0.1499  +0.0000i       -0.3032  -0.0000i

   0.6581  -0.0000i       -0.7997  +0.0000i        0.6821  -0.0000i

  -0.3435  -0.0000i       -0.1065  -0.0000i        0.5219  -0.0000i

  -0.7682  +0.0000i        0.7656  -0.0000i       -0.5148  +0.0000i

   0.0971  +0.0000i        0.3520  +0.0000i       -0.6870  +0.0000i

   0.7994  -0.0000i       -0.6527  -0.0000i        0.2945  -0.0000i

   0.1592  +0.0000i       -0.5612  +0.0000i        0.7814  +0.0000i

  -0.7484  +0.0000i        0.4728  +0.0000i       -0.0440  +0.0000i

  -0.3991  +0.0000i        0.7128  -0.0000i       -0.7955  +0.0000i

   0.6204  -0.0000i       -0.2442  -0.0000i       -0.2111  -0.0000i

   0.5980  -0.0000i       -0.7911  -0.0000i        0.7278  -0.0000i

  -0.4287  +0.0000i       -0.0094  -0.0000i        0.4445  -0.0000i

  -0.7354  +0.0000i        0.7881  +0.0000i       -0.5853  -0.0000i

   0.1929  -0.0000i        0.2621  -0.0000i       -0.6321  +0.0000i

   0.7973  -0.0000i       -0.7041  +0.0000i        0.3826  +0.0000i

   0.0628  -0.0000i       -0.4878  +0.0000i        0.7548  +0.0000i

  -0.7772  -0.0000i

(2) 时域抽取的FFT和IFFT算法

#include

#include

#include

#define pi 3.14159265

#define N 

typedef std::complex complex;

/*----时域抽取的FFT算法----*/

void fft(complex f[])

{

    int i,j,k,m,n,l,r,M;

    int la,lb,lc;

    complex t;

    /*----计算分解的级数M=log2(N)----*/

    for(i=N,M=1;(i=i/2)!=1;M++); 

    

    /*----按照倒位序重新排列原信号----*/

for(i=1,j=N/2;i<=N-2;i++)

    {

        if(i        {

            t=f[j];

            f[j]=f[i];

            f[i]=t;

        }

        k=N/2;

        while(k<=j)

        {

            j=j-k;

            k=k/2;

        }

        j=j+k;

    }

    /*----FFT算法----*/

for(m=1;m<=M;m++)

    {

        la=pow(2,m); //la=2^m代表第m级每个分组所含节点数        

        lb=la/2;    //lb代表第m级每个分组所含碟形单元数

                   //同时它也表示每个碟形单元上下节点之间的距离

        /*----碟形运算----*/

        for(l=1;l<=lb;l++)

        {

            r=(l-1)*pow(2,M-m);

            for(n=l-1;n            {

                lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     

                t=f[lc]*complex(cos(2*pi*r/N),-sin(2*pi*r/N));

                f[lc]=f[n]-t;

                f[n]=f[n]+t;

            }

        }

    }

}

/*----时域抽取的IFFT算法----*/

void ifft(complex f[])

{

    int i,j,k,m,n,l,r,M;

    int la,lb,lc;

    complex t;

    /*----计算分解的级数M=log2(N)----*/

    for(i=N,M=1;(i=i/2)!=1;M++); 

    /*----将信号乘以1/N----*/

for(i=0;i    /*----IFFT算法----*/

for(m=1;m<=M;m++)

    {

        la=pow(2,M+1-m); //la=2^m代表第m级每个分组所含节点数    

        lb=la/2;    //lb代表第m级每个分组所含碟形单元数

                   //同时它也表示每个碟形单元上下节点之间的距离

        /*----碟形运算----*/

        for(l=1;l<=lb;l++)

        {

            r=(l-1)*pow(2,m-1);

            for(n=l-1;n            {

                lc=n+lb;  //n,lc分别代表一个碟形单元的上、下节点编号     

                t=f[n]+f[lc];

                f[lc]=(f[n]-f[lc])*complex(cos(2*pi*r/N),sin(2*pi*r/N));

                f[n]=t;

            }

        }

    }

    /*----按照倒位序重新排列变换后信号----*/

for(i=1,j=N/2;i<=N-2;i++)

    {

        if(i        {

            t=f[j];

            f[j]=f[i];

            f[i]=t;

        }

        k=N/2;

        while(k<=j)

        {

            j=j-k;

            k=k/2;

        }

        j=j+k;

    }

}

/*----显示信号数据----*/

void display(complex f[])

{

    int i;

for(i=0;i    {

        cout.width(9);

        cout.setf(ios::fixed);

        cout.precision(4);

        cout<        cout.flags(ios::showpos);

        cout.width(9);

        cout.setf(ios::fixed);

        cout.precision(4);

        cout<        cout.unsetf(ios::showpos);

        if((i+1)%3==0) cout<    }

}

/*----主函数----*/

void main()

{

    int i;

    complex f[N];

for(i=0;i        f[i]=complex(0.8*sin(2*pi*3265*i/8000),0.0);

    cout<    display(f);

    fft(f);

    cout<    display(f);

    ifft(f);

    cout<    display(f);

}

运行结果为:

原信号:

   0.0000  +0.0000i        0.4366  +0.0000i       -0.7317  +0.0000i

   0.77  +0.0000i       -0.5917  +0.0000i        0.2020  +0.0000i

   0.2532  +0.0000i       -0.6263  +0.0000i        0.79  +0.0000i

  -0.7085  +0.0000i        0.3909  +0.0000i        0.0534  +0.0000i

  -0.4803  +0.0000i        0.7516  +0.0000i       -0.7793  +0.0000i

   0.5545  +0.0000i       -0.1499  +0.0000i       -0.3032  +0.0000i

   0.6581  +0.0000i       -0.7997  +0.0000i        0.6821  +0.0000i

  -0.3435  +0.0000i       -0.1065  +0.0000i        0.5219  +0.0000i

  -0.7682  +0.0000i        0.7656  +0.0000i       -0.5148  +0.0000i

   0.0971  +0.0000i        0.3520  +0.0000i       -0.6870  +0.0000i

   0.7994  +0.0000i       -0.6527  +0.0000i        0.2945  +0.0000i

   0.1592  +0.0000i       -0.5612  +0.0000i        0.7814  +0.0000i

  -0.7484  +0.0000i        0.4728  +0.0000i       -0.0440  +0.0000i

  -0.3991  +0.0000i        0.7128  +0.0000i       -0.7955  +0.0000i

   0.6204  +0.0000i       -0.2442  +0.0000i       -0.2111  +0.0000i

   0.5980  +0.0000i       -0.7911  +0.0000i        0.7278  +0.0000i

  -0.4287  +0.0000i       -0.0094  +0.0000i        0.4445  +0.0000i

  -0.7354  +0.0000i        0.7881  +0.0000i       -0.5853  +0.0000i

   0.1929  +0.0000i        0.2621  +0.0000i       -0.6321  +0.0000i

   0.7973  +0.0000i       -0.7041  +0.0000i        0.3826  +0.0000i

   0.0628  +0.0000i       -0.4878  +0.0000i        0.7548  +0.0000i

  -0.7772  +0.0000i

FFT变换后的信号:

  -0.2416  +0.0000i       -0.2415  -0.0146i       -0.2413  -0.0294i

  -0.2409  -0.0443i       -0.2402  -0.0595i       -0.2394  -0.0751i

  -0.2384  -0.0911i       -0.2371  -0.1078i       -0.2355  -0.1253i

  -0.2336  -0.1438i       -0.2314  -0.1634i       -0.2286  -0.1844i

  -0.2253  -0.2072i       -0.2214  -0.2322i       -0.2165  -0.2600i

  -0.2106  -0.2911i       -0.2032  -0.3268i       -0.1939  -0.3683i

  -0.1818  -0.4177i       -0.1658  -0.4784i       -0.1439  -0.5557i

  -0.1124  -0.6588i       -0.03  -0.8062i        0.0168  -1.0398i

   0.1783  -1.4797i        0.6372  -2.6746i        8.8462 -23.4497i

  -1.6196  +2.9360i       -0.9624  +1.2195i       -0.7711  +0.6680i

  -0.6881  +0.3740i       -0.6501  +0.1707i       -0.63  +0.0000i

  -0.6501  -0.1707i       -0.6881  -0.3740i       -0.7711  -0.6680i

  -0.9624  -1.2195i       -1.6196  -2.9360i        8.8462 +23.4497i

   0.6372  +2.6746i        0.1783  +1.4797i        0.0168  +1.0398i

  -0.03  +0.8062i       -0.1124  +0.6588i       -0.1439  +0.5557i

  -0.1658  +0.4784i       -0.1818  +0.4177i       -0.1939  +0.3683i

  -0.2032  +0.3268i       -0.2106  +0.2911i       -0.2165  +0.2600i

  -0.2214  +0.2322i       -0.2253  +0.2072i       -0.2286  +0.1844i

  -0.2314  +0.1634i       -0.2336  +0.1438i       -0.2355  +0.1253i

  -0.2371  +0.1078i       -0.2384  +0.0911i       -0.2394  +0.0751i

  -0.2402  +0.0595i       -0.2409  +0.0443i       -0.2413  +0.0294i

  -0.2415  +0.0146i

IFFT变换后的信号:

   0.0000  -0.0000i        0.4366  +0.0000i       -0.7317  -0.0000i

   0.77  +0.0000i       -0.5917  -0.0000i        0.2020  -0.0000i

   0.2532  +0.0000i       -0.6263  -0.0000i        0.79  -0.0000i

  -0.7085  -0.0000i        0.3909  -0.0000i        0.0534  +0.0000i

  -0.4803  -0.0000i        0.7516  +0.0000i       -0.7793  +0.0000i

   0.5545  -0.0000i       -0.1499  +0.0000i       -0.3032  -0.0000i

   0.6581  -0.0000i       -0.7997  -0.0000i        0.6821  +0.0000i

  -0.3435  +0.0000i       -0.1065  -0.0000i        0.5219  +0.0000i

  -0.7682  +0.0000i        0.7656  +0.0000i       -0.5148  +0.0000i

   0.0971  -0.0000i        0.3520  +0.0000i       -0.6870  -0.0000i

   0.7994  -0.0000i       -0.6527  +0.0000i        0.2945  -0.0000i

   0.1592  +0.0000i       -0.5612  -0.0000i        0.7814  +0.0000i

  -0.7484  +0.0000i        0.4728  -0.0000i       -0.0440  +0.0000i

  -0.3991  -0.0000i        0.7128  -0.0000i       -0.7955  -0.0000i

   0.6204  -0.0000i       -0.2442  +0.0000i       -0.2111  -0.0000i

   0.5980  +0.0000i       -0.7911  +0.0000i        0.7278  -0.0000i

  -0.4287  +0.0000i       -0.0094  -0.0000i        0.4445  +0.0000i

  -0.7354  -0.0000i        0.7881  -0.0000i       -0.5853  +0.0000i

   0.1929  -0.0000i        0.2621  +0.0000i       -0.6321  +0.0000i

   0.7973  +0.0000i       -0.7041  +0.0000i        0.3826  -0.0000i

   0.0628  +0.0000i       -0.4878  -0.0000i        0.7548  -0.0000i

  -0.7772  +0.0000i 下载本文

显示全文
专题