ウェーブレット逆変換やっとこさ成功
疲れた．原因は基本数列のインデックスの入れ替え間違い．言い訳じゃないですけど，もらった文書に詳しく書いてなかったので試行錯誤を繰り返しました．
ソースは見せるほどもないものなので自粛．もうっちょと進化したら・・・．んーでもこのブログを日記として残しておくなら，ソースコードも書いておきたいところですね．紛失する可能性がかなり高いので．
そんことはまあいいとして，自分の理解度確認のためにウェーブレットについての説明．以下でウェーブレット変換っていうのは離散ウェーブレット変換のことです，念のため．
ウェーブレットについて何となく分かるためには，単なる数列を考えると良いです．ウェーブレット変換の対象となる入力は，以下のようにyで与えられます．
y ---> (1.5, 1.2, 1.8, 0.9, 0.5, 1.1, 1.6, 1.4)
今のところ入力として受け付けるのは，1次元のデータで，しかもデータ長が2の累乗になってるものだけです．実際は，2の累乗に足りない部分を0で埋めればいいでしょう．するってーと音楽ファイルみたいのも変換できそうです．
ウェーブレット変換は，このyを数学的な計算でいじくって，どんどん分解していこうってもんです．分解した後に再構成することも可能で，元のデータを取り戻すことができます．
ウェーブレット変換にはウェーブレット関数が必要です．難しい話を抜きにすると，これも関数じゃなくて数列だと思えば簡単です．以下の4つの数列が基本数列になります．
pp ---> (0.683013, 1.18301, 0.316987, -0.183013)
qq ---> (-0.183013, -0.316987, 1.18301, -0.683013)
hh ---> (0.482963, 0.836516, 0.224144, -0.12941)
gg ---> (-0.12941, -0.224144, 0.836516, -0.482963)
これはDaubechies2（ドブシー2階）というウェーブレットに関する基本数列です．他にも色々な数列がありますが，普通はそれらを自由に選ぶことができます．この数列がウェーブレット変換の重要な部分で（私の考えですが），こういった数列は簡単には得られません．上の数列はドブシーっていう偉い人が考えてくれたものです．ドブシー2階は長さが4の数列なんです．これはもう決まってることなので，ただその式をプログラムに書けばいいだけです．
既に書いたように，基本数列は上の4つなんですが，今回のプログラムではhhとggしか使っていませんので，これ以降はそれらについてだけ書いていきます．上の4つの数列は互いに深く関連しているので，1つ求めればインデックスをずらしたり，ルート2をかけたり割ったりすると他の数列も求めることができます．そんなわけで，本当は上の数列はインデックスが-1から始まってたりするんですが，その辺の意味はまだよく分かってません．とりあえずそのことは置いておいて．
さて，上の数列が基本になるのですが，実際は入力の長さと同じ長さの数列が必要になります．今，入力の長さが2^3=8だとすると，長さが8の数列が必要になります．長さ4の基本数列の，足りない部分を0で埋めて，長さ8の数列を作ります．
h ---> (0, 0.482963, 0.836516, 0.224144, -0.12941, 0, 0, 0)
g ---> (-0.224144, 0.836516, -0.482963, 0, 0, 0, 0, -0.12941)
これが実際に変換で使用する数列です．インデックスがちょっと変わってますが，これも決まってることです．
ここからが問題です．厄介です．イメージだけで話すと，ウェーブレット変換するっていうと，普通は何段階かの変換を行います．1段階の変換で，データの長さが半分になります．最終的に長さが1のデータを得るまでウェーブレット変換を行うので，入力のデータ長が2^3だと，3回のウェーブレット変換を行うことになります．1回のウェーブレット変換で，入力の長さが半分になってしまうので，先ほどのhとgもそれに合わせて半分にしてやらなくてはなりません．
そこで，上のhとgをそれぞれh0，g0として，h1,h2,h3とg1,g2,g3という数列を作ることになります．これらの数列の作り方は非常に簡単で，長さiの数列のちょうど半分で列をぶったぎって，長さi/2の2つの数列を作り，それらを重ね合わせることで長さがi/2の数列が1つ得られます．先ほどのh（h0）を例にすると，h1は，(0 + -0.12941, 0.482963 + 0, 0.836516 + 0, 0.224144 + 0)となります．
同様にして全てのhiとgiを作ります．
h[0] ---> (0, 0.482963, 0.836516, 0.224144, -0.12941, 0, 0, 0)
h[1] ---> (-0.12941, 0.482963, 0.836516, 0.224144)
h[2] ---> (0.707107, 0.707107)
h[3] ---> (1.41421)
g[0] ---> (-0.224144, 0.836516, -0.482963, 0, 0, 0, 0, -0.12941)
g[1] ---> (-0.224144, 0.836516, -0.482963, -0.12941)
g[2] ---> (-0.707107, 0.707107)
g[3] ---> (5.42101e-20)
ここでhとgが正しく作られたかを確かめるのに，各hiとgiの数列の和をとる方法があります．やってみるとわかりますが，hiの各要素の合計は全て1.41421くらいで，giの方はほぼ0になります．1.41421はルート2なのですが，必ず，sum(hi)=ルート2，sum(gi)=0となります．こうなっていなかったら，数列の生成に失敗したということです．
hとgを全部求めたら，やっとウェーブレット変換を行うことができます．ウェーブレット変換自体もイメージで説明すると，数列c0をウェーブレット変換すると，長さがc0の半分の2つの数列，c1とd1が得られます．ciは低周波成分，diは高周波成分となるのですが，簡単に言えば，元の数列の特徴をよく表したものがciになります．それ以外のどうでもいいっちゃどうでもいい部分がdiになります．そうやって出来たc1をウェーブレット変換して，c2とd2が得られます．あとはciのデータ長が1になるまで（もちろんその手前まででもおｋ）ウェーブレット変換を繰り返します．ここで，ciは捨てて構いませんが，diはとっておきます．再構成のときに必要になるからです．
ウェーブレット変換自体の計算は非常に単純で，プログラムだと2重の繰り返し文と掛け算で書けます．私の書いたプログラムだと，インデント込みで30行弱です．マジで簡単です．そんな簡単なのにすごいことになります．アドレナリンが出ます．うほほ．
ここでは便宜上，c0=yとします．yは入力でした．さきほどのhとgを総動員してウェーブレット変換を行っていくと，c0から以下のような数列がどんどん得られます．c[0]からc[1]とd[1]，c[1]からc[2]とd[2]って感じです．1つ変換が進むごとに，データ長が半分になっているのが分かると思います．
c[0] ---> (1.5, 1.2, 1.8, 0.9, 0.5, 1.1, 1.6, 1.4)
c[1] ---> (2.22231, 0.892428, 1.98937, 1.96696)
d[1] ---> (-0.382903, -0.0473672, -0.0811132, -0.0543022)c[2] ---> (2.24845, 2.75155)
d[2] ---> (-0.966923, 0.0107051)c[3] ---> (3.53553)
d[3] ---> (0.355752)
ここでちょっとc0に関して補足しておくと，本来なら，c0はyの補間係数の列となります．そのため，c0を求めるのに補間係数を求める計算が必要になるのですが，これがちょっと面倒です．そこで，「c0=yとしちゃえ」っていう考えがあります．これをマラーの便法というそうです．本気でやるなら補間係数を求める計算が必要ですが（本気っていやオレはいつでも本気なわけだが），この程度のカスプログラムならマラーさんにお世話になります．
あっさりしてますが，ウェーブレット変換はこれでお終いです．上のcとdのうち，dをいじくったりして圧縮できます．先ほども書いたように，dはあんまり意味ないデータだと考えることができるので，極端な話，dを全部捨ててしまっても，ウェーブレット逆変換をすると，元の入力データに割りと近いデータが得られたりします．
最後に再構成，つまりウェーブレット逆変換を行います．ウェーブレット逆変換は，そのままウェーブレット変換の逆を行います．cとdで説明すると，ciとdiからci-1を作るのがウェーブレット逆変換です．計算式は，こっちも簡単です．あっさり結果を出すと次のようになります．
out[0] ---> (2.24845, 2.75155)
out[1] ---> (2.22231, 0.892428, 1.98937, 1.96696)
out[2] ---> (1.5, 1.2, 1.8, 0.9, 0.5, 1.1, 1.6, 1.4)
プログラム的な理由で，outの添え字が進むほど，逆変換が進んでいることになっています．したがって，c[3]とd[3]からout[0]（＝c[2]）が得られ，out[0]とd[2]からout[1]（＝c[1]）が得られ，out[2]とd[1]からout[0]，つまりc[0]，つまりy，つまり入力が出てきます．すごいっすね，ほんと．なんでこうなるのかっていうのは分かってません．
ちなみに，何度か書いたように，dを全部0にしてから再構成すると以下のようになりました．
result ---> (1.38795, 1.51417, 1.60446, 0.929127, 0.458942, 1.16083, 1.54865, 1.39587)
yと結構近いことが分かります．それだけですけど．
私自身，ここまで来るのにそれなりに苦労してしまいましたが，これからは結構楽しそうです．見通しが大分良くなったので，次はとりあえず，音楽とかもっと大規模の周波数データを解析して，できれば早めに2次元に拡張したいですね．ICPCと院試が邪魔なんですけど．
