javaを用いた短時間フーリエ変換

package application;

import java.io.BufferedReader;
import java.io.File;
import java.io.FileInputStream;
import java.io.IOException;
import java.io.InputStreamReader;
import java.nio.ByteBuffer;
import java.nio.ByteOrder;
import java.util.Arrays;

import javax.sound.sampled.AudioFormat;
import javax.sound.sampled.AudioInputStream;
import javax.sound.sampled.AudioSystem;

import javafx.application.Application;
import javafx.scene.Node;
import javafx.scene.Scene;
import javafx.scene.chart.LineChart;
import javafx.scene.chart.NumberAxis;
import javafx.scene.chart.XYChart;
import javafx.scene.chart.XYChart.Series;
import javafx.scene.layout.HBox;
import javafx.scene.layout.VBox;
import javafx.stage.Stage;
/**
 * 音声(wav)データの波形を見るプログラム
 * ただし、Wav(PCM・リトルエディアン)形式で保存された
 * ファイルのチャンネル１のみ出力
 *
 * @author karura
 *
 */
public class Main extends Application
{
    // 定数
    private final String    fileName    = "/onnsei/0119.wav";        // チャートに表示する音声ファイルへのパス
    private final double    sec         = 1;                     // チャートに表示する期間(s)
    private final double    cut         = 0.05;
    // 取得する音声情報用の変数
    private AudioFormat     format                  = null;
    private double[]       valuesActual            = null;  // 音声データ
    private double[]        spectrumActual          = null;
    private double[]        spectrumImaginal        = null;
    private double[]        spectrumpower           = null;

    public static void main(String[] args) {
        launch(args);
    }

