真実の楽譜(フルスコア)

プログラム関係の忘備録になるはず

Open JTalkをWindowsでビルドして動かしてみる

オープンソースの日本語の音声合成エンジンOpen JTalk
Open JTalk - HMM-based Text-to-Speech System
日本語合成エンジンなのに公式解説は英語だらけだし、
Windowsに導入する際にはソースファイルからビルドが必要だったりと、
一癖あったので忘備録として導入手順をご紹介。

準備

必要なもの

  • Open JTalk本体, open_jtalk-1.xx (記述時1.08)
  • 解析用辞書ファイル, open_jtalk_dic-1.xx (記述時1.08)
  • 音声合成エンジン, hts_engine_API-1.xx (記述時1.09)
  • 音響モデルデータ, hts_voice_nitech_jp_atr503_m001-1.xx (記述時1.05)
  • Visual Studio Community 2013 (VC++12.0)

ファイルのダウンロード

http://sourceforge.net/projects/hts-engine/から以下のファイルをダウンロード

  • hts_engine_API-1.09.tar.gz

http://sourceforge.net/projects/open-jtalk/から以下のファイルをダウンロード。

  • open_jtalk-1.08.tar.gz
  • open_jtalk_dic_shift_jis-1.08.tar.gz
  • hts_voice_nitech_jp_atr503_m001-1.05.tar.gz

辞書ファイル(open_jtalk_dic)は文字コード毎に用意されている模様で、
今回はWindows環境なのでshift_jisを選択。
他の文字コードだとどうなるか不明。

ビルド

hts_engine_APIのビルド

hts_engine_APIがないとOpen JTalkをビルドできないので、
先にhts_engine_APIを以下の手順でビルドします。

  1. ダウンロートした「hts_engine_API-1.09.tar.gz」を解凍
  2. コマンドプロンプトを起動し解凍先のフォルダへ移動後以下のコマンドを入力

call "C:\Program Files (x86)\Microsoft Visual Studio 12.0\VC\vcvarsall.bat"
nmake -f Makefile.mak
nmake -f Makefile.mak install

ビルドに成功すると以下のパスにhts_engine_APIが作成されます。

C:\hts_engine_API

