離散ウェーブレット変換の実装

ウェーブレット変換は、時系列データの時間ごとの周波数成分を解析するための手法です。
以前にもウェーブレット変換は やってたのだけど、今回は計算の軽い離散ウェーブレット変換をやってみます。

計算としては、隣り合う2項目の移動差分を値として使い、移動平均をオクターブ下の解析に使うという感じ。
結果、こうなりました。

ところで、解説書としてこれを読んでたのだけど、今は絶版なんですね。

やり直しのための通信数学―フーリエ変換からウェーブレット変換へ (ディジタル信号処理シリーズ)

8要素の数列のウェーブレット変換の手順が書いてあって、すごく具体的にわかりやすくていいのだけど。これ書名がよくないですよね。「通信数学」って、なんか通信教育っぽくて、本屋でみても、まさかウェーブレットの解説本だとはだれも思わない気がします。

コードはこんな感じ。MP3の読み込みにはMP3SPIが必要なのでcom.googlecode.soundlibs:mp3spi:1.9.5.4あたりをdependencyに突っ込んでおく必要があります。

import java.awt.Color;
import java.awt.Graphics2D;
import java.awt.geom.GeneralPath;
import java.awt.image.BufferedImage;
import java.io.File;
import java.util.*;
import java.util.stream.Collectors;
import javax.sound.sampled.*;
import javax.swing.*;

public class DiscreteWavelet {
    public static void main(String[] args) throws Exception {
        AudioInputStream ais = AudioSystem.getAudioInputStream(new File(
                "C:\\Music\\Kiko Loureiro\\No Gravity\\"
                        + "08 - Moment Of Truth.mp3"));//適当なファイルに
        AudioFormat format = ais.getFormat();
        
        AudioFormat decodedFormat = new AudioFormat(
                AudioFormat.Encoding.PCM_SIGNED,
                format.getSampleRate(),
                16,
                format.getChannels(),
                format.getFrameSize(),
                format.getFrameRate(),
                false); // little endian
        AudioInputStream decoded = AudioSystem.getAudioInputStream(decodedFormat, ais);
        
        double[] data = new double[1024];
        byte[] buf = new byte[4];
        
        // 4秒とばす(曲の冒頭はおもしろくない)
        for(int i = 0; i < decodedFormat.getSampleRate() * 4
                && decoded.read(buf, 0, buf.length) != -1; ++i) {}

        // 読み込み
        for (int i = 0; i < data.length && decoded.read(buf, 0, buf.length) != -1; ++i) {
            int left = buf[1] * 256 + buf[0];
            data[i] = left / 32_768.;
        }
        double max = Arrays.stream(data).map(Math::abs).max().orElse(1);
        double[] normalized = Arrays.stream(data).map(d -> d / max).toArray();
        
        // 離散ウェーブレット変換
        List<List<Double>> wavlets = new ArrayList<>();
        List<Double> v = Arrays.stream(normalized).mapToObj(Double::valueOf)
                .collect(Collectors.toList());
        while(v.size() > 1) {
            List<Double> wavlet = new ArrayList<>();
            List<Double> next = new ArrayList<>();
            for (int i = 0; i < v.size(); i += 2) {
                wavlet.add((v.get(i) - v.get(i + 1)) / 2); // 移動差分がwavlet値
                next.add((v.get(i) + v.get(i + 1)) / 2); // 移動平均をつぎにまわす
            }
            wavlets.add(wavlet);
            v = next;
        }
        wavlets.add(v);
        
        // 描画
        BufferedImage image = new BufferedImage(600, 400, BufferedImage.TYPE_INT_RGB);
        Graphics2D g = image.createGraphics();
        g.setColor(Color.WHITE);
        g.fillRect(0, 0, 600, 400);
        g.setColor(Color.BLACK);
        
        // 元データ
        GeneralPath org = new GeneralPath();
        org.moveTo(30, 100);
        for (int i = 0; i < data.length; ++i) {
            int x = i / 2 + 30;
            org.lineTo(x, -data[i] * 100 + 100);
        }
        g.draw(org);
        g.drawLine(30, 100, data.length / 2 + 30, 100);
        
        // ウェーブレット
        double wmax = wavlets.stream().flatMap(List::stream)
                .mapToDouble(Double::doubleValue)
                .map(Math::abs)
                .max().orElse(1);
        
        for (int y = 0; y < wavlets.size(); ++y) {
            List<Double> wavlet = wavlets.get(y);
            int width = data.length / wavlet.size();
            for (int i = 0; i < wavlet.size(); ++i) {
                float d = (float)(wavlet.get(i) / wmax / 2 + .5);
                float b = 1 - d;
                float gr = (d < .5) ? d * 2 : (1 - d) * 2;
                float r = d;
                g.setColor(new Color(r, gr, b));
                g.fillRect(i * width / 2 + 30, y * 20 + 180, width / 2 , 20);
            }
        }
        
        // ウィンドウ表示
        JFrame f = new JFrame("MP3 wave");
        f.setDefaultCloseOperation(JFrame.EXIT_ON_CLOSE);
        f.setSize(600, 450);
        f.add(new JLabel(new ImageIcon(image)));
        f.setVisible(true);
    }
}