    @Override
    public void start(Stage primaryStage) throws Exception
    {
    	  // 音声ストリームを取得
        File                file    = new File( fileName );
        AudioInputStream    is      = AudioSystem.getAudioInputStream( file );

        // メタ情報の取得
        format = is.getFormat();




        // 取得する標本数を計算
        // 1秒間で取得した標本数がサンプルレートであることから計算
        int mount   = (int) ( format.getSampleRate() *sec );
        int cutmount   = (int) ( format.getSampleRate() *cut );
        System.out.println(mount);
        System.out.println(cutmount);
        int samp   = (int) ( format.getSampleRate() );
        int hsamp = samp/2;
        int t = samp / cutmount;
        System.out.println(t);
        System.out.println(samp);

        // フォント色がおかしくなることへの対処
        System.setProperty( "prism.lcdtext" , "false" );

        // シーングラフの作成
        HBox        root    = new HBox();
        VBox        box1    = new VBox();
        VBox        box2    = new VBox();
        VBox        box3    = new VBox();
        root.getChildren().addAll( box1 , box2 ,box3);
        System.out.println( "loading wav data..." );
       initialize();


       int count = 0;
       double goukei = 0;
       double[] P       = new double[ valuesActual.length ];//短時間フーリエ変換後の一時値
       
       spectrumActual       = new double[ cutmount ];
       spectrumImaginal     = new double[ cutmount ];
       spectrumpower     = new double[ cutmount ];
      int n = (mount+cutmount-1)/cutmount;
      int m = 0;
      int o = cutmount;
      System.out.println("ループ回数"+n);


      for(int k =0; k<n; k++){
     if(o > valuesActual.length){

    	 double[] punc =Arrays.copyOfRange(valuesActual, m,valuesActual.length );
    	 DFT(punc, spectrumActual , spectrumImaginal , spectrumpower,false );

    	    for(int i = 0; i < cutmount; i++){
    	   P[i] = spectrumpower[i] *(t/hsamp);
     }}


     else{
      double[] punc =Arrays.copyOfRange(valuesActual, m,o );
      DFT(punc, spectrumActual , spectrumImaginal , spectrumpower,false );
      for(int s = 0; s < cutmount; s++){
    	P[k] = P[k] + spectrumpower[s] * ((double)t / hsamp) ;

      }






     }
       m = m + cutmount;
       o = o + cutmount;
      count ++;

        // 離散変換後の波形をチャート表



      }
      int plength = 0;
      for(int x =0; x < P.length;x++ ){
    		if(P[x] != 0){
    		plength ++;	
    		}


  }
      
      double[] POWER  = new double[ plength ];//
      for(int x =0; x < plength;x++ ){
  		if(P[x] != 0){
  			POWER[x]= P[x];
  		}}
    //powerの最大値，最小値，平均値の算出
      
      
      
      
      
      
      
      
     //powerの最大値，最小値，平均値の算出
	     double[] heikin = new double[20000];
	     double[] max = new double[20000];
	     double[] min = new double[20000];
	     double minh = 10000;
	     double maxh = 0;
	     int counts =0;
	     int countt=0;
	     double goukeis=0;


	     for(int i = 0; i< POWER.length;){
	    		
	    	 if(POWER[i] <= 1.75){//閾値1.75として単語ごとに分割
	    		 if(i != 0){
	    		 heikin[countt] = goukeis /counts;
	    		 max[countt] = maxh;
	    		 min[countt] = minh;
	    		 countt++;
	    		 }
	    	 while(POWER[i] <= 1.75){
	    		 i++;
	    		 if(i == POWER.length-1){
	    			 break;
	    		 }
	    		 }
	    	 goukeis =0;
	    	 counts =0;
	    	 maxh = 0;
	    	 }
	    	 if(maxh < POWER[i]){
	    		 maxh = POWER[i];
	    	 }
	    	 if(minh > POWER[i]){
	    		 minh = POWER[i];
	    	 }
	    	 goukeis += POWER[i];
	    	 counts++;
	    	 
	    	 
	    	 if(i == POWER.length-1){
	    		 max[countt] = maxh;
	    		 min[countt] = minh;
	    		 heikin[countt] = goukeis /counts; 
	    		 System.out.println(heikin[countt]);
	    		 break;
	    	 }
	    	 i++;
	    	
	    	
	    	
	     }
	     
	     for(int u = 0;u<=countt; u++){
	    	 System.out.println("平均値"+heikin[u]);
	    	 System.out.println("最大値"+max[u]);
	    	 System.out.println("最小値"+min[u]);
	     }
	     
	     
	     
	     
	     
	

 for(int y =0; y < POWER.length;y++ ){
	
		 System.out.println(POWER[y]);


 }
 
 
 

        box1.getChildren().add( createLineChart( "音声波形" , valuesActual ) );
        // 離散フーリエ変換後のスペクトルをチャート表示
        box2.getChildren().add( createLineChart( "スペクトル(実数部)" , spectrumActual ) );            // 折れ線グラフの追加
        box2.getChildren().add( createLineChart( "スペクトル(虚数部)" , spectrumImaginal ) );          // 折れ線グラフの追加
        box1.getChildren().add( createLineChartnew( "パワースペクトル" ,  P,n) );

        // シーンの作成
        Scene       scene   = new Scene( root , 800 , 400 );


        // ウィンドウ表示
        primaryStage.setScene( scene );
        primaryStage.show();

    }

    /**
     * 音声ファイルを読み込み、メタ情報とサンプリング・データを取得
     * @throws Exception
     */
    protected void initialize() throws Exception
    {
        // 音声ストリームを取得
        File                file    = new File( fileName );
        AudioInputStream    is      = AudioSystem.getAudioInputStream( file );

        // メタ情報の取得
        format = is.getFormat();



        // 取得する標本数を計算
        // 1秒間で取得した標本数がサンプルレートであることから計算
        int mount   = (int) ( format.getSampleRate() *sec );

        int samp   = (int) ( format.getSampleRate() );
        int t = samp / mount;
        System.out.println(t);

        // 音声データの取得
        valuesActual    = new double[ mount ];

        for( int i=0 ; i<mount ; i++ )
        {
            // 1標本分の値を取得
            int     size        = format.getFrameSize();
            byte[]  data        = new byte[ size ];
            int     readedSize  = is.read(data);

            // データ終了でループを抜ける
            if( readedSize == -1 ){ break; }

            // 1標本分の値を取得
            switch( format.getSampleSizeInBits() )
            {
                case 8:
                    valuesActual[i]   = (int) data[0];
                    break;
                case 16:
                    valuesActual[i]   = (int) ByteBuffer.wrap( data ).order( ByteOrder.LITTLE_ENDIAN ).getShort();
                    break;
                default:
            }
        }

        // 音声ストリームを閉じる
        is.close();
    }

