C#で離散ウェーブレット変換 - 真実の楽譜（フルスコア）

Scilabで離散ウェーブレット変換
 http://truthfullscore.hatenablog.com/entry/2014/01/17/172746

をC#で実装してみた。

ドベシィウェーブレット変換プログラム

DaubechiesWavelet.cs

using System;

namespace Wavelet
{
  class DaubechiesWavelet
  {
    //Daubechies係数を定義
    public static double[][] DaubechiesCoefficients =
    { 
      new double[] {1.0/Math.Sqrt(2.0),1.0/Math.Sqrt(2.0)},
      new double[] {0.6830127/Math.Sqrt(2.0),1.1830127/Math.Sqrt(2.0),0.3169873/Math.Sqrt(2.0),-0.1830127/Math.Sqrt(2.0)},
      new double[] {0.47046721/Math.Sqrt(2.0),1.14111692/Math.Sqrt(2.0),0.650365/Math.Sqrt(2.0),-0.19093442/Math.Sqrt(2.0),-0.12083221/Math.Sqrt(2.0),0.0498175/Math.Sqrt(2.0)},
      new double[] {0.32580343/Math.Sqrt(2.0),1.01094572/Math.Sqrt(2.0),0.8922014/Math.Sqrt(2.0),-0.03957503/Math.Sqrt(2.0),
              -0.26450717/Math.Sqrt(2.0),0.0436163/Math.Sqrt(2.0),0.0465036/Math.Sqrt(2.0),-0.01498699/Math.Sqrt(2.0)},
      new double[] {0.22641898/Math.Sqrt(2.0),0.85394354/Math.Sqrt(2.0),1.02432694/Math.Sqrt(2.0),0.19576696/Math.Sqrt(2.0),
              -0.34265671/Math.Sqrt(2.0),-0.04560113/Math.Sqrt(2.0),0.10970265/Math.Sqrt(2.0),-0.0088268/Math.Sqrt(2.0),
              -0.01779187/Math.Sqrt(2.0),4.72E-03/Math.Sqrt(2.0)}
    };

    static void Main(string[] args)
    {
      int totalLength = 1024;  //データ数（2のべき乗）
      int maxLevel = 2;   //分解レベル
      int N = 1;          //次数
      int kMax = 2 * N;
      int daubechiesIndex = N - 1;

      double[] inputSignal = new double[totalLength];
      double[] outputSignal = new double[totalLength];

      //確認用に入力データをCSVファイル出力
      System.IO.StreamWriter swriter = new System.IO.StreamWriter("input.csv");
      //初期化
      for (int signalIndex = 0; signalIndex < totalLength; signalIndex++)
      {
        inputSignal[signalIndex] = 10.0 * Math.Sin(signalIndex * 2.0 * Math.PI / (totalLength / 2.0))
                      + 1.0 * Math.Sin(signalIndex * 2.0 * Math.PI / (totalLength / 200.0));
        outputSignal[signalIndex] = 0.0;
        swriter.WriteLine(inputSignal[signalIndex]);
      }
      swriter.Close();

      for (int level = 1; level <= maxLevel; level++)
      {
        int scale = 1;
        for (int i = 0; i < level; i++)
        {
          scale *= 2;
        }

        //ダウンサンプリング
        for (int index = 0; index < totalLength / scale; index++)
        {
          int setApproximationIndex = index;
          int setDetailIndex = index + totalLength / scale;

          outputSignal[setApproximationIndex] = 0.0;
          outputSignal[setDetailIndex] = 0.0;

          for (int k = 0; k < kMax; k++)
          {
            int getInputIndex = 2 * index + k;

            if (getInputIndex >= 2 * totalLength / scale)
            {
              //配列外にデータにアクセスする際に折り返しを行う
              int over = getInputIndex % (2 * totalLength / scale) + 1;
              getInputIndex = (2 * totalLength / scale) - over;
            }

            //フィルタリング
            outputSignal[setApproximationIndex] += inputSignal[getInputIndex] * DaubechiesCoefficients[daubechiesIndex][k];
            outputSignal[setDetailIndex] += inputSignal[getInputIndex] * Math.Pow(-1.0, k) * DaubechiesCoefficients[daubechiesIndex][kMax - 1 - k];
          }
        }

        //配列コピー
        Array.Copy(outputSignal, inputSignal, totalLength);

        //結果出力
        swriter = new System.IO.StreamWriter("signalLevel"+level.ToString()+".csv");
        for (int signalIndex = 0; signalIndex < totalLength / scale; signalIndex++)
        {
          swriter.WriteLine(outputSignal[signalIndex] + "," + outputSignal[signalIndex + totalLength / scale]);
        }
        swriter.Close();
      }
    }
  }
}

