dofft.java

自己编写的

package fft;

/**
 * 此处插入类型描述。
 * 创建日期：(2003-9-27 20:01:02)
 * @author：Administrator
 */
import java.awt.*;
import java.io.*;

public class DoFFT {
    private double inFFT[][], outFFT[][];
    private double output_pow[];
    
    public double input[], output[];
    public int N = 64;


    public DoFFT(String strInputFile) {

        FileReader(strInputFile);

    } //end of constructor

    public void FileReader(String strFileName) {
        try {
            FileInputStream fis;
            DataInputStream dis;
            
            input = new double[N];

            //xha
            fis = new FileInputStream(strFileName);
            dis = new DataInputStream(fis);


            for (int i = 0; i < N; i++) {
                input[i] = dis.readFloat();
                //System.out.println(i+" "+input[i]);
            }

            dis.close();
            
        } catch (IOException e) {
	        System.err.println("File not opened properly\n"+e);
	        System.exit(1);
        }
    }

    /*The array length must be a power of two. The array size is [L][2],
    where each sample is complex; array[n][0] is the real part, array[n][1] is the imaginary part of sample n.
    */
    public double[][] fft_1d(double[][] array) {
        double u_r, u_i, w_r, w_i, t_r, t_i;
        int ln, nv2, k, l, le, le1, j, ip, i, n;

        n = array.length;
        ln = (int) (Math.log((double) n) / Math.log(2) + 0.5);
        nv2 = n / 2;
        //输入值换位
        j = 1;
        for (i = 1; i < n; i++) {
            if (i < j) {
                t_r = array[i - 1][0];
                t_i = array[i - 1][1];
                array[i - 1][0] = array[j - 1][0];
                array[i - 1][1] = array[j - 1][1];
                array[j - 1][0] = t_r;
                array[j - 1][1] = t_i;
            }
            k = nv2; // k = nv2 = n/2
            while (k < j) {
                j = j - k;
                k = k / 2;
            }
            j = j + k;
        } //end for i=1~n

        for (l = 1; l <= ln; l++) /* loops thru stages */ {
            le = (int) (Math.exp((double) l * Math.log(2)) + 0.5);
            le1 = le / 2;
            u_r = 1.0;
            u_i = 0.0;
            w_r = Math.cos(Math.PI / (double) le1);
            w_i = -Math.sin(Math.PI / (double) le1);
            for (j = 1; j <= le1; j++) /* loops thru 1/2 twiddle values per stage */ {
                for (i = j; i <= n; i += le) /* loops thru points per 1/2 twiddle */ {
                    ip = i + le1;
                    t_r = array[ip - 1][0] * u_r - u_i * array[ip - 1][1];
                    t_i = array[ip - 1][1] * u_r + u_i * array[ip - 1][0];

                    array[ip - 1][0] = array[i - 1][0] - t_r;
                    array[ip - 1][1] = array[i - 1][1] - t_i;

                    array[i - 1][0] = array[i - 1][0] + t_r;
                    array[i - 1][1] = array[i - 1][1] + t_i;
                }
                t_r = u_r * w_r - w_i * u_i;
                u_i = w_r * u_i + w_i * u_r;
                u_r = t_r;
            }
        } //end of  for loops thru stages

        /*
        for(i=0;i<n;i++){
	        array[i][0]=array[i][0]/(double)N;
	        array[i][1]=array[i][1]/(double)N;
		   } //2003.12.9
         */
        
        return array;
    } /* end of FFT_1d method. */


    //FFT变换，时域数值输入，经过fft_ld计算得出复值，转变为幅值数组
    public double[] fft_sy(double[] input_array) {

        inFFT = new double[N][2];
        outFFT = new double[N][2];
        output=new double[N];
        output_pow = new double[N];

        for (int i = 0; i < N; i++) {
            inFFT[i][0] = input_array[i];
            inFFT[i][1] = 0;
            
        }

        outFFT = fft_1d(inFFT); //use fft method to get the complex after invert

        
        for(int i=0;i<N;i++){
	       
	        output_pow[i] =Math.pow(outFFT[i][0],2)+Math.pow(outFFT[i][1],2); //平方和	        
	        output[i]=Math.sqrt(output_pow[i]);//开根号
	     }

         
       /*
        double a=0,b=0;   
         for(int i=0;i<(N/2);i++){
	        
	         if(output[i]>a){a=output[i];  b=i;}
	         
	         if(output[i]<0.0000000001){
		         System.out.println(i + "  " + "0.0");
		         }
	         else
	         
	         System.out.println(i + "  " + output[i]);
	         }  //2003.12.9
         System.out.println(" the max pl is "+b + "  " + a);
             */
         
         return output;
    }

    

//周期图法功率谱计算，时域数值输入，输出幅值功率谱   
    public double[] powerspetrum(double[] input_array) { 

        inFFT = new double[N][2];
        outFFT = new double[N][2];
        output=new double[N];
        output_pow = new double[N];

        for (int i = 0; i < N; i++) {
            inFFT[i][0] = input_array[i];
            inFFT[i][1] = 0;
        }

        outFFT = fft_1d(inFFT); //use fft method to get the complex after invert
        
        for(int i=0;i<N;i++){
	        output_pow[i] =Math.pow(outFFT[i][0],2)+Math.pow(outFFT[i][1],2); //平方和	        
	        output[i]=output_pow[i]/N;
	     }
      
         return output;
    }




    //倒频谱计算，时域数值输入，输出？？     
    public double[] cepstrum(double[] input_array) {

        inFFT = new double[N][2];
        outFFT = new double[N][2];
        output=new double[N];
        output_pow = new double[N];

        for (int i = 0; i < N; i++) {
            inFFT[i][0] = input_array[i];
            inFFT[i][1] = 0;
        }

        outFFT = fft_1d(inFFT); //use fft method to get the complex after invert
        
        for(int i=0;i<N;i++){
	        output_pow[i]=Math.pow(outFFT[i][0],2)+Math.pow(outFFT[i][1],2); //平方和	        
	        output_pow[i]=output_pow[i]/N;
	        output_pow[i]=Math.log(output_pow[i])/Math.log(10);

	        inFFT[i][0] = output_pow[i];
            inFFT[i][1] = 0;
	     }

         outFFT = fft_1d(inFFT);
  
        for(int i=0;i<N;i++){
	        output[i]=outFFT[i][0];
	        }
         return output;
    }

}