    /**
     * 離散フーリエ変換
     * @param in フーリエ変換を行う実数配列
     * @param outActual 計算結果の実数部配列
     * @param outImaginal 計算結果の虚数部配列
     * @param winFlg 窓関数の使用フラグ
     */
    protected void DFT( double[] in , double[] outActual , double[] outImaginal ,double[] power, boolean winFlg )
    {
        // 配列初期化
        int  length             = in.length;

        // 離散フーリエ変換
        for( int k=0 ; k<length ; k++ )
        {
            // 初期化

            outActual[k]    = 0.0d;
            outImaginal[k]  = 0.0d;

            // 計算
            for( int n=0 ; n<length ; n++ )
            {
                // 入力値に窓関数を適用
                double normal   = ( !winFlg )? in[n]  : hanWindow( in[n] , n , 0 , length );

                // k次高周波成分を計算
                outActual[k]    +=        normal * Math.cos( 2.0 * Math.PI * (double)n * (double)k / (double)length );
                outImaginal[k]  += -1.0 * normal * Math.sin( 2.0 * Math.PI * (double)n * (double)k / (double)length );
            }



            // 残りの計算
            outActual[k]    /= length;
            outImaginal[k]  /= length;
            power[k] = Math.sqrt(Math.pow(outActual[k],2)+Math.pow(outImaginal[k],2));
        }
    }



    /**
     * 窓関数（ハン窓）
     * @param in 変換する値
     * @param i 配列中のインデックス
     * @param minIndex 配列の最小インデックス
     * @param maxIndex 配列の最大インデックス
     * @return
     */
    protected double hanWindow( double in , double i , double minIndex , double maxIndex )
    {
        // 入力値の正規化
        double normal   = i / ( maxIndex - minIndex );


        // ハン窓関数の値を取得
        double  han     =  0.5 - 0.5 * Math.cos( 2.0 * Math.PI * normal );

        return in * han;
    }



    /**
     * 折れ線グラフで波形表示
     * @param title グラフのタイトル文字
     * @param values グラフに出力するデータ配列
     * @return 作成したグラフノード
     */
    @SuppressWarnings("unchecked")
    protected Node createLineChart( String title , double[] values )
    {
        // 折れ線グラフ
        NumberAxis                  xAxis   = new NumberAxis();
        NumberAxis                  yAxis   = new NumberAxis();
        LineChart<Number, Number>   chart   = new LineChart<Number, Number>( xAxis , yAxis );
        chart.setMinWidth( 600 );
        chart.setMinHeight( 400 );

        // データを作成
        Series< Number , Number > series1    = new Series<Number, Number>();
        series1.setName( title  );
        for( int i=1 ; i<values.length ; i++ )
        {

        	double m = ((i /sec));
            series1.getData().add( new XYChart.Data<Number, Number>( m , values[i] ) );
        }

        // データを登録
        chart.getData().addAll( series1 );

        // 見た目を調整
        chart.setCreateSymbols(false);                                                          // シンボルを消去
        series1.getNode().lookup(".chart-series-line").setStyle("-fx-stroke-width: 0.75px;");   // 線を細く

        return chart;
    }
    protected Node createLineChartnew( String title , double[] values ,int t)
    {

        // 折れ線グラフ
        NumberAxis                  xAxis   = new NumberAxis();
        NumberAxis                  yAxis   = new NumberAxis();
        LineChart<Number, Number>   chart   = new LineChart<Number, Number>( xAxis , yAxis );
        chart.setMinWidth( 600 );
        chart.setMinHeight( 400 );

        // データを作成
        Series< Number , Number > series1    = new Series<Number, Number>();
        series1.setName( title  );
        double time = 0;

        for(  int i =0; i< t ; i++ )
        {


            series1.getData().add( new XYChart.Data<Number, Number>(time , values[i] ) );
           time = time + cut;

        }

        // データを登録
        chart.getData().addAll( series1 );

        // 見た目を調整
        chart.setCreateSymbols(false);                                                          // シンボルを消去
        series1.getNode().lookup(".chart-series-line").setStyle("-fx-stroke-width: 0.75px;");   // 線を細く

        return chart;
    }}

このサイトを元にJavaで離散フーリエ変換を少し書き加え短時間フーリエ変換を行うようにした。
パワースペクトルを算出し，時間ごとの音声のパワーの変化をグラフ表示できるようになっている。
少し実行時間がかかってしまうのが難点である。
OSはiOSを用いている。

参照サイト: Javaで周波数分析をしてみる｜軽Lab