解説

入力信号作成

前回の記事同様に、
振幅10, 2Hzの正弦波と、
振幅1 , 100Hzの正弦波を足し合わせた信号を作成。

      System.IO.StreamWriter swriter = new System.IO.StreamWriter("input.csv");
      //初期化
      for (int signalIndex = 0; signalIndex < totalLength; signalIndex++)
      {
        inputSignal[signalIndex] = 10.0 * Math.Sin(signalIndex * 2.0 * Math.PI / (totalLength / 2.0))
                      + 1.0 * Math.Sin(signalIndex * 2.0 * Math.PI / (totalLength / 200.0));
        outputSignal[signalIndex] = 0.0;
        swriter.WriteLine(inputSignal[signalIndex]);
      }
      swriter.Close();
||

フィルタリングとダウンサンプリング

離散ウェーブレット変換は、元の信号x(n)を下の図のように、
ハイパスとローパスの線形フィルタリングを実行。
その後フィルタリング処理データを1点おきに取得して、
信号をダウンサンプリング。
f:id:s_sikisya:20140122131835p:plain
このときに表れた低周波信号をApproximation
もう一方の高周波信号をDetailと呼ぶ。

そして指定した分解レベルになるまで、
低周波信号に対して再帰的に
フィルタリングとダウンサンプリングを実行。
f:id:s_sikisya:20140122130900p:plain

最終的に得られた信号列

Approximation Level N, Detail Level Detail N, Level N-1, … , Detail Level 2, Detail Level 1

がScilabの以下の式で求めたウェーブレット係数Cに該当します。

[C,L]=wavedec(wave,2,"haar");

フィルタリング式

これで離散ウェーブレット変換は入力信号から、
ApproximationとDetailを求めるフィルタリング式さえあれば、
順次求めることができそう。

フィルタリング式は、
DaubechiesWaveletの場合それぞれ、

${ \displaystyle g_n = \sum_{k=0}^{2N-1} a_k x(n-k) }$

${ \displaystyle h_n = \sum_{m=0}^{2N-1} b_k x(n-k) }$

となります。

Daubechies係数

また係数には以下の関係式が成り立つ。

${ \displaystyle b_k = {(-1)}^{m} a_{2N-1-k} }$

${ \displaystyle \sum_{m=0}^{2N-1} a_k = \sqrt{2} }$

${ \displaystyle \sum_{m=0}^{2N-1} a_k^2 = 1 }$

このDaubechies係数を求める方法があるようですが、
一般的には既に求められているのを使う模様。
英語版Wikipediaに次数N=10の場合までの係数が用意されています。

a_k	D2 (Haar)	D4	D6	D8	D10
0	1	0.6830127	0.47046721	0.32580343	0.22641898
1	1	1.1830127	1.14111692	1.01094572	0.85394354
2		0.3169873	0.650365	0.8922014	1.02432694
3		-0.1830127	-0.19093442	-0.03957503	0.19576696
4			-0.12083221	-0.26450717	-0.34265671
5			0.0498175	0.0436163	-0.04560113
6				0.0465036	0.10970265
7				-0.01498699	-0.0088268
8					-0.01779187
9					0.004717428

Daubechies wavelet - Wikipedia, the free encyclopediaより一部(N=5まで)引用

注意
この表は係数の二乗和が２で正規化されているので、
全てルート２で割った値を使用すること。

結果出力

上記プログラムは分解レベル毎に
ApproximationとDetailをCSVファイルで保存しています。
ウェーブレット変換結果をExcelでプロットとした図。

レベル１
f:id:s_sikisya:20140205171527p:plain
レベル２
f:id:s_sikisya:20140205171546p:plain

以前行ったScilabで離散ウェーブレット変換
での結果同様にLevel2の段階でDetailに高周波成分が出てきています。
f:id:s_sikisya:20140117131947p:plain
サンプリングができない点の処理（プログラムでは折り返し）が、
異なるのかデータの端だけScilabと値が異なるけど今回は気にしない方向で。

参考
離散ウェーブレット変換 - Wikipedia
Daubechies wavelet - Wikipedia, the free encyclopedia
直交ウェーブレット変換について