ウェーブレット変換で時間ごとの周波数解析

フーリエ変換では、信号すべての周波数特性を得ます。時系列での周波数特性をみるのは得意ではありません。
そこで、今回は時間ごとの周波数特性をみるためにウェーブレット変換をやってみました。
赤いのが強い部分、青いのが弱い部分で、上にいくほど周波数が高くなります。

ウェーブレット変換では、前後を絞った正弦波など解析の種になる波を使って、周波数ごと時間ごとの特性を解析していきます。
このような、解析につかう波を、マザーウェーブレットといいます。ウェーブレット解析では、周波数ごとに絞りの長さや振幅を変えて、マザーウェーブレットをずらしながら時間ごとに解析していきます。
今回はマザーウェーブレットとして、ガボールのマザーウェーブレットと呼ばれるものを使いました。
解析する信号は、MP3ファイルを読み込んで、4秒後からのデータを使います。

//Wavelet.java
import java.awt.*;
import java.awt.geom.GeneralPath;
import java.awt.image.BufferedImage;
import java.io.InputStream;
import javax.sound.sampled.*;
import javax.swing.*;
import static java.lang.Math.*;

public class Wavelet {
    public static void main(String[] args) throws Exception{
        
        //MP3読み込みの準備
        InputStream is = Mp3.class.getResourceAsStream("/test.mp3");

        AudioInputStream in = null;
        AudioInputStream din = null;
        in = AudioSystem.getAudioInputStream(is);
        AudioFormat format = in.getFormat();

        AudioFormat decodedFormat = new AudioFormat(
                AudioFormat.Encoding.PCM_SIGNED,
                format.getSampleRate(),
                16,
                format.getChannels(),
                format.getChannels() * 2,
                format.getSampleRate(),
                false);

        din = AudioSystem.getAudioInputStream(decodedFormat, in);
        
        double[] data = new double[512];
        byte[] buf = new byte[4];
        //4秒飛ばす
        for(int i = 0; i < format.getSampleRate() * 4
                && din.read(buf, 0, buf.length) != -1; ++i){}
        //データ読み込み
        double max = 0;
        for(int i = 0; i < data.length && din.read(buf, 0, buf.length) != -1; ++i){
            int left = buf[1] * 256 + buf[0];
            //int right = buf[3] * 256 + buf[2];//今回は使わない
            data[i] = left / 32768.;
            if(max < abs(data[i])) max = abs(data[i]);
        }
        //最大値で正規化
        for(int i = 0; i < data.length; ++i){
            data[i] /= max;
        }
        
        double[][] p = new double[25][data.length];
        max = 0;
        for(int fi = 0; fi < p.length; ++fi){
            double a = data.length / pow(2, fi / 3.);//1/3オクターブごと
            //ガボールのウェーブレットを作成
            double sigma = 2;
            double winrate = 0.01;
            int offset = (int)(sqrt(-2 * sigma * sigma * a * a * log(winrate)));
            int gaborlen = offset * 2 + 1;
            double[] gaborre = new double[gaborlen];
            double[] gaborim = new double[gaborlen];
            for(int i = 0; i < gaborlen; ++i){
                double x = (i - offset) / a;
                double d = exp(-x * x /  (2 * sigma * sigma)) 
                        / sqrt(2 * PI * sigma * sigma);
                gaborre[i] = d * cos(PI * x) / sqrt(a);
                gaborim[i] = d * sin(PI * x) / sqrt(a);
            }
            //解析
            for(int i = 0; i < data.length; ++i){
                double re = 0;
                double im = 0;
                for(int gb = 0; gb < gaborlen; ++gb){
                    int gg = gb - gaborlen / 2 + i;
                    double d = (gg < 0 || gg >= data.length) ? 0 : data[gg];
                    re += gaborre[gb] * d;
                    im += gaborim[gb] * d;
                }
                p[fi][i] = (sqrt(re * re + im * im));
                if(p[fi][i] > max) max = p[fi][i];
            }
        }
        
        //描画
        BufferedImage img = new BufferedImage(400, 300, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = (Graphics2D) img.getGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 400, 300);
        g.setColor(Color.BLACK);

        //元データ
        GeneralPath gp = new GeneralPath();
        gp.moveTo(30, 50);
        for(int i = 0; i < data.length; ++i){
            int x = i * 5 / 8 + 30;
            gp.lineTo(x, -data[i] * 40 + 50);
        }
        g.draw(gp);
        
        //周波数分布の描画
        for(int fi = 0; fi < p.length; ++fi){
            int y = 300 - fi * 200 / p.length;
            for(int i = 0; i < p[fi].length; ++i){
                float d = (float) (p[fi][i] / max);
                float r = 1 - d;
                float gr = (d < .5) ? d * 2 : (1 - d) * 2;
                float b = d;
                g.setColor(new Color(b, gr, r));       
                int x = i * 5 / 8 + 30;
                g.fillRect(x, y, 1, 8);
            }
        }
        
        //表示
        JFrame f = new JFrame("ウェーブレット変換");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(400, 350);
        f.add(new JLabel(new ImageIcon(img)));
        f.setVisible(true);       
    }
}