FFT（快速傅立叶算法　for java）

た入场券 2022-09-19 00:16 190阅读 0赞

public class FFT \{     
    public static final int FFT\_N\_LOG = 10; // FFT\_N\_LOG <= 13     
    public static final int FFT\_N = 1 << FFT\_N\_LOG;     
    private static final float MINY = (float) ((FFT\_N << 2) \* Math.sqrt(2)); //(\*)     
    private float\[\] real, imag, sintable, costable;     
    private int\[\] bitReverse;

public FFT() \{

real = new float\[FFT\_N\];     
        imag = new float\[FFT\_N\];     
        sintable = new float\[FFT\_N >> 1\];     
        costable = new float\[FFT\_N >> 1\];     
        bitReverse = new int\[FFT\_N\];     
    
        int i, j, k, reve;     
        for (i = 0; i < FFT\_N; i++) \{     
            k = i;     
            for (j = 0, reve = 0; j != FFT\_N\_LOG; j++) \{     
                reve <<= 1;     
                reve |= (k & 1);     
                k >>>= 1;     
            \}     
            bitReverse\[i\] = reve;     
        \}     
    
        double theta, dt = 2 \* 3.14159265358979323846 / FFT\_N;     
        for (i = 0; i < (FFT\_N >> 1); i++) \{     
            theta = i \* dt;     
            costable\[i\] = (float) Math.cos(theta);     
            sintable\[i\] = (float) Math.sin(theta);     
        \}     
    \}     
    
    /\*\*    
     \* 用于频谱显示的快速傅里叶变换    
     \* @param realIO 输入FFT\_N个实数，也用它暂存fft后的FFT\_N/2个输出值(复数模的平方)。    
     \*/    
    public void calculate(float\[\] realIO) \{     
        int i, j, k, ir, exchanges = 1, idx = FFT\_N\_LOG - 1;     
        float cosv, sinv, tmpr, tmpi;     
        for (i = 0; i != FFT\_N; i++) \{     
            real\[i\] = realIO\[bitReverse\[i\]\];     
            imag\[i\] = 0;     
        \}     
    
        for (i = FFT\_N\_LOG; i != 0; i--) \{     
            for (j = 0; j != exchanges; j++) \{     
                cosv = costable\[j << idx\];     
                sinv = sintable\[j << idx\];     
                for (k = j; k < FFT\_N; k += exchanges << 1) \{     
                    ir = k + exchanges;     
                    tmpr = cosv \* real\[ir\] - sinv \* imag\[ir\];     
                    tmpi = cosv \* imag\[ir\] + sinv \* real\[ir\];     
                    real\[ir\] = real\[k\] - tmpr;     
                    imag\[ir\] = imag\[k\] - tmpi;     
                    real\[k\] += tmpr;     
                    imag\[k\] += tmpi;     
                \}     
            \}     
            exchanges <<= 1;     
            idx--;     
        \}     
    
        j = FFT\_N >> 1;     
        /\*    
         \* 输出模的平方(的FFT\_N倍):    
         \* for(i = 1; i <= j; i++)    
         \*  realIO\[i-1\] = real\[i\] \* real\[i\] +  imag\[i\] \* imag\[i\];    
         \*     
         \* 如果FFT只用于频谱显示,可以"淘汰"幅值较小的而减少浮点乘法运算. MINY的值    
         \* 和Spectrum.Y0,Spectrum.logY0对应.    
         \*/    
        sinv = MINY;     
        cosv = -MINY;     
        for (i = j; i != 0; i--) \{     
            tmpr = real\[i\];     
            tmpi = imag\[i\];     
            if (tmpr > cosv && tmpr < sinv && tmpi > cosv && tmpi < sinv)     
                realIO\[i - 1\] = 0;     
            else    
                realIO\[i - 1\] = tmpr \* tmpr + tmpi \* tmpi;     
        \}     
    \}     
    
\}