守望者--AIR技术交流

 找回密码
 立即注册

QQ登录

只需一步,快速开始

扫一扫,访问微社区

搜索
热搜: ANE FlasCC 炼金术
查看: 712|回复: 3

[音频分析] FFT(快速傅立叶算法 各种编程语言实现)

[复制链接]
  • TA的每日心情
    擦汗
    11 小时前
  • 签到天数: 433 天

    [LV.9]以坛为家II

    1739

    主题

    2091

    帖子

    12万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    457
    贡献
    24
    金币
    50924
    钢镚
    1419

    开源英雄守望者

    发表于 2015-9-10 18:02:26 | 显示全部楼层 |阅读模式
    Java


    1. package com.test.test2;

    2. public class FFT {
    3.     public static final int FFT_N_LOG = 10; // FFT_N_LOG <= 13
    4.     public static final int FFT_N = 1 << FFT_N_LOG;
    5.     private static final float MINY = (float) ((FFT_N << 2) * Math.sqrt(2)); // (*)
    6.     private final float[] real, imag, sintable, costable;
    7.     private final int[] bitReverse;

    8.     public FFT() {
    9.         real = new float[FFT_N];
    10.         imag = new float[FFT_N];
    11.         sintable = new float[FFT_N >> 1];
    12.         costable = new float[FFT_N >> 1];
    13.         bitReverse = new int[FFT_N];

    14.         int i, j, k, reve;
    15.         for (i = 0; i < FFT_N; i++) {
    16.             k = i;
    17.             for (j = 0, reve = 0; j != FFT_N_LOG; j++) {
    18.                 reve <<= 1;
    19.                 reve |= (k & 1);
    20.                 k >>>= 1;
    21.             }
    22.             bitReverse[i] = reve;
    23.         }

    24.         double theta, dt = 2 * 3.14159265358979323846 / FFT_N;
    25.         for (i = 0; i < (FFT_N >> 1); i++) {
    26.             theta = i * dt;
    27.             costable[i] = (float) Math.cos(theta);
    28.             sintable[i] = (float) Math.sin(theta);
    29.         }
    30.     }

    31.     /**
    32.      * 用于频谱显示的快速傅里叶变换
    33.      *
    34.      * @param realIO
    35.      *            输入FFT_N个实数,也用它暂存fft后的FFT_N/2个输出值(复数模的平方)。
    36.      */
    37.     public void calculate(float[] realIO) {
    38.         int i, j, k, ir, exchanges = 1, idx = FFT_N_LOG - 1;
    39.         float cosv, sinv, tmpr, tmpi;
    40.         for (i = 0; i != FFT_N; i++) {
    41.             real[i] = realIO[bitReverse[i]];
    42.             imag[i] = 0;
    43.         }

    44.         for (i = FFT_N_LOG; i != 0; i--) {
    45.             for (j = 0; j != exchanges; j++) {
    46.                 cosv = costable[j << idx];
    47.                 sinv = sintable[j << idx];
    48.                 for (k = j; k < FFT_N; k += exchanges << 1) {
    49.                     ir = k + exchanges;
    50.                     tmpr = cosv * real[ir] - sinv * imag[ir];
    51.                     tmpi = cosv * imag[ir] + sinv * real[ir];
    52.                     real[ir] = real[k] - tmpr;
    53.                     imag[ir] = imag[k] - tmpi;
    54.                     real[k] += tmpr;
    55.                     imag[k] += tmpi;
    56.                 }
    57.             }
    58.             exchanges <<= 1;
    59.             idx--;
    60.         }

    61.         j = FFT_N >> 1;
    62.         /*
    63.          * 输出模的平方(的FFT_N倍):
    64.          * for(i = 1; i <= j; i++)
    65.          * realIO[i-1] = real[i] * real[i] + imag[i] * imag[i];
    66.          *
    67.          * 如果FFT只用于频谱显示,可以"淘汰"幅值较小的而减少浮点乘法运算. MINY的值
    68.          * 和Spectrum.Y0,Spectrum.logY0对应.
    69.          */
    70.         sinv = MINY;
    71.         cosv = -MINY;
    72.         for (i = j; i != 0; i--) {
    73.             tmpr = real[i];
    74.             tmpi = imag[i];
    75.             if (tmpr > cosv && tmpr < sinv && tmpi > cosv && tmpi < sinv)
    76.                 realIO[i - 1] = 0;
    77.             else
    78.                 realIO[i - 1] = tmpr * tmpr + tmpi * tmpi;
    79.         }
    80.     }

    81.     public static void main(String[] args) {
    82.         FFT fft2 = new FFT();
    83.         float[] realIo = { 1, 2 };
    84.         fft2.calculate(realIo);
    85.     }
    86. }
    复制代码


    守望者AIR技术交流社区(www.airmyth.com)
    回复

    使用道具 举报

  • TA的每日心情
    擦汗
    11 小时前
  • 签到天数: 433 天

    [LV.9]以坛为家II

    1739

    主题

    2091

    帖子

    12万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    457
    贡献
    24
    金币
    50924
    钢镚
    1419

    开源英雄守望者

     楼主| 发表于 2015-9-10 18:29:22 | 显示全部楼层
    1. #include<graphics.h>
    2. #include<stdio.h>
    3. #include<math.h>
    4. #include<conio.h>
    5. #include<dos.h>

    6. float ar[2048],ai[2048];
    7. float a[4100];

    8. void testdata()
    9. {
    10. int i;
    11.    for(i=0;i<512;i++)
    12.    {
    13.        ar[i]=50*cos(i*3.14159/32)+60*cos(i*3.14159/16)+53*cos(i*3.14159/64)+24*cos(i*3.14159/8)+10*cos(i*3.14159/128);
    14.        ai[i]=0;
    15.    }
    16. }

    17. void  fft(int nn,int t)
    18. {
    19.       int n1,n2,i,j,k,l,m,s,l1;
    20.       float t1,t2,x,y;
    21.       float w1,w2,u1,u2,z;
    22.       float pi=3.14159;
    23.       switch(nn)
    24.       {
    25.           case 2048:  s=11; break;
    26.           case 1024:  s=10; break;
    27.           case 512:   s=9;  break;
    28.           case 256:   s=8;  break;
    29.       }

    30.       n1=nn/2;  n2=nn-1;
    31.       j=1;
    32.       for(i=1;i<=nn;i++)
    33.       {
    34.                 a[2*i]=ar[i-1];
    35.                 a[2*i+1]=ai[i-1];
    36.       }
    37.       for(l=1;l<n2;l++)
    38.       {
    39.        if(l<j)
    40.        {
    41.         t1=a[2*j];
    42.         t2=a[2*j+1];
    43.         a[2*j]=a[2*l];
    44.         a[2*j+1]=a[2*l+1];
    45.         a[2*l]=t1;
    46.         a[2*l+1]=t2;
    47.        }
    48.        k=n1;
    49.        while (k<j)
    50.        {
    51.         j=j-k;
    52.         k=k/2;
    53.        }
    54.        j=j+k;
    55.      }
    56.      for(i=1;i<=s;i++)
    57.      {
    58.         u1=1;
    59.         u2=0;
    60.         m=(1<<i);
    61.         k=m>>1;
    62.         w1=cos(pi/k);
    63.         w2=t*sin(pi/k);
    64.         for(j=1;j<=k;j++)
    65.         {
    66.          for(l=j;l<nn;l=l+m)
    67.          {
    68.                 l1=l+k;
    69.                 t1=a[2*l1]*u1-a[2*l1+1]*u2;
    70.                 t2=a[2*l1]*u2+a[2*l1+1]*u1;
    71.                 a[2*l1]=a[2*l]-t1;
    72.                 a[2*l1+1]=a[2*l+1]-t2;
    73.                 a[2*l]=a[2*l]+t1;
    74.                 a[2*l+1]=a[2*l+1]+t2;
    75.          }
    76.          z=u1*w1-u2*w2;
    77.          u2=u1*w2+u2*w1;
    78.          u1=z;
    79.         }
    80.      }
    81.      for(i=1;i<=nn/2;i++)
    82.      {
    83.         ar[i]=a[2*i+2]/nn;
    84.         ai[i]=-a[2*i+3]/nn;
    85.         a[i]=4*sqrt(ar[i]*ar[i]+ai[i]*ai[i]);
    86.      }
    87. }

    88. int main ()
    89. {
    90.    int i;
    91.    int gdriver=DETECT,gmode;
    92.    initgraph(&gdriver,&gmode,"");
    93.    setfillstyle(SOLID_FILL,WHITE);
    94.    bar(0,0,639,479);
    95. //模拟测试数据   
    96.    testdata();
    97. //波形显示   
    98.    setcolor(BLACK);
    99.    line(10,119,550,119); // x轴
    100.    line(10,10,10,239); // y轴
    101.    setcolor(BLUE);
    102.    for(i=0;i<511;i++)
    103.      line(i+10,119-ar[i],i+10+1,119-ar[i+1]);
    104. //FFT   
    105.    fft(512,-1);
    106.    
    107. //频谱显示   
    108.    setcolor(BLACK);
    109.    line(10,459,550,459); // x轴
    110.    line(10,259,10,459); // y轴
    111.    setcolor(RED);
    112.    for(i=0;i<256;i++)
    113.      line(2*i+10,459-a[i],2*i+10,459);
    114.    getch();
    115.    closegraph();
    116.    return 0;
    117. }
    复制代码


    C++  
    守望者AIR技术交流社区(www.airmyth.com)
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    擦汗
    11 小时前
  • 签到天数: 433 天

    [LV.9]以坛为家II

    1739

    主题

    2091

    帖子

    12万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    457
    贡献
    24
    金币
    50924
    钢镚
    1419

    开源英雄守望者

     楼主| 发表于 2015-9-11 09:55:43 | 显示全部楼层
    2的整次幂长度FFT算法的C++实现
    头文件:
    1. /*
    2. * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
    3. *
    4. * This program is free software; you can redistribute it and/or modify it
    5. * under the terms of the GNU General Public License as published by the
    6. * Free Software Foundation, either version 2 or any later version.
    7. *
    8. * Redistribution and use in source and binary forms, with or without
    9. * modification, are permitted provided that the following conditions are met:
    10. *
    11. * 1. Redistributions of source code must retain the above copyright notice,
    12. *    this list of conditions and the following disclaimer.
    13. *
    14. * 2. Redistributions in binary form must reproduce the above copyright
    15. *    notice, this list of conditions and the following disclaimer in the
    16. *    documentation and/or other materials provided with the distribution.
    17. *
    18. * This program is distributed in the hope that it will be useful, but WITHOUT
    19. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    20. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
    21. * more details. A copy of the GNU General Public License is available at:
    22. * http://www.fsf.org/licensing/licenses
    23. */


    24. /*****************************************************************************
    25. *                                fftmr.h
    26. *
    27. * Fast Fourier Transform with Mixed Radix Algorithm
    28. *
    29. * This class is designed for calculating discrete Fourier transform and
    30. * inverse discrete Fourier transform of 1D signals by using Mixed-Radix
    31. * Algorithms. The length of signals should equal to powers of 2.
    32. *
    33. * The algorithm is modified from "kube-gustavson-fft.c". Which is a pretty
    34. * fast FFT algorithm. Not super-optimized lightning-fast, but very good
    35. * considering its age and relative simplicity.
    36. *
    37. * Zhang Ming, 2010-04, Xi'an Jiaotong University.
    38. *****************************************************************************/


    39. #ifndef FFTMR_H
    40. #define FFTMR_H


    41. #include <vector.h>


    42. namespace splab
    43. {

    44.     /**
    45.      * complex node converted from "C" version for this algorithm
    46.      */
    47.     template<typename Type>
    48.     struct Complex
    49.     {
    50.         Type re;
    51.         Type im;
    52.     };


    53.     /**
    54.      * two routines frequently used in FFT algorithm
    55.      */
    56.     bool    isPower2( int );
    57.     int     fastLog2( int );


    58.     /**
    59.      * Radix based FFT class
    60.      */
    61.         template<typename Type>
    62.         class FFTMR
    63.         {

    64.         public:

    65.                 FFTMR();
    66.                 ~FFTMR();

    67.         void fft( Vector< complex<Type> > &xn );
    68.         void ifft( Vector< complex<Type> > &Xk );
    69.                 void fft( const Vector<Type> &xn, Vector< complex<Type> > &Xk );
    70.         void ifft( const Vector< complex<Type> > &Xk, Vector<Type> &xn );

    71.         private:

    72.         void bitReverse( Vector<int> &bitrev );
    73.                 void radix2( int nthpo, Complex<Type> *c0, Complex<Type> *c1 );
    74.                 void radix4( int nthpo, Complex<Type> *c0, Complex<Type> *c1,
    75.                            Complex<Type> *c2, Complex<Type> *c3 );
    76.                 void radix8( int nxtlt, int nthpo, int length,
    77.                      Complex<Type> *cc0, Complex<Type> *cc1,
    78.                      Complex<Type> *cc2, Complex<Type> *cc3,
    79.                      Complex<Type> *cc4, Complex<Type> *cc5,
    80.                      Complex<Type> *cc6, Complex<Type> *cc7 );
    81.                 void dft( int direction, int n, Complex<Type> *b );
    82.         };
    83.         //        class FFTMR


    84.         #include <fftmr-impl.h>

    85. }
    86. // namespace splab


    87. #endif
    88. // FFTMR_H
    复制代码
    实现文件:
    1. /*
    2. * Copyright (c) 2008-2011 Zhang Ming (M. Zhang), zmjerry@163.com
    3. *
    4. * This program is free software; you can redistribute it and/or modify it
    5. * under the terms of the GNU General Public License as published by the
    6. * Free Software Foundation, either version 2 or any later version.
    7. *
    8. * Redistribution and use in source and binary forms, with or without
    9. * modification, are permitted provided that the following conditions are met:
    10. *
    11. * 1. Redistributions of source code must retain the above copyright notice,
    12. *    this list of conditions and the following disclaimer.
    13. *
    14. * 2. Redistributions in binary form must reproduce the above copyright
    15. *    notice, this list of conditions and the following disclaimer in the
    16. *    documentation and/or other materials provided with the distribution.
    17. *
    18. * This program is distributed in the hope that it will be useful, but WITHOUT
    19. * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    20. * FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for
    21. * more details. A copy of the GNU General Public License is available at:
    22. * http://www.fsf.org/licensing/licenses
    23. */


    24. /*****************************************************************************
    25. *                               fftmr-impl.h
    26. *
    27. * Implementation for FFTMR class.
    28. *
    29. * Zhang Ming, 2010-04, Xi'an Jiaotong University.
    30. *****************************************************************************/


    31. /**
    32. * To determine the input integer is or isn't the power of positive
    33. * integer of 2. If it is, return true, or return false.
    34. */
    35. bool isPower2( int n )
    36. {
    37.     int bits = 0;
    38.     while( n )
    39.     {
    40.         bits += n & 1;
    41.         n >>= 1;
    42.     }

    43.     return ( bits == 1 );
    44. }


    45. /**
    46. * Get binary log of integer argument -- exact if n is a power of 2.
    47. */
    48. int fastLog2( int n )
    49. {
    50.     int log = -1;
    51.     while( n )
    52.     {
    53.         log++;
    54.         n >>= 1;
    55.     }

    56.     return log;
    57. }


    58. /**
    59. * constructors and destructor
    60. */
    61. template<typename Type>
    62. FFTMR<Type>::FFTMR()
    63. {
    64. }

    65. template<typename Type>
    66. FFTMR<Type>::~FFTMR()
    67. {
    68. }


    69. /**
    70. * Generate the bit reversed array.
    71. */
    72. template<typename Type>
    73. void FFTMR<Type>::bitReverse( Vector<int> &bitrev )
    74. {
    75.     int i, j;
    76.         bitrev[0] = 0;

    77.         for( j=1; j<bitrev.size(); j<<=1 )
    78.                 for( i=0; i<j; ++i )
    79.                 {
    80.                         bitrev[i] <<= 1;
    81.                         bitrev[i+j] = bitrev[i]+1;
    82.                 }
    83. }


    84. /**
    85. * Radix-2 iteration subroutine.
    86. */
    87. template<typename Type>
    88. void FFTMR<Type>::radix2( int nthpo, Complex<Type> *c0, Complex<Type> *c1 )
    89. {
    90.     int     k, kk;
    91.     Type    *cr0, *ci0, *cr1, *ci1, r1, fi1;

    92.     cr0 = &(c0[0].re);
    93.     ci0 = &(c0[0].im);
    94.     cr1 = &(c1[0].re);
    95.     ci1 = &(c1[0].im);

    96.     for( k=0; k<nthpo; k+=2 )
    97.     {
    98.         kk = k*2;

    99.         r1 = cr0[kk] + cr1[kk];
    100.         cr1[kk] = cr0[kk] - cr1[kk];
    101.         cr0[kk] = r1;
    102.         fi1 = ci0[kk] + ci1[kk];
    103.         ci1[kk] = ci0[kk] - ci1[kk];
    104.         ci0[kk] = fi1;
    105.     }
    106. }


    107. /**
    108. * Radix-4 iteration subroutine.
    109. */
    110. template<typename Type>
    111. void FFTMR<Type>::radix4( int nthpo, Complex<Type> *c0, Complex<Type> *c1,
    112.                           Complex<Type> *c2, Complex<Type> *c3 )
    113. {
    114.     int     k, kk;
    115.     Type    *cr0, *ci0, *cr1, *ci1, *cr2, *ci2, *cr3, *ci3;
    116.     Type    r1, r2, r3, r4, i1, i2, i3, i4;

    117.     cr0 = &(c0[0].re);
    118.     cr1 = &(c1[0].re);
    119.     cr2 = &(c2[0].re);
    120.     cr3 = &(c3[0].re);
    121.     ci0 = &(c0[0].im);
    122.     ci1 = &(c1[0].im);
    123.     ci2 = &(c2[0].im);
    124.     ci3 = &(c3[0].im);

    125.     for( k=1; k<=nthpo; k+=4 )
    126.     {
    127.         kk = (k-1)*2;

    128.         r1 = cr0[kk] + cr2[kk];
    129.         r2 = cr0[kk] - cr2[kk];
    130.         r3 = cr1[kk] + cr3[kk];
    131.         r4 = cr1[kk] - cr3[kk];
    132.         i1 = ci0[kk] + ci2[kk];
    133.         i2 = ci0[kk] - ci2[kk];
    134.         i3 = ci1[kk] + ci3[kk];
    135.         i4 = ci1[kk] - ci3[kk];
    136.         cr0[kk] = r1 + r3;
    137.         ci0[kk] = i1 + i3;
    138.         cr1[kk] = r1 - r3;
    139.         ci1[kk] = i1 - i3;
    140.         cr2[kk] = r2 - i4;
    141.         ci2[kk] = i2 + r4;
    142.         cr3[kk] = r2 + i4;
    143.         ci3[kk] = i2 - r4;
    144.     }
    145. }


    146. /**
    147. * Radix-8 iteration subroutine.
    148. */
    149. template<typename Type>
    150. void FFTMR<Type>::radix8( int nxtlt, int nthpo, int length,
    151.                           Complex<Type> *cc0, Complex<Type> *cc1,
    152.                           Complex<Type> *cc2, Complex<Type> *cc3,
    153.                           Complex<Type> *cc4, Complex<Type> *cc5,
    154.                           Complex<Type> *cc6, Complex<Type> *cc7 )
    155. {
    156.     int     j, k, kk;
    157.     Type    scale, arg, tr, ti;

    158.     Type    c1, c2, c3, c4, c5, c6, c7,
    159.             s1, s2, s3, s4, s5, s6, s7;
    160.     Type    ar0, ar1, ar2, ar3, ar4, ar5, ar6, ar7,
    161.             ai0, ai1, ai2, ai3, ai4, ai5, ai6, ai7;
    162.     Type    br0, br1, br2, br3, br4, br5, br6, br7,
    163.             bi0, bi1, bi2, bi3, bi4, bi5, bi6, bi7;

    164.     Type    *cr0, *cr1, *cr2, *cr3, *cr4, *cr5, *cr6, *cr7;
    165.     Type    *ci0, *ci1, *ci2, *ci3, *ci4, *ci5, *ci6, *ci7;

    166.     cr0 = &(cc0[0].re);
    167.     cr1 = &(cc1[0].re);
    168.     cr2 = &(cc2[0].re);
    169.     cr3 = &(cc3[0].re);
    170.     cr4 = &(cc4[0].re);
    171.     cr5 = &(cc5[0].re);
    172.     cr6 = &(cc6[0].re);
    173.     cr7 = &(cc7[0].re);

    174.     ci0 = &(cc0[0].im);
    175.     ci1 = &(cc1[0].im);
    176.     ci2 = &(cc2[0].im);
    177.     ci3 = &(cc3[0].im);
    178.     ci4 = &(cc4[0].im);
    179.     ci5 = &(cc5[0].im);
    180.     ci6 = &(cc6[0].im);
    181.     ci7 = &(cc7[0].im);

    182.     scale = Type(TWOPI/length);

    183.     for( j=1; j<=nxtlt; j++ )
    184.     {
    185.         arg = (j-1)*scale;
    186.         c1 = cos(arg);
    187.         s1 = sin(arg);
    188.         c2 = c1*c1 - s1*s1;
    189.         s2 = c1*s1 + c1*s1;
    190.         c3 = c1*c2 - s1*s2;
    191.         s3 = c2*s1 + s2*c1;
    192.         c4 = c2*c2 - s2*s2;
    193.         s4 = c2*s2 + c2*s2;
    194.         c5 = c2*c3 - s2*s3;
    195.         s5 = c3*s2 + s3*c2;
    196.         c6 = c3*c3 - s3*s3;
    197.         s6 = c3*s3 + c3*s3;
    198.         c7 = c3*c4 - s3*s4;
    199.         s7 = c4*s3 + s4*c3;

    200.         for( k=j; k<=nthpo; k+=length )
    201.               {
    202.             kk = (k-1)*2; /* index by twos; re & im alternate */

    203.             ar0 = cr0[kk] + cr4[kk];
    204.             ar1 = cr1[kk] + cr5[kk];
    205.             ar2 = cr2[kk] + cr6[kk];
    206.             ar3 = cr3[kk] + cr7[kk];
    207.             ar4 = cr0[kk] - cr4[kk];
    208.             ar5 = cr1[kk] - cr5[kk];
    209.             ar6 = cr2[kk] - cr6[kk];
    210.             ar7 = cr3[kk] - cr7[kk];
    211.             ai0 = ci0[kk] + ci4[kk];
    212.             ai1 = ci1[kk] + ci5[kk];
    213.             ai2 = ci2[kk] + ci6[kk];
    214.             ai3 = ci3[kk] + ci7[kk];
    215.             ai4 = ci0[kk] - ci4[kk];
    216.             ai5 = ci1[kk] - ci5[kk];
    217.             ai6 = ci2[kk] - ci6[kk];
    218.             ai7 = ci3[kk] - ci7[kk];
    219.             br0 = ar0 + ar2;
    220.             br1 = ar1 + ar3;
    221.             br2 = ar0 - ar2;
    222.             br3 = ar1 - ar3;
    223.             br4 = ar4 - ai6;
    224.             br5 = ar5 - ai7;
    225.             br6 = ar4 + ai6;
    226.             br7 = ar5 + ai7;
    227.             bi0 = ai0 + ai2;
    228.             bi1 = ai1 + ai3;
    229.             bi2 = ai0 - ai2;
    230.             bi3 = ai1 - ai3;
    231.             bi4 = ai4 + ar6;
    232.             bi5 = ai5 + ar7;
    233.             bi6 = ai4 - ar6;
    234.             bi7 = ai5 - ar7;
    235.             cr0[kk] = br0 + br1;
    236.             ci0[kk] = bi0 + bi1;

    237.             if( j > 1 )
    238.                   {
    239.                 cr1[kk] = c4*(br0-br1) - s4*(bi0-bi1);
    240.                 cr2[kk] = c2*(br2-bi3) - s2*(bi2+br3);
    241.                 cr3[kk] = c6*(br2+bi3) - s6*(bi2-br3);
    242.                 ci1[kk] = c4*(bi0-bi1) + s4*(br0-br1);
    243.                 ci2[kk] = c2*(bi2+br3) + s2*(br2-bi3);
    244.                 ci3[kk] = c6*(bi2-br3) + s6*(br2+bi3);
    245.                 tr = Type(IRT2)*(br5-bi5);
    246.                 ti = Type(IRT2)*(br5+bi5);
    247.                 cr4[kk] = c1*(br4+tr) - s1*(bi4+ti);
    248.                 ci4[kk] = c1*(bi4+ti) + s1*(br4+tr);
    249.                 cr5[kk] = c5*(br4-tr) - s5*(bi4-ti);
    250.                 ci5[kk] = c5*(bi4-ti) + s5*(br4-tr);
    251.                 tr = -Type(IRT2)*(br7+bi7);
    252.                 ti = Type(IRT2)*(br7-bi7);
    253.                 cr6[kk] = c3*(br6+tr) - s3*(bi6+ti);
    254.                 ci6[kk] = c3*(bi6+ti) + s3*(br6+tr);
    255.                 cr7[kk] = c7*(br6-tr) - s7*(bi6-ti);
    256.                 ci7[kk] = c7*(bi6-ti) + s7*(br6-tr);
    257.                   }
    258.             else
    259.                   {
    260.                 cr1[kk] = br0 - br1;
    261.                 cr2[kk] = br2 - bi3;
    262.                 cr3[kk] = br2 + bi3;
    263.                 ci1[kk] = bi0 - bi1;
    264.                 ci2[kk] = bi2 + br3;
    265.                 ci3[kk] = bi2 - br3;
    266.                 tr = Type(IRT2)*(br5-bi5);
    267.                 ti = Type(IRT2)*(br5+bi5);
    268.                 cr4[kk] = br4 + tr;
    269.                 ci4[kk] = bi4 + ti;
    270.                 cr5[kk] = br4 - tr;
    271.                 ci5[kk] = bi4 - ti;
    272.                 tr = -Type(IRT2)*(br7+bi7);
    273.                 ti = Type(IRT2)*(br7-bi7);
    274.                 cr6[kk] = br6 + tr;
    275.                 ci6[kk] = bi6 + ti;
    276.                 cr7[kk] = br6 - tr;
    277.                 ci7[kk] = bi6 - ti;
    278.                   }
    279.               }
    280.     }
    281. }


    282. /**
    283. * This routine replaces the input Complex<Type> vector by its
    284. * finite discrete Fourier transform if direction == FORWARD.
    285. * It replaces the input Complex<Type> vector by its finite discrete
    286. * inverse Fourier transform if direction == INVERSE.
    287. */
    288. template<typename Type>
    289. void FFTMR<Type>::dft( int direction, int n, Complex<Type> *b )
    290. {
    291.     int     i, j;
    292.     int     n2pow, n8pow, nthpo, ipass, nxtlt, length;
    293.     Complex<Type> tmp;

    294.     n2pow = fastLog2(n);
    295.     nthpo = n;

    296.     // Conjugate the input
    297.     if( direction == FORWARD )
    298.         for( i=0; i<n; i++ )
    299.             b[i].im = -b[i].im;

    300.     n8pow = n2pow/3;
    301.     if( n8pow )
    302.     {
    303.         // radix 8 iterations
    304.         for( ipass=1; ipass<=n8pow; ipass++ )
    305.               {
    306.             nxtlt = 0x1 << ( n2pow - 3*ipass );
    307.             length = 8*nxtlt;
    308.             radix8( nxtlt, nthpo, length,
    309.                     b, b+nxtlt, b+2*nxtlt, b+3*nxtlt,
    310.                     b+4*nxtlt, b+5*nxtlt, b+6*nxtlt, b+7*nxtlt );
    311.               }
    312.     }

    313.     // A final radix 2 or radix 4 iteration is needed.
    314.     if( n2pow%3 == 1 )
    315.             radix2( nthpo, b, b+1 );
    316.     if( n2pow%3 == 2 )
    317.         radix4( nthpo, b, b+1, b+2, b+3 );

    318.      // scale outputs
    319.     if( direction == FORWARD )
    320.         for( i=0; i<n; i++ )
    321.             b[i].im *= -1;

    322.     if( direction == INVERSE )
    323.     {
    324.         for( i=0; i<nthpo; i++ )
    325.         {
    326.             b[i].re /= n;
    327.             b[i].im /= n;
    328.         }
    329.     }

    330.     // bit reverse
    331.     Vector<int> bitrev(n);
    332.     bitReverse( bitrev );
    333.         for( i=0; i<n; ++i )
    334.         {
    335.                 j = bitrev[i];
    336.                 if( i < j )
    337.                 {
    338.             tmp = b[i];
    339.             b[i] = b[j];
    340.             b[j] = tmp;
    341.                 }
    342.         }
    343. }


    344. /**
    345. * Forward Transform
    346. */
    347. template<typename Type>
    348. void FFTMR<Type>::fft( Vector<complex<Type> > &xn )
    349. {
    350.     int nFFT = xn.size();
    351.     if( isPower2(nFFT)  )
    352.         dft( FORWARD, nFFT, reinterpret_cast<Complex<Type>*>(xn.begin()) );
    353.     else
    354.         cerr << "The length of signal must abe power of 2!" ;
    355. }

    356. template<typename Type>
    357. void FFTMR<Type>::fft( const Vector<Type> &xn, Vector< complex<Type> > &Xk )
    358. {
    359.     int nFFT = xn.size();
    360.     if( isPower2(nFFT)  )
    361.     {
    362.         if( Xk.size() != nFFT )
    363.             Xk.resize(nFFT);
    364.         for( int i=0; i<nFFT; ++i )
    365.             Xk[i] = xn[i];

    366.         dft( FORWARD, nFFT, reinterpret_cast<Complex<Type>*>(Xk.begin()) );
    367.     }
    368.     else
    369.         cerr << "The length of signal must abe power of 2!" ;
    370. }


    371. /**
    372. * Inverse Transform
    373. */
    374. template<typename Type>
    375. void FFTMR<Type>::ifft( Vector<complex<Type> > &Xk )
    376. {
    377.     int nFFT = Xk.size();
    378.     if( isPower2(nFFT)  )
    379.         dft( INVERSE, nFFT, reinterpret_cast<Complex<Type>*>(Xk.begin()) );
    380.     else
    381.         cerr << "The length of signal must abe power of 2!" ;
    382. }

    383. template<typename Type>
    384. void FFTMR<Type>::ifft( const Vector< complex<Type> > &Xk, Vector<Type> &xn )
    385. {
    386.     int nFFT = Xk.size();
    387.     Vector< complex<Type> > tmp(Xk);
    388.     if( isPower2(nFFT)  )
    389.     {
    390.         dft( INVERSE, nFFT, reinterpret_cast<Complex<Type>*>(tmp.begin()) );

    391.         if( xn.size() != nFFT )
    392.             xn.resize(nFFT);
    393.         for( int i=0; i<nFFT; ++i )
    394.             xn[i] = tmp[i].real();
    395.     }
    396.     else
    397.         cerr << "The length of signal must abe power of 2!" ;
    398. }
    复制代码
    测试代码:
    1. /*****************************************************************************
    2. *                               fftmr_test.cpp
    3. *
    4. * Mixed Radix Algorithm FFT testing.
    5. *
    6. * Zhang Ming, 2010-04, Xi'an Jiaotong University.
    7. *****************************************************************************/


    8. #include <iostream>
    9. #include <iomanip>
    10. #include <fftmr.h>


    11. using namespace std;
    12. using namespace splab;


    13. typedef double  Type;
    14. const   int     LENGTH = 32;


    15. int main()
    16. {
    17.     int     i, j, index, rows = LENGTH/4;

    18.     Vector<Type> xn(LENGTH);
    19.     Vector< complex<Type> > yn(LENGTH),
    20.                               Xk(LENGTH);
    21.     FFTMR<Type> Fourier;

    22.     cout << "The original signal is: " << endl;
    23.     for( i=0; i<rows; i++ )
    24.     {
    25.         cout << endl;
    26.         for( j=0; j<3; j++ )
    27.         {
    28.             index = 3*i+j;
    29.             xn[index] = i+j;
    30.             cout << setiosflags(ios::fixed) << setprecision(6);
    31.             cout << "\t" << xn[index];
    32.         }
    33.     }
    34.     cout << endl << endl;

    35.     Fourier.fft( xn, Xk );

    36.     cout << "The Fourier transform of original signal is:" << endl;
    37.     for( i=0; i<rows; i++ )
    38.     {
    39.         cout << endl;
    40.         for( j=0; j<3; j++ )
    41.         {
    42.             index = 3*i+j;
    43.             cout << setiosflags(ios::fixed) << setprecision(6);
    44.             cout << "\t" << Xk[index];
    45.         }
    46.     }
    47.     cout << endl << endl;

    48.     Fourier.ifft( Xk, xn );
    49.     cout << "The inverse Fourier transform is" << endl;
    50.     for( i=0; i<rows; i++ )
    51.     {
    52.         cout << endl;
    53.         for( j=0; j<3; j++ )
    54.         {
    55.             index = 3*i+j;
    56.             cout << setiosflags(ios::fixed) << setprecision(6);
    57.             cout << "\t" << xn[index];
    58.         }
    59.     }
    60.     cout << endl << endl;

    61.     cout << "The original signal is: " << endl;
    62.     for( i=0; i<rows; i++ )
    63.     {
    64.         cout << endl;
    65.         for( j=0; j<3; j++ )
    66.         {
    67.             index = 3*i+j;
    68.             yn[index] = complex<double>(i,j);
    69.             cout << setiosflags(ios::fixed) << setprecision(6);
    70.             cout << "\t" << yn[index];
    71.         }
    72.     }
    73.     cout << endl << endl;

    74.     Fourier.fft( yn );
    75.     cout << "The Fourier transform of original signal is:" << endl;
    76.     for( i=0; i<rows; i++ )
    77.     {
    78.         cout << endl;
    79.         for( j=0; j<3; j++ )
    80.         {
    81.             index = 3*i+j;
    82.             cout << setiosflags(ios::fixed) << setprecision(6);
    83.             cout << "\t" << yn[index];
    84.         }
    85.     }
    86.     cout << endl << endl;

    87.     Fourier.ifft( yn );
    88.     cout << "The inverse Fourier transform is" << endl;
    89.     for( i=0; i<rows; i++ )
    90.     {
    91.         cout << endl;
    92.         for( j=0; j<3; j++ )
    93.         {
    94.             index = 3*i+j;
    95.             cout << setiosflags(ios::fixed) << setprecision(6);
    96.             cout << "\t" << yn[index];
    97.         }
    98.     }

    99.     cout << endl << endl;
    100.     return 0;
    101. }
    复制代码
    运行结果:
    1. The original signal is:

    2.         0.000000        1.000000        2.000000
    3.         1.000000        2.000000        3.000000
    4.         2.000000        3.000000        4.000000
    5.         3.000000        4.000000        5.000000
    6.         4.000000        5.000000        6.000000
    7.         5.000000        6.000000        7.000000
    8.         6.000000        7.000000        8.000000
    9.         7.000000        8.000000        9.000000

    10. The Fourier transform of original signal is:

    11.         (108.000000,0.000000)   (-52.827446,1.358153)   (-0.613126,-23.640092)
    12.         (13.310056,1.647944)    (-4.000000,9.656854)    (-8.902567,-4.353284)
    13.         (3.082392,-7.681941)    (5.303682,2.676653)     (-4.000000,4.000000)
    14.         (-5.314558,-4.659704)   (0.917608,-8.368233)    (0.220781,11.701395)
    15.         (-4.000000,1.656854)    (-0.343999,-3.954231)   (4.613126,-0.326383)
    16.         (0.554051,4.364941)     (-4.000000,-0.000000)   (0.554051,-4.364941)
    17.         (4.613126,0.326383)     (-0.343999,3.954231)    (-4.000000,-1.656854)
    18.         (0.220781,-11.701395)   (0.917608,8.368233)     (-5.314558,4.659704)

    19. The inverse Fourier transform is

    20.         0.000000        1.000000        2.000000
    21.         1.000000        2.000000        3.000000
    22.         2.000000        3.000000        4.000000
    23.         3.000000        4.000000        5.000000
    24.         4.000000        5.000000        6.000000
    25.         5.000000        6.000000        7.000000
    26.         6.000000        7.000000        8.000000
    27.         7.000000        8.000000        9.000000

    28. The original signal is:

    29.         (0.000000,0.000000)     (0.000000,1.000000)     (0.000000,2.000000)
    30.         (1.000000,0.000000)     (1.000000,1.000000)     (1.000000,2.000000)
    31.         (2.000000,0.000000)     (2.000000,1.000000)     (2.000000,2.000000)
    32.         (3.000000,0.000000)     (3.000000,1.000000)     (3.000000,2.000000)
    33.         (4.000000,0.000000)     (4.000000,1.000000)     (4.000000,2.000000)
    34.         (5.000000,0.000000)     (5.000000,1.000000)     (5.000000,2.000000)
    35.         (6.000000,0.000000)     (6.000000,1.000000)     (6.000000,2.000000)
    36.         (7.000000,0.000000)     (7.000000,1.000000)     (7.000000,2.000000)

    37. The Fourier transform of original signal is:

    38.         (84.000000,24.000000)   (-42.542528,1.020494)   (5.034128,-18.695144)
    39.         (13.685581,5.361763)    (-4.000000,9.656854)    (-6.244002,-4.826961)
    40.         (6.192124,-5.705118)    (6.023013,5.306175)     (-4.000000,4.000000)
    41.         (-1.277237,-6.299534)   (13.722858,-4.086928)   (-8.776842,-2.726593)
    42.         (-4.000000,1.656854)    (-1.716734,-2.445797)   (2.729646,0.149298)
    43.         (-0.485914,4.057658)    (-4.000000,-0.000000)   (0.861334,-3.324976)
    44.         (4.137445,2.209863)     (-1.852433,5.326967)    (-4.000000,-1.656854)
    45.         (14.648770,-2.703773)   (-3.363697,-4.437017)   (-3.674728,0.622383)

    46. The inverse Fourier transform is

    47.         (-0.000000,0.000000)    (0.000000,1.000000)     (0.000000,2.000000)
    48.         (1.000000,-0.000000)    (1.000000,1.000000)     (1.000000,2.000000)
    49.         (2.000000,0.000000)     (2.000000,1.000000)     (2.000000,2.000000)
    50.         (3.000000,-0.000000)    (3.000000,1.000000)     (3.000000,2.000000)
    51.         (4.000000,-0.000000)    (4.000000,1.000000)     (4.000000,2.000000)
    52.         (5.000000,-0.000000)    (5.000000,1.000000)     (5.000000,2.000000)
    53.         (6.000000,-0.000000)    (6.000000,1.000000)     (6.000000,2.000000)
    54.         (7.000000,0.000000)     (7.000000,1.000000)     (7.000000,2.000000)


    55. Process returned 0 (0x0)   execution time : 0.109 s
    56. Press any key to continue.
    复制代码



    本文来自:http://my.oschina.net/zmjerry/blog/3649

    守望者AIR技术交流社区(www.airmyth.com)
    回复 支持 反对

    使用道具 举报

  • TA的每日心情
    擦汗
    11 小时前
  • 签到天数: 433 天

    [LV.9]以坛为家II

    1739

    主题

    2091

    帖子

    12万

    积分

    超级版主

    Rank: 18Rank: 18Rank: 18Rank: 18Rank: 18

    威望
    457
    贡献
    24
    金币
    50924
    钢镚
    1419

    开源英雄守望者

     楼主| 发表于 2015-9-14 10:47:17 | 显示全部楼层
    Java

    1. /**Calculates the DFT using a split-radix Fast Fourier Transform algorithm.
    2. **/

    3. package sinica.music.tag.mfcc;

    4. public final class FFT {
    5.   public int logm;
    6.   final  int  MAXLOGM=20;     /* max FFT length      2^MAXLOGM */
    7.   final  double TWOPI=6.28318530717958647692;
    8.   final  double SQHALF=0.707106781186547524401;
    9.   int brseed[]= new int[4048];
    10.   float tab[][];

    11.   public FFT(int nlength ) {
    12.         double dtemp = Math.log(nlength) / Math.log(2);
    13.         if ( (dtemp - (int) dtemp) != 0.0) {
    14.             throw new Error("FFT length must be a power of 2.");
    15.         } else {
    16.             this.logm = (int) dtemp;
    17.         }
    18.         if (logm>=4) {
    19.             creattab(logm);
    20.         }
    21.   }

    22. /** Calculates the magnitude spectrum of a real signal.
    23.   * The returned vector contains only the positive frequencies.
    24.   */
    25.   public float[] calculateFFTMagnitude(float x[]) {
    26.         int i,n;
    27.     n=1<<this.logm;

    28.     if (x.length > n) {
    29.         throw new Error("Tried to use a " + n + "-points FFT for a vector with " +
    30.         x.length + " samples!");
    31.     }




    32.         rsfft(x);

    33.         float[] mag = new float[n/2 + 1];
    34.         mag[0] = x[0]; //DC frequency must be positive always

    35.         if (n==1) {
    36.             return mag;
    37.         }
    38.         mag[n/2] = (float) Math.abs(x[n/2]); //pi (meaning: fs / 2)

    39.         //System.out.println("FFT before magnitude");
    40.         //IO.DisplayVector(x);

    41.         for (i=1;i<n/2;i++) {
    42.             mag[i] = (float) Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
    43.             //System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
    44.         }

    45.         //IO.DisplayVector(mag);
    46.         return mag;
    47.   }

    48. /** Calculates the magnitude spectrum of a real signal.
    49.   * The returned vector contains only the positive frequencies.
    50.   */
    51.   public double[] calculateFFTMagnitude(double inputData[]) {
    52.         int i,n;
    53.         n=1<<this.logm;
    54.         if (inputData.length > n) {
    55.             throw new Error("Tried to use a " + n + "-points FFT for a vector with " +
    56.             inputData.length + " samples!");
    57.         }

    58.         //System.out.println("magnitude test");
    59.         //double[] dtest = DSP.DFTMagnitude(inputData);
    60.         //IO.DisplayVector(dtest);

    61.         float[] x = new float[n];
    62.         for (i=0; i<inputData.length; i++) {
    63.             x[i] = (float) inputData[i];
    64.         }

    65.         rsfft(x);

    66.         //System.out.println("FFT before magnitude");
    67.         //IO.DisplayVector(x);

    68.         double[] mag = new double[n/2 + 1];
    69.         mag[0] = x[0]; //DC frequency must be positive always

    70.         if (n==1) {
    71.             return mag;
    72.         }
    73.         mag[n/2] = Math.abs(x[n/2]); //pi (meaning: fs / 2)

    74.         for (i=1;i<n/2;i++) {
    75.             mag[i] = Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
    76.             //System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
    77.         }

    78.         //IO.DisplayVector(mag);
    79.         return mag;
    80.   }

    81. /** Calculates the power (magnitude squared) spectrum of a real signal.
    82.   * The returned vector contains only the positive frequencies.
    83.   */
    84.   public double[] calculateFFTPower(double inputData[]) {
    85.         int i,n;
    86.         n=1<<this.logm;

    87.         //System.out.println("magnitude test");
    88.         //double[] dtest = DSP.DFTMagnitude(inputData);
    89.         //IO.DisplayVector(dtest);

    90.         float[] x = new float[n];
    91.         for (i=0; i<inputData.length; i++) {
    92.             x[i] = (float) inputData[i];
    93.         }

    94.         rsfft(x);

    95.         //System.out.println("FFT before magnitude");
    96.         //IO.DisplayVector(x);

    97.         double[] mag = new double[n/2 + 1];
    98.         mag[0] = x[0]; //DC frequency must be positive always

    99.         if (n==1) {
    100.             return mag;
    101.         }
    102.         mag[n/2] = Math.abs(x[n/2]); //pi (meaning: fs / 2)

    103.         for (i=1;i<n/2;i++) {
    104.             mag[i] = x[i]*x[i]+x[n-i]*x[n-i];
    105.             //mag[i] = Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
    106.             //System.out.println(mag[i] + " " + x[i] + " " + x[n-i]);
    107.         }

    108.         //IO.DisplayVector(mag);
    109.         return mag;
    110.   }

    111.   /**In place calculation of FFT magnitude.
    112.    */
    113. public void FFTMagnitude(float x[])
    114. { int i,n;
    115.   rsfft(x);
    116.   n=1<<this.logm;
    117.   if (n==1) return;
    118.   for (i=1;i<n/2;i++)
    119.   {x[i]=(float)Math.sqrt(x[i]*x[i]+x[n-i]*x[n-i]);
    120.    x[n-i]=x[i];
    121.   }
    122.   //Aldebaro modification:
    123.     x[n/2] = Math.abs(x[n/2]);
    124. }

    125.   void rsfft(float x[])
    126.   {
    127.       /*creat table*/
    128.   //  if(logm>=4) creattab(logm);
    129.       /* Call recursive routine */

    130.       rsrec(x,logm);

    131.     /* Output array unshuffling using bit-reversed indices */
    132.     if (logm > 1) {
    133.         BR_permute(x, logm);
    134.         return ;
    135. }
    136.   }

    137. /* -------------------------------------------------------------------- *
    138.    *   Inverse  transform  for  real  inputs                              *
    139.    *--------------------------------------------------------------------  */

    140. void  rsifft(float x[])
    141. {
    142.    int       i, m;
    143.    float     fac;
    144.    int         xp;

    145.    /* Output array unshuffling using bit-reversed indices */
    146.    if (logm > 1) {
    147.       BR_permute(x, logm);
    148.    }
    149.    x[0] *= 0.5;
    150.    if (logm > 0) x[1] *= 0.5;
    151. /*creat table*/
    152.     //if(logm>=4) creattab(logm);
    153.    /*  Call  recursive  routine */

    154.    rsirec(x, logm);

    155.    /* Normalization */
    156.    m = 1 << logm;
    157.    fac = (float)2.0 / m;
    158.    xp = 0;

    159.    for (i = 0; i < m; i++) {
    160.        x[xp++] *= fac;
    161.    }
    162. }

    163. /* -------------------------------------------------------------------- *
    164. *     Creat multiple fator table                                       *
    165. * -------------------------------------------------------------------- */

    166. void   creattab(int logm)
    167. { int       m, m2, m4, m8, nel, n, rlogm;
    168.    int       cn, spcn, smcn, c3n, spc3n, smc3n;
    169.    double    ang, s, c;
    170.    tab=new float [logm-4+1][6*((1<<logm)/4-2)];
    171.   for(rlogm=logm; rlogm>=4;rlogm--)

    172.     {m = 1 << rlogm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2; nel=m4-2;
    173. /* Initialize pointers */

    174.          cn =0; spcn = cn + nel;  smcn = spcn + nel;c3n = smcn + nel; spc3n = c3n + nel; smc3n = spc3n + nel;


    175.     /* Compute tables */
    176.         for (n = 1; n < m4; n++) {
    177.            if (n == m8) continue;
    178.              ang = n * TWOPI / m;
    179.              c = Math.cos(ang);  s = Math.sin(ang);
    180.              tab[rlogm-4][cn++] = (float)c;  tab[rlogm-4][spcn++] = (float)(- (s + c)); tab[rlogm-4][smcn++] =(float)( s - c);

    181.              ang = 3 * n * TWOPI / m;
    182.              c = Math.cos(ang);  s = Math. sin(ang);
    183.              tab[rlogm-4][c3n++] = (float)c; tab[rlogm-4][spc3n++] = (float)(- (s + c)); tab[rlogm-4][smc3n++] = (float)(s - c);
    184.         }
    185.     }
    186. }

    187. /* -------------------------------------------------------------------- *
    188. *     Recursive part of the RSFFT algorithm.       Not externally      *
    189. *     callable.                                                        *
    190. * -------------------------------------------------------------------- */

    191.   void  rsrec(float x[],int logm)
    192. {
    193.      int       m, m2, m4, m8, nel, n;
    194.      int       x0=0;
    195.      int       xr1, xr2, xi1;
    196.      int       cn=0;
    197.      int       spcn=0;
    198.      int       smcn=0;
    199.      float     tmp1, tmp2;
    200.      double    ang, c, s;




    201.      /* Check range   of logm */
    202.      try{ if ((logm < 0) || (logm > MAXLOGM)) {
    203.            System.err.println("FFT length m is too big: log2(m) = "+logm+"is out of bounds ["+0+","+MAXLOGM+"]");

    204.         throw new OutofborderException(logm);
    205.           }}
    206.      catch( OutofborderException e)
    207.      {throw new OutOfMemoryError();}

    208.      /* Compute trivial cases */

    209.      if (logm < 2) {
    210.           if (logm == 1) {    /* length m = 2 */
    211.             xr2  = x0 + 1;
    212.             tmp1 = x[x0] + x[xr2];
    213.             x[xr2] = x[x0] - x[xr2];
    214.             x[x0]   =  tmp1;
    215.             return;
    216.           }
    217.         else if (logm == 0) return;      /* length m = 1 */
    218.     }

    219.      /* Compute a few constants */
    220.      m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;

    221.      /* Build tables of butterfly coefficients, if necessary */
    222.      //if ((logm >= 4) && (tab[logm-4][0] == 0)) {

    223.          /* Allocate memory for tables */
    224.        //  nel = m4 - 2;

    225.     /*if ((tab[logm-4] = (float *) calloc(3 * nel, sizeof(float)))
    226.        == NULL) {
    227.        printf("Error : RSFFT : not enough memory for cosine tables.\n");
    228.        error_exit();
    229.     }*/


    230.     /* Initialize pointers */
    231.     //tabi=logm-4;
    232.     //   cn  =0; spcn = cn + nel;  smcn = spcn + nel;

    233.     /* Compute tables */
    234.     /*for (n = 1; n  m4; n++) {
    235.        if (n == m8) continue;
    236.        ang = n * (float)TWOPI / m;
    237.        c = Math.cos(ang);  s = Math.sin(ang);
    238.        tab[tabi][cn++] = (float)c;  tab[tabi][spcn++] = (float)(- (s + c)); tab[tabi][smcn++] =(float)( s - c);
    239.    }
    240. }

    241. /*  Step  1 */
    242. xr1 = x0;  xr2 = xr1 + m2;
    243. for (n = 0; n < m2; n++) {
    244.     tmp1 = x[xr1] + x[xr2];
    245.     x[xr2] = x[xr1] - x[xr2];
    246.     x[xr1] = tmp1;
    247.     xr1++; xr2++;
    248. }

    249. /*  Step  2        */
    250. xr1 = x0 + m2 + m4;
    251. for (n = 0; n < m4; n++) {
    252.     x[xr1] = - x[xr1];
    253.     xr1++;
    254. }

    255. /*  Steps 3 &  4 */
    256. xr1 = x0 + m2; xi1 = xr1 + m4;
    257. if (logm >= 4) {
    258.     nel = m4 - 2;
    259.     cn  = 0; spcn = cn + nel;  smcn = spcn + nel;
    260. }

    261. xr1++; xi1++;
    262. for (n = 1; n < m4; n++) {
    263.     if (n == m8) {
    264.        tmp1 = (float)( SQHALF * (x[xr1] + x[xi1]));
    265.        x[xi1]  = (float)(SQHALF * (x[xi1] - x[xr1]));
    266.        x[xr1]  = tmp1;
    267.     }  else {//System.out.println ("logm-4="+(logm-4));
    268.        tmp2 = tab[logm-4][cn++] * (x[xr1] + x[xi1]);
    269.        tmp1 = tab[logm-4][spcn++] * x[xr1] + tmp2;
    270.        x[xr1] = tab[logm-4][smcn++] * x[xi1] + tmp2;
    271.        x[xi1] = tmp1;
    272.        }
    273.     //System.out.println ("logm-4="+(logm-4));
    274.     xr1++; xi1++;
    275. }

    276.     /*  Call rsrec again with half DFT length */
    277.     rsrec(x,logm-1);

    278.     /* Call complex DFT routine, with quarter DFT length.
    279.         Constants have to be recomputed, because they are static! */
    280.     m = 1 << logm; m2 = m / 2; m4 = 3 * (m / 4);
    281.     srrec(x,x0 + m2, x0 + m4, logm-2);

    282.     /* Step 5: sign change & data reordering */
    283.     m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;
    284.     xr1 = x0 + m2 + m4;
    285.     xr2 = x0 + m - 1;
    286.     for (n = 0; n < m8; n++) {
    287.        tmp1   = x[xr1];
    288.        x[xr1++] = - x[xr2];
    289.        x[xr2--] = - tmp1;
    290.     }
    291.     xr1 = x0 + m2 + 1;
    292.     xr2 = x0 + m - 2;
    293.     for (n = 0; n < m8; n++) {
    294.         tmp1   =   x[xr1];
    295.         x[xr1++] = - x[xr2];
    296.         x[xr2--] =   tmp1;
    297.         xr1++;
    298.         xr2--;
    299.     }
    300.     if (logm == 2) x[3] = -x[3];
    301. }
    302. /* --------------------------------------------------------------------- *
    303. *  Recursive part of the inverse RSFFT algorithm.  Not externally       *
    304. *  callable.                                                            *
    305. *  -------------------------------------------------------------------- */

    306. void  rsirec(float  x[],  int   logm)
    307. {
    308.      int       m, m2, m4, m8, nel, n;
    309.      int       xr1, xr2, xi1;
    310.      int       x0=0;
    311.      int       cn, spcn, smcn;
    312.      float     tmp1, tmp2;
    313.      cn=0;smcn=0;spcn=0;

    314.      /* Check  range  of logm */
    315.       try{ if ((logm < 0) || (logm > MAXLOGM)) {
    316.         System.err.println("FFT length m is too big: log2(m) = "+logm+"is out of bounds ["+0+","+MAXLOGM+"]");
    317.         throw new OutofborderException(logm);
    318.           }}
    319.      catch( OutofborderException e)
    320.      {throw new OutOfMemoryError();}

    321.      /*  Compute  trivial  cases */
    322.      if (logm < 2) {
    323.         if (logm == 1) {     /* length m = 2 */
    324.            xr2  = x0 + 1;
    325.            tmp1 = x[x0] + x[xr2];
    326.            x[xr2] = x[x0] - x[xr2];
    327.            x[0]= tmp1;
    328.            return;
    329.         }
    330.        else if (logm == 0) return;       /* length m = 1 */
    331.     }

    332.      /* Compute a few constants */
    333.      m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;

    334.      /* Build tables of butterfly    coefficients, if necessary */
    335.      // if((logm >= 4) && (tab[logm-4] == NULL)) {

    336.        /* Allocate memory for tables */
    337.       /*el = m4 - 2;
    338.        if ((tab[logm-4] = (float *) calloc(3 * nel, sizeof(float)))
    339.            == NULL) {
    340.            printf("Error : RSFFT : not enough memory for cosine tables.\n");
    341.            error_exit();
    342.        }

    343.        /*  Initialize   pointers */
    344.        //cn  = tab[logm-4] ; spcn = cn + nel; smcn = spcn + nel;

    345.      /*  Compute  tables */
    346.      /* (n = 1; n  m4; n++) {
    347.          if (n == m8) continue;
    348.          ang = n * TWOPI / m;
    349.          c = cos(ang); s = sin(ang);
    350.          *cn++ = c; *spcn++ = - (s + c); *smcn++ = s - c;
    351.     }
    352. }
    353. /* Reverse Step 5: sign change & data reordering */
    354. if (logm == 2) x[3] = -x[3];
    355. xr1 = x0+ m2 + 1;
    356. xr2 = x0+ m - 2;
    357. for (n = 0; n < m8; n++) {
    358.      tmp1   =   x[xr1];
    359.      x[xr1++] =   x[xr2];
    360.      x[xr2--] = - tmp1;
    361.      xr1++;
    362.      xr2--;
    363. }
    364. xr1 = x0 + m2 + m4;
    365. xr2 = x0 + m - 1;
    366. for (n = 0; n < m8; n++) {
    367.      tmp1   =   x[xr1];
    368.      x[xr1++] = - x[xr2];
    369.      x[xr2--] = - tmp1;
    370. }
    371. /*  Call   rsirec again with half DFT length */
    372. rsirec(x, logm-1);

    373. /* Call complex DFT routine, with quarter DFT length.
    374.      Constants have to be recomputed, because they are static! */

    375. /*Now in Java version, we set the multiple Constant to be global*/
    376. m = 1 << logm; m2 = m / 2; m4 = 3 * (m / 4);
    377. srrec(x,x0 + m4, x0 + m2, logm-2);

    378. /* Reverse Steps 3 & 4 */
    379. m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;
    380. xr1 = x0 + m2; xi1 = xr1 + m4;
    381. if (logm >= 4) {
    382.      nel = m4 - 2;
    383.      cn  = 0; spcn = cn + nel; smcn = spcn + nel;
    384. }
    385. xr1++; xi1++;
    386. for (n = 1; n < m4; n++) {
    387.      if (n == m8) {
    388.          tmp1 = (float)(SQHALF * (x[xr1] - x[xi1]));
    389.          x[xi1] = (float)(SQHALF * (x[xi1] + x[xr1]));
    390.          x[xr1] = tmp1;
    391.      }    else {
    392.             tmp2 = tab[logm-4][cn++] * (x[xr1] + x[xi1]);
    393.             tmp1 = tab[logm-4][smcn++] * x[xr1] + tmp2;
    394.             x[xr1] = tab[logm-4][spcn++] * x[xi1] + tmp2;
    395.             x[xi1] = tmp1;
    396.          }
    397.         xr1++; xi1++;
    398. }

    399. /* Reverse Step 2 */
    400. xr1 = x0 + m2 + m4;
    401. for (n = 0; n < m4; n++) {
    402.      x[xr1] = - x[xr1];
    403.      xr1++;
    404. }

    405. /* Reverse  Step  1 */
    406. xr1 = x0; xr2 = xr1 + m2;
    407. for (n = 0; n < m2; n++) {
    408.      tmp1 = x[xr1] + x[xr2];
    409.      x[xr2] = x[xr1] - x[xr2];
    410.      x[xr1] = tmp1;
    411.      xr1++; xr2++;
    412. }
    413. }


    414.   /* -------------------------------------------------------------------- *
    415. *      Recursive part of the SRFFT algorithm.                          *
    416. * -------------------------------------------------------------------- */

    417. void srrec(float x[],int xr, int xi, int logm)
    418. {
    419.    int        m, m2, m4, m8, nel, n;
    420.   // int        x0=0;
    421.    int        xr1, xr2, xi1, xi2;
    422.    int        cn, spcn, smcn, c3n, spc3n, smc3n;
    423.    float      tmp1, tmp2;
    424.   cn=0; spcn=0; smcn=0; c3n=0; spc3n=0; smc3n=0;




    425. /* Check range of logm */
    426.    try
    427.    {if ((logm < 0) || (logm > MAXLOGM)) {
    428.     System.err.println("FFT length m is too big: log2(m) = "+logm+"is out of bounds ["+0+","+MAXLOGM+"]");

    429.         throw new OutofborderException(logm) ;
    430.     }
    431.    }
    432.     catch ( OutofborderException e)
    433.     {throw new OutOfMemoryError();}

    434. /*  Compute trivial cases */
    435. if (logm < 3) {
    436.       if (logm == 2) {  /* length m = 4 */
    437.        xr2  = xr + 2;
    438.        xi2  = xi + 2;
    439.        tmp1 = x[xr] + x[xr2];
    440.        x[xr2] = x[xr] - x[xr2];
    441.        x[xr]  = tmp1;
    442.        tmp1 = x[xi] + x[xi2];
    443.        x[xi2] = x[xi] - x[xi2];
    444.        x[xi]  = tmp1;
    445.        xr1  = xr + 1;
    446.        xi1  = xi + 1;
    447.        xr2++;
    448.        xi2++;
    449.        tmp1 = x[xr1] + x[xr2];
    450.        x[xr2] = x[xr1] - x[xr2];
    451.        x[xr1] = tmp1;
    452.        tmp1 = x[xi1] + x[xi2];
    453.        x[xi2] = x[xi1] - x[xi2];
    454.        x[xi1] = tmp1;
    455.        xr2  = xr + 1;
    456.        xi2  = xi + 1;
    457.        tmp1 = x[xr] + x[xr2];
    458.        x[xr2] = x[xr] - x[xr2];
    459.        x[xr]  = tmp1;
    460.        tmp1 = x[xi] + x[xi2];
    461.        x[xi2] = x[xi] - x[xi2];
    462.        x[xi]  = tmp1;
    463.        xr1  = xr + 2;
    464.        xi1  = xi + 2;
    465.        xr2  = xr + 3;
    466.        xi2  = xi + 3;
    467.        tmp1 = x[xr1] + x[xi2];
    468.        tmp2 = x[xi1] + x[xr2];
    469.        x[xi1] = x[xi1] - x[xr2];
    470.        x[xr2] = x[xr1] - x[xi2];
    471.        x[xr1] =tmp1;
    472.        x[xi2] =tmp2;
    473.        return;
    474. }

    475.     else  if (logm == 1) { /* length m = 2 */
    476.        xr2   = xr +  1;
    477.        xi2   = xi +  1;
    478.        tmp1  = x[xr] + x[xr2];
    479.        x[xr2]  = x[xr] - x[xr2];
    480.        x[xr]   = tmp1;
    481.        tmp1  = x[xi] + x[xi2];
    482.        x[xi2]  = x[xi] - x[xi2];
    483.        x[xi]   = tmp1;
    484.        return;
    485.     }
    486.     else if (logm == 0) return;     /* length m = 1*/
    487. }

    488. /* Compute a few constants */
    489. m = 1 << logm; m2 = m / 2; m4 = m2 / 2; m8 = m4 / 2;

    490. /* Build tables of butterfly coefficients, if necessary */
    491. //if ((logm >= 4) && (tab[logm-4] == NULL)) {

    492.     /* Allocate memory for tables */
    493.   /*  nel = m4 - 2;
    494.     if ((tab[logm-4] = (float *) calloc(6 * nel, sizeof(float)))
    495.        == NULL) {
    496.        exit(1);
    497.      }
    498.     /* Initialize pointers */

    499.     /*cn  = tab[logm-4]; spcn = cn + nel;  smcn = spcn + nel;
    500.     c3n = smcn + nel; spc3n = c3n + nel; smc3n = spc3n + nel;


    501.     /* Compute tables */
    502.     /*for (n = 1; n  m4; n++) {
    503.        if (n == m8) continue;
    504.        ang = n * TWOPI / m;
    505.        c = cos(ang); s = sin(ang);
    506.        *cn++ = c; *spcn++ = - (s + c); *smcn++ = s - c;
    507.        ang = 3 * n * TWOPI / m;
    508.        c = cos(ang); s = sin(ang);
    509.        *c3n++ = c; *spc3n++ = - (s + c); *smc3n++ = s - c;
    510.    }
    511. }


    512. /*  Step 1 */
    513. xr1 = xr;  xr2 = xr1  +  m2;
    514. xi1 = xi;  xi2 = xi1  +  m2;

    515. for (n = 0; n < m2; n++) {
    516.    tmp1 = x[xr1] + x[xr2];
    517.    x[xr2] = x[xr1] - x[xr2];
    518.    x[xr1] = tmp1;
    519.    tmp2 = x[xi1] + x[xi2];
    520.    x[xi2] = x[xi1] - x[xi2];
    521.    x[xi1] = tmp2;
    522.    xr1++;  xr2++;  xi1++;  xi2++;
    523. }
    524. /*   Step 2  */
    525. xr1 = xr + m2; xr2 = xr1 + m4;
    526. xi1 = xi + m2; xi2 = xi1 + m4;
    527. for (n = 0; n < m4; n++) {
    528.    tmp1 = x[xr1] + x[xi2];
    529.    tmp2 = x[xi1] + x[xr2];
    530.    x[xi1] = x[xi1] - x[xr2];
    531.    x[xr2] = x[xr1] - x[xi2];
    532.    x[xr1] = tmp1;
    533.    x[xi2] = tmp2;
    534.    xr1++;  xr2++;  xi1++;  xi2++;
    535. }

    536. /*   Steps  3 & 4 */
    537. xr1 = xr + m2; xr2 = xr1 + m4;
    538. xi1 = xi + m2; xi2 = xi1 + m4;
    539. if (logm >= 4) {
    540.    nel = m4 - 2;
    541.    cn  = 0; spcn  = cn + nel;  smcn  = spcn + nel;
    542.    c3n = smcn + nel;  spc3n = c3n + nel; smc3n = spc3n + nel;
    543. }
    544. xr1++; xr2++; xi1++; xi2++;
    545. for (n = 1; n < m4; n++) {
    546.    if (n == m8) {
    547.        tmp1 = (float)(SQHALF * (x[xr1] + x[xi1]));
    548.        x[xi1] = (float)(SQHALF * (x[xi1] - x[xr1]));
    549.        x[xr1] = tmp1;
    550.        tmp2 = (float)(SQHALF * (x[xi2] - x[xr2]));
    551.        x[xi2] = (float)(-SQHALF * (x[xr2] + x[xi2]));
    552.        x[xr2] = tmp2;
    553. }     else {
    554.        tmp2 = tab[logm-4][cn++] * (x[xr1] + x[xi1]);
    555.        tmp1 = tab[logm-4][spcn++] * x[xr1] + tmp2;
    556.        x[xr1] = tab[logm-4][smcn++] * x[xi1] + tmp2;
    557.        x[xi1] = tmp1;
    558.        tmp2 = tab[logm-4][c3n++] * (x[xr2] + x[xi2]);
    559.        tmp1 = tab[logm-4][spc3n++] * x[xr2] + tmp2;
    560.        x[xr2] = tab[logm-4][smc3n++] * x[xi2] + tmp2;
    561.        x[xi2] = tmp1;
    562. }
    563. // System.out.println ("logm-4="+(logm-4));
    564.    xr1++; xr2++; xi1++; xi2++;
    565. }
    566.    /* Call ssrec again with half DFT length  */
    567.    srrec(x,xr, xi, logm-1);

    568.    /* Call ssrec again twice with one quarter DFT length.
    569.      Constants have to be recomputed, because they are static!*/
    570.    m = 1 << logm; m2 = m / 2;
    571.    srrec(x,xr + m2, xi + m2, logm-2);
    572.    m = 1 << logm; m4 = 3 * (m / 4);
    573.    srrec(x,xr + m4, xi + m4, logm-2);
    574. }

    575. /* -------------------------------------------------------------------- *
    576. *    Data unshuffling according to bit-reversed indexing.              *
    577. *                                                                      *
    578. *                                                                      *
    579. *    Bit reversal is done using Evan's algorithm (Ref: D. M. W.        *
    580. *    Evans, "An improved digit-reversal permutation algorithm...",     *
    581. *    IEEE Trans.  ASSP, Aug. 1987, pp. 1120-1125).                     *
    582. * -------------------------------------------------------------------- */

    583. //static    int   brseed[256];     /* Evans' seed table */
    584. //static    int     brsflg;         /* flag for table building */


    585. void creatbrseed( int logm)
    586. {int lg2,n;
    587. lg2 = logm >> 1;
    588. n =1 << lg2;
    589. if (logm!=(logm>>1)<<1) lg2++;
    590. brseed[0] = 0;
    591.        brseed[1] = 1;
    592.        for  (int j=2; j <= lg2; j++) {
    593.           int imax = 1 << (j-1);
    594.           for (int i = 0; i < imax; i++) {
    595.              brseed[i] <<= 1;
    596.              brseed[i + imax] = brseed[i] + 1;
    597.           }
    598.        }
    599. }
    600. void BR_permute(float x[], int logm)
    601. {
    602.    int       i, j, imax, lg2, n;
    603.    int       off, fj, gno;
    604.    float     tmp;
    605.    int       xp, xq, brp;
    606.    int       x0=0;

    607.    lg2 =  logm >> 1;
    608.    n = 1  << lg2;
    609.    if  (logm !=(logm>>1)<<1) lg2++;

    610.    /*  Create seed table if not yet built */
    611.   /* if  (brsflg != logm) {
    612.        brsflg = logm;
    613.        brseed[0] = 0;
    614.        brseed[1] = 1;
    615.        for  (j=2; j = lg2; j++) {
    616.           imax = 1 < (j-1);
    617.           for (i = 0; i  imax; i++) {
    618.              brseed[i] <= 1;
    619.              brseed[i + imax] = brseed[i] + 1;
    620.           }
    621.        }
    622.    }*/
    623.   creatbrseed(logm);

    624.     /*  Unshuffling   loop */
    625.     for (off = 1; off < n; off++) {
    626.          fj = n * brseed[off]; i = off; j = fj;
    627.          tmp = x[i]; x[i] = x[j]; x[j] = tmp;
    628.          xp = i;
    629.          brp = 1;

    630.          for (gno = 1; gno < brseed[off]; gno++) {
    631.              xp += n;
    632.              j = fj + brseed[brp++];
    633.              xq = x0 + j;
    634.              tmp = x[xp]; x[xp] = x[xq]; x[xq] = tmp;
    635.          }
    636.     }

    637.   }
    638.     class OutofborderException extends Exception {
    639.         public OutofborderException (int logm)
    640.         {super();
    641.         }
    642.     }
    643. }
    复制代码


    守望者AIR技术交流社区(www.airmyth.com)
    回复 支持 反对

    使用道具 举报

    您需要登录后才可以回帖 登录 | 立即注册

    本版积分规则

    
    关闭

    站长推荐上一条 /4 下一条

    QQ|手机版|Archiver|网站地图|小黑屋|守望者 ( 京ICP备14061876号

    GMT+8, 2017-10-19 22:30 , Processed in 1.203125 second(s), 32 queries .

    守望者AIR

    守望者AIR技术交流社区

    本站成立于 2014年12月31日

    快速回复 返回顶部 返回列表