環境変数のPATHにダブルクォーテーション["]が含まれていると、
 「vcvarsall.bat」を呼び出した時点でエラーが発生してビルドが行われません。
 環境変数設定のGUI上では入力できないのですが、自分の環境ではDirectX関連のパスが
 何故かダブルクォーテーションで囲まれていたのでご注意を。

Open JTalkのビルド

続いてOpen JTalkを同様にビルドします。

  1. ダウンロートした「open_jtalk-1.08.tar.gz」を解凍
  2. コマンドプロンプトを起動し解凍先のフォルダへ移動後以下のコマンドを入力

call "C:\Program Files (x86)\Microsoft Visual Studio 12.0\VC\vcvarsall.bat"
nmake -f Makefile.mak
nmake -f Makefile.mak install

ビルドに成功すると以下のパスにOpen JTalkの実行ファイルが作成されます。

C:\open_jtalk

音声合成

辞書・音響モデルの配置

Open JTalkはコマンドライン実行プログラムであり、
コマンドライン引数として辞書データファイルと音響モデルファイルを指定する必要があるため、
使いやすくするために予め実行ファイルディレクトリに配置しておきます。

  1. ダウンロートした「open_jtalk_dic_shift_jis-1.08.tar.gz」を解凍
  2. 解凍したフォルダ名を「dic」にして「C:\open_jtalk\bin」直下に配置
  3. ダウンロートしたhts_voice_nitech_jp_atr503_m001-1.05.tar.gzを解凍
  4. 解凍したフォルダ内の「nitech_jp_atr503_m001.htsvoice」を「C:\open_jtalk\bin」直下に配置
入力テキストファイルの作成

テキストエディタを開き音声合成したい文字を入力し、
shift_jis文字コードで「C:\open_jtalk\bin\Input.txt」に保存します。

input.txt例

こんにちは、お元気ですか

データ配置例

f:id:s_sikisya:20150529162309j:plain
f:id:s_sikisya:20150529162319j:plain

実行

コマンドプロンプトを起動し以下のコマンドを入力

open_jtalk.exe -m nitech_jp_atr503_m001.htsvoice -x dic -ow output.wav input.txt

同一フォルダ内に「output.wav」の音声合成した音声ファイルの出来上がり〜
音声のモデルはどこぞやの男性っぽいですね。


声質を変化させる場合は、コマンドライン引数のパラメータを弄ればいいけど、
別の人の声に変えるにはhtsvoiceファイルを用意しないと無理そうです。
名工大の音声案内端末として使われいるメイちゃん声は落とせるようだけど、
好きな声を使うにはモデル生成する必要があるようでまだまだ謎が多いです。

C#でボックスミューラー法による正規分布に従う乱数生成

.NET Frameworkでは一様乱数を生成するクラスとして
System.Randomクラスが提供されています。
Random クラス (System)

データにノイズを加えるときに一様乱数よりも、
正規分布に従うノイズ(ホワイトガウス)を加えたいというケースがありますが、
.NET Frameworkではそこまでは提供してくれていないので、
ボックスミューラー法による正規分布に従う乱数生成クラスを作ってみました。

ボックスミューラー法

ボックスミューラー法のアルゴリズムは単純で、
(0~1]の範囲を取る一様乱数を2つ用意し(R1,R2)、
以下の式に当てはめるだけで、
N(0.0,1.0)の正規分布に従う乱数が1組生成できます。

{ \displaystyle
Norm_1 = \sqrt{-2 \log R_1} \cos{2 \pi R_2}
}
{ \displaystyle
Norm_2 = \sqrt{-2 \log R_1} \sin{2 \pi R_2}
}

また平均値μ、分散σ^2に従う正規分布N(μ,σ^2)での乱数は、
得られた値にσを乗算して、平均値を加えることで求まります。

{ \displaystyle
Norm_1 = \sigma \sqrt{-2 \log R_1} \cos{2 \pi R_2} + \mu
}
{ \displaystyle
Norm_2 =  \sigma \sqrt{-2 \log R_1} \sin{2 \pi R_2} + \mu
}

参考
http://ja.wikipedia.org/wiki/乱数列
http://ja.wikipedia.org/wiki/ボックス=ミュラー法

プログラムソース

RandomBoxMuller.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace RandomBoxMuller
{
  public class RandomBoxMuller
  {
    private Random random;

    public RandomBoxMuller()
    {
      random = new Random(Environment.TickCount);
    }

    public RandomBoxMuller(int seed)
    {
      random = new Random(seed);
    }

    public double next(double mu = 0.0, double sigma = 1.0, bool getCos = true)
    {
      if (getCos)
      {
        double rand = 0.0;
        while ((rand = random.NextDouble()) == 0.0) ;
        double rand2 = random.NextDouble();
        double normrand = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Cos(2.0 * Math.PI * rand2);
        normrand = normrand * sigma + mu;
        return normrand;
      }
      else
      {
        double rand;
        while ((rand = random.NextDouble()) == 0.0) ;
        double rand2 = random.NextDouble();
        double normrand = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Sin(2.0 * Math.PI * rand2);
        normrand = normrand * sigma + mu;
        return normrand;
      }
    }

    public double[] nextPair(double mu = 0.0, double sigma = 1.0)
    {
      double[] normrand = new double[2];
      double rand = 0.0;
      while ((rand = random.NextDouble()) == 0.0) ;
      double rand2 = random.NextDouble();
      normrand[0] = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Cos(2.0 * Math.PI * rand2);
      normrand[0] = normrand[0] * sigma + mu;
      normrand[1] = Math.Sqrt(-2.0 * Math.Log(rand)) * Math.Sin(2.0 * Math.PI * rand2);
      normrand[1] = normrand[1] * sigma + mu;
      return normrand;
    }
  }
}

実行例

作成したクラスはRandom関数と同様に、
コンストラクタにシード値を指定させてオブジェクトを初期化し、
(未指定の場合はシステムの起動時刻をシード値として利用)
next関数の引数に正規分布のμとσを指定することで乱数が求まります。

Program.cs
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace RandomBoxMuller
{
  class Program
  {
    static void Main(string[] args)
    {
      var random = new RandomBoxMuller();

      using (var sr = new System.IO.StreamWriter("output.csv"))
      {
        for (int i = 0; i < 10000; i++)
        {
          sr.WriteLine(random.next());
        }
      }
    }
  }
}


上記プログラムではN(0.0,1.0)の乱数を1万個生成しています。
得られたデータのヒストグラムは以下のようになりました。

f:id:s_sikisya:20140603132934p:plain

今回-1.0~1.0のデータは6816個得られており、
Z=1.0のときの正規分布データ点数は6826個(0.3413)
なのでほぼ近い値が得られていることが確認できました。

参考
http://www.koka.ac.jp/morigiwa/sjs/standard_normal_distribution.htm

MPU-9150をArduinoに接続して9軸モーションデータを取得する

MPU-9150

MPU-9150はGoogle Grassにも使用されている(らしい)
9軸(加速度・角速度・地磁気)モーションセンサー。
ストロベリーリナックス等でモジュールが購入できます。
https://strawberry-linux.com/catalog/items?code=12150

Arduinoに繋いでセンサーと仲良くなるためのメモ。

Arduinoを使用したデータ取得

単純に9軸データを取得するだけなら、
SparkfunがGitHubで公開している以下のソースを使うと楽。
https://github.com/sparkfun/MPU-9150_Breakout

MPU-9150_Breakoutを使った9軸データ取得

1.GitHubからリポジトリのマスターダウンロード

ページ内のDownload.Zipボタンをクリックして、
MPU-9150_Breakout-master.zip」をダウンロード。

2.ダウンロードしたファイルを解凍しライブラリに追加

MPU-9150_Breakout-master.zip」を解凍し解凍したフォルダから、

firmware\MPU6050
firmware\I2Cdev

の2つをArduino IDEのライブラリフォルダ「Arduino\libraries」にコピー。

3.Arduino IDEでスケッチを送信

MPU-9150_Breakout-master.zip」を解凍し解凍したフォルダから、

firmware\MPU6050\Examples\MPU9150_raw\MPU9150_raw.ino

Arduinoスケッチを書き込めば、
シリアル通信で9軸のデータ値が取得できます。

4.物理量として利用するには

このサンプルでは全てのセンサー値が符号付き16ビットで格納されていますが、
値磁気センサーのみ分解能が13ビットなので注意が必要です。
デフォルトでは以下のようにセンサー情報が表されています。

  • 加速度:±2G
  • 角速度:±250(degree/sec)
  • 地磁気:±1200(μT)

また物理量に変換するなら次の式を利用します。

加速度:{ \displaystyle
Acc = 2 * SensorValue / 32768
}
角速度:{ \displaystyle
{\omega} = 250 * SensorValue  / 32768
}
地磁気{ \displaystyle
Mag = 1200 * SensorValue / 4096
}

ちなみに地磁気は周囲の影響を受けまくっているので、
ぐりぐりと8の字に動かしてこの値からオフセットを求めて
補正することが必要になります。

5.もうちょっとできる子なんですが

実はこのセンサー実はもうちょっとできる子です。
センサー公式に概要が書かれていますがセンサーの中に、
Digital Motion Processing™ (DMP™) エンジンが搭載されており、
センサーフュージョンによる姿勢計算をして出力してくれます。
MEMS Gyro-Accel | Gyroscope | Accelerometer | Processing - MPU-9150 9軸(ジャイロ+加速度+コンパス)MEMS MotionTracking™ デバイス

MPU-9150_Breakout-master内の「MPU6050_DMP6.ino」は
DMPエンジンを利用したサンプルスケッチなんですが、
値磁気センサーが組み込まれていない下位モデルのMPU-6050用で、
このままだと6軸のみで地磁気のデータが取れてこない。

そこでまた色々と試行錯誤したわけですがそれは書く機会があれば。

Android Developer Toolkit が 22.6.1 にアップデートできないときの対処

Android SDKマネージャでAndroid SDK Toolsのアップデートを行ったら、
Android Developer Toolkitも更新するように表示されたので、
更新の確認からEclipseプラグインの更新を行おうとしたところ

No repository found containing: org.eclipse.update.feature,com.android.ide.eclipse.ddms,22.6.1.v201403111859-1066720

とエラーが発生し、更新できない&Androidプロジェクト全てがビルドエラー状態に。
サーバー側にデータが置かれてないのか、パスのスペルミスなのかは不明ですが、
新規にAndroid Developer Toolkitをインストールさせることで解決できた。

手順

  1. [ヘルプ]→[新規ソフトウェアのインストール]をクリック
  2. [追加]ボタンをクリック
  3. ロケーションに「https://dl-ssl.google.com/android/eclipse/」を入力し[OK]をクリック
  4. リスト内に表示された項目にチェックを入れ、[次へ]をクリック
  5. インストールダイアログの指示に従いインストール

ReadOnlyの時にチカチカしないNumericUpDownを作る

NumericUpDown

GUIで数値を上下に変化させて指定することができる、
NumericUpDownというFormコントロールがあります。
NumericUpDown クラス (System.Windows.Forms)

アプリケーションでユーザーに数値を入力してもらう際に、
ボタンを押せば一定間隔で指定ができるので、
開発側としても変な入力がされにくい便利なコントロール。

ReadOnly時のチカチカがいらないよね

f:id:s_sikisya:20140310185141p:plain
ユーザが直接値を指定して変更させないために、
ReadOnlyのプロパティが用意されていますが、
ReadOnly時でもフォーカスに入るとカーソルがチカチカと点滅。
これじゃユーザが入力できると勘違いする恐れがあります。

チカチカを消してしまうには

調べてみるとこのチカチカするヤツ、
caret(キャレット)って名前らしく、
これを消し去るにはWindows APIが必要とのこと。

参考
[C#] 独自に作成したウィンドウコントロールのコンポーネントにキャレットを表示する

チカチカしないNumericUpDownを作る

今回はSystem.Windows.Forms.NumericUpDownクラスを継承して、
ReadOnly時にキャレットを表示させないコントロールを実装してみました。

NumericUpDown継承クラスの作成

NumericUpDownクラスを継承した以下のクラスを作成します。

NumericUpDownNoBlinking.cs
using System;
using System.Runtime.InteropServices;
using System.Windows.Forms;
using System.ComponentModel;

namespace NumericUpDownNoBlinking
{
  public partial class NumericUpDownNoBlinking : NumericUpDown
  {

    [DllImport("user32.dll")]
    static extern bool HideCaret(IntPtr hWnd);

    //キャレット表示の有無を指定
    private bool caretVisible = true;

    //ReadOnlyメソッドを上書き
    public new bool ReadOnly
    {
      get
      {
        return base.ReadOnly;
      }
      set
      {
        base.ReadOnly = value;

        //キャレットの有無をReadOnly属性に応じて変更
        caretVisible = !value;

        //ReadOnly有効時、NumericUpDownのTextBoxコントロールにカーソルを変化させない。
        foreach (Control cntr in this.Controls)
        {
          if (cntr is TextBox)
          {
            if (caretVisible)
            {
              cntr.Cursor = Cursors.IBeam;
            }
            else
            {
              cntr.Cursor = Cursors.Arrow;
            }
          }
        }
      }
    }

    public NumericUpDownNoBlinking()
    {
      InitializeComponent();
    }

    public NumericUpDownNoBlinking(IContainer container)
    {
      container.Add(this);

      InitializeComponent();
    }

    //フォーカス取得時にキャレットを消去
    protected override void OnGotFocus(EventArgs e)
    {
      base.OnGotFocus(e);

      if (!caretVisible)
      {
        //NumericUpDownのTextBoxコントロールのキャレットを非表示にする
        foreach (Control cntr in this.Controls)
        {
          if (cntr is TextBox)
          {
            HideCaret(cntr.Handle);
          }
        }
      }
    }
  }
}

解説

Windwos APIのHideCaret関数を使えるように、
宣言しておきます。

[DllImport("user32.dll")]
static extern bool HideCaret(IntPtr hWnd);

コントロールのOnGotFocusメソッドをオーバーライドし、
フォーカス取得時にControlsプロパティに格納されている
TextBoxのHandleをHideCaret関数の引数に渡せばOK。

protected override void OnGotFocus(EventArgs e)
{
  base.OnGotFocus(e);

  if (!caretVisible)
  {
    //NumericUpDownのTextBoxコントロールのキャレットを非表示にする
    foreach (Control cntr in this.Controls)
    {
      if (cntr is TextBox)
      {
        HideCaret(cntr.Handle);
      }
    }
  }
}

NumericUpDownのHandleではなくTextBoxを指定することがポイント!

また上記プログラムでは、
テキストボックス上にマウスカーソルが移動しても、
マウスポインタの形状を矢印のままにしています。
こちらもControlsプロパティに格納されているTextBoxに対して、
プロパティを変更することで対処しています。

参考
ネコキス: NumericUpDown のキャレットを消す
[C#] 独自に作成したウィンドウコントロールのコンポーネントにキャレットを表示する

C# Fromで動画再生するには

WPFやWinRTならXAMLベースのMediaElementクラスを使うことが、
MSDNでも書かれているように鉄板だと思います。
http://msdn.microsoft.com/ja-jp/library/aa970915(v=vs.110).aspx

しかし現在自分が今行っている案件はWindows Form開発。
Formアプリケーションではどうするのって事になったので、
調べてみた。

WindowsMediaPlayerコントロールを使う

一番簡単に実装するならコレが手っ取り早い模様。
再生、停止、スクロール等のUI操作は、
全部このコントロールにおまかせできます。

コントロール追加方法
  1. VisualStudioから[ツール]→[ツールボックス アイテムの選択]をクリック。
  2. [COMコンポーネント]のタブを選択
  3. [Windows Media Player]にチェックを入れOKをクリック
  4. ツールボックスに追加された[Windows Media Player]をFormに貼り付け

f:id:s_sikisya:20140303095357p:plain
f:id:s_sikisya:20140303095857p:plain

これだけでプレイヤーが画面に追加できます。
再生時にはコントロールのURLプロパティに、
ファイルパスを指定するだけでOK。

axWindowsMediaPlayer1.URL = "moviefile.wmv";

ちょっと気になるところ
  • コマ送り再生ができない
  • WindowsMediaPlayerの外観使ってます感バリバリ
  • UI部分のカスタマイズができない
  • 一時停止中にシークしてもサムネイルが移動しない

これらの点が気に食わないので、
楽せずにやる方法で現在は実装中。
具体的にはC#でDirectShowを使ってます。
DirectShowについては今後まとめていく予定。

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まで)引用

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

結果出力

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

レベル1
f:id:s_sikisya:20140205171527p:plain
レベル2
f:id:s_sikisya:20140205171546p:plain

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

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