高速Haar Wavelet変換と逆Wavelet変換
Fourier変換との比較

Masahiro MIZUTANI
Nov 27, Nov 22 2013
29 Nov 2012,  11 Nov  2011,  22 Dec 2011

データ数が2のスケール J 乗 L=(2^J)の標本リスト signal={(Subscript[c, 0]^J), (Subscript[c, 1]^J), ..., (Subscript[c, L-1]^J)} を与えて、J回のWavelet変換を繰り返して得られる係数リストを高速に計算する。

高速Haar wavelet変換

[高速Wavelet変換アルゴリズム] Mathematicaの配列指定は 1 から始まることを考慮している：signal=

L=(2^J)の標本データ signal=((Subscript[c, 0]^J), (Subscript[c, 1]^J), (Subscript[c, 2]^J),...., (Subscript[c, L-1]^J)) を与えて、J回wavelet変換を繰り返して得られるリスト ((Subscript[c, 0]^((J))), (Subscript[d, 0]^((J-1))) , (Subscript[d, (Subscript[k, 2])]^((2))),...., (Subscript[d, (Subscript[k, i])]^((n)))) を返す
 
L=(2^J)
data = signal 
for j = J to 1
	for k = 0 to (2^(j-1))-1
 		tmp = data[(2^(J+1-j))k + 1]
 		data[(2^(J+1-j))k + 1]= (1/2) (tmp + data[(2^(J+1-j))k + (2^(J-j)) + 1])
  		data[(2^(J+1-j))k + (2^(J-j)) + 1] = (1/2) (tmp - data[(2^(J+1-j))k + (2^(J-j)) + 1])
dataを返り値とする

FirstHaarWaveletTransform[signal_List]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j≥1,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
data
]

u={3,1,0,4,8,6,9,9};

haartransf=FirstHaarWaveletTransform[u]

wavelet変換レベル毎のスケール係数

関数 f のスケール J の標本デー(Subscript[タf, J])（(2^J)(Subscript[個の関数値c, k]^J)）に対し、Haar近似（(Subscript[c, k]^j)  から (Subscript[c, k]^(j-1)) を求めること）を J-level 回（J≥level≥0）を繰り返して得られる近似関数 (Subscript[f, level])の係数リスト (Subscript[c, k]^level)を求める（デー((タ数は2)^J)から (2^level)個に減少）

HaarScalingCoefficient[signal_List,level_Integer]:=Module[{scale,tmp,data,j,k,rec={}},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j>level,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
For[k=0, k≤2^level-1,k++,
AppendTo[rec,data[[2^(scale-level)  k+1]]];
];
rec
]

HaarScalingCoefficient[u,3]

HaarScalingCoefficient[u,2]

HaarScalingCoefficient[u,1]

HaarScalingCoefficient[u,0]

w={3,1,0,4,8,6,9,9,4,6,8,2,1,5,9,3};

Length[w]

FirstHaarWaveletTransform[w]

Haar wavelet逆変換

[高速逆Wavelet変換]
L=(2^J)の標本データ signal=((Subscript[c, 0]^J), (Subscript[c, 1]^J), (Subscript[c, 2]^J),...., (Subscript[c, L-1]^J)) をJ回wavelet変換を繰り返して得られるリスト ((Subscript[c, 0]^((J))), (Subscript[d, 0]^((J-1))) , (Subscript[d, (Subscript[k, 2])]^((2))),...., (Subscript[d, (Subscript[k, i])]^((n)))) が与えられたとして、逆変換を繰り返して「元データ」を返す。 Mathematicaの配列指定は 1 から始まることを考慮している：

L=(2^J)
data = waveletdata
for j = 0 to J-1
	 for k = 0 to  (2^j)-1
 		tmp = data[(2^(J-j))k + 1]
 		data[(2^(J-j))k + 1] =tmp + data[(2^(J-j))k + (2^(J-j-1))+1]
  		data[(2^(J-j))k + (2^(J-j-1))+1] = tmp - data[(2^(J-j))k+ (2^(J-j-1))+1]
dataを返り値とする

FirstHaarInverseWaveletTransform[waveletdata_List]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[waveletdata]];
data=waveletdata;
For[j=0,j<scale,j++,
For[k=0,k<=2^j -1,k++,
tmp=data[[2^(scale-j) k+1]];
data[[2^(scale-j) k+1]]=tmp + data[[2^(scale-j) k+2^(scale-j-1)+1]];
data[[2^(scale-j) k+2^(scale-j -1)+1]]=tmp-data[[2^(scale-j) k+2^(scale-j-1)+1]];
]
];
data
]

FirstHaarInverseWaveletTransform[FirstHaarWaveletTransform[u]]

計算時間の測定

長さ (2^J)のリスト datalist に対して、高速Harr Wavelet変換を使うとどの程度計算時間がかかるかを計測してみる。

(2^maxpower)の長さの乱数リストを生成

FirstHaarWaveletTransform[datalist=Table[RandomReal[],{2^4}]]

AbsoluteTiming[FirstHaarWaveletTransform[datalist]]

First[AbsoluteTiming[FirstHaarWaveletTransform[datalist]]]

(2^1)((から2)^maxpower)の長さのHaar wavelet変換の計算時間を計測する

TimingHaarWaveletTransf[maxpower_Integer]:=Module[{timelist={},rlonglist,i,rlist,n,usedtime},
rlonglist=Table[RandomReal[],{2^maxpower}];
For[n=1,n≤ maxpower,n++,
rlist=Take[rlonglist,2^n];
usedtime=First[AbsoluteTiming[FirstHaarWaveletTransform[rlist]]];
AppendTo[timelist,usedtime];
];
timelist
]

高速Wavelet変換/Wavelet逆変換を使うと、(2^18)=262144個程度のリスト長のWavelet変換が可能となる。
(2^18)個は、CDのサンプリング周波数44100Hzのモノラルデータで換算すると、(2^18)/44100=5.94秒分に相当する。

timinglist=TimingHaarWaveletTransf[18]

ListPlot[timinglist,Filling->Axis]

レベル指定のwavelet変換

本来は、スケール J=scale からJ-1,..., 0 まで変換を繰り返す

レベルpまでのWavelet変換

((長さを2)^scale)のスケールscaleの信号 (Subscript[f, scale])（信号係数{(Subscript[c, k]^scale)}）をレベル scaleとする。レベル scale≥p≥0((を指定して得られる2)^scale)個のHaar scaling/wavelet係数のリスト{(Subscript[c, k]^j),(Subscript[d, k]^j)}を求める。

HaarWaveletTransfLevel[signal_List,p_Integer]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j>p,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
data
]

u={3,1,0,4,8,6,9,9};

HaarWaveletTransfLevel[u,3]

HaarWaveletTransfLevel[u,2]

HaarWaveletTransfLevel[u,1]

v={4,6,10,12,8,6,5,5};

HaarWaveletTransfLevel[v,0]

レベルpからの逆wavelet変換

スケール scale=J ((の信号の長さを2)^scale)としたとき、元信号{(Subscript[c, k]^scale)}をレベル scale として,レベル scale≥p≥0((までの2)^scale)個のHaar scaling/wavelet係数のリスト{(Subscript[c, k]^j),(Subscript[d, k]^j)}から、元信号{(Subscript[c, k]^scale)}を求める

InvHaarWaveletTransfLevel[waveletleveldata_List,p_Integer]:=Module[{data,scale,j,k,tmp},
scale=Log[2,Length[waveletleveldata]];
data=waveletleveldata;
For[j=p,j<scale,j++,
For[k=0,k<=2^j -1,k++,
tmp=data[[2^(scale-j) k+1]];
data[[2^(scale-j) k+1]]=tmp + data[[2^(scale-j) k+2^(scale-j-1)+1]];
data[[2^(scale-j) k+2^(scale-j -1)+1]]=tmp-data[[2^(scale-j) k+2^(scale-j-1)+1]];
]
];
data
]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,0],0]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,1],1]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,2],2]

InvHaarWaveletTransfLevel[HaarWaveletTransfLevel[u,3],3]

Level pのHaar近似値

HaarScalingCoefficient[signal_List, level_Integer]で得られる。

レベルpのwavelet変換係数を並べる

スケー((ルscaleの長さを2)^scale)の信号 signal={(Subscript[c, k]^scale)} が与えられている。
まず、高速wavelet変換を使ってレベル p (0≤p≤scale)までの (2^scale)個のHaar scaling/wavelet係数のリストを求める。
このリストから、
	リストの前半（左側）にレベルp ((の2)^p)個の平均化データ {(Subscript[c, k]^p)}
	リストの後半（右側）に平均化データを取り除いた（順序は変えない）詳細データ  {{(Subscript[d, k]^p)},{(Subscript[d, k]^(p+1))},..,{(Subscript[d, k]^(scale-1))}}
を並べる。

リスト lst からその部分リスト rmlst を取り除く removeList[lst, rmlst]

リスト lst から、部分リスト（lstに含まれる取り除きたい要素の集まりからなるリスト） rmlst を取り除く

removeList[lst_List,rmlst_List]:=Delete[lst,Flatten[Map[Position[lst,#]&,rmlst],1]]

removeList[{a,b,c,d,e,f},{b,d,f}]

Flatten[Map[Position[{a,b,c,d,e,f},#]&,{a,b}],1]

Delete[{a,b,c,d,e,f},Flatten[Map[Position[{a,b,c,d,e,f},#]&,{a,b}],1]]

arrangeHaarWaveletTransfList[signal_List,p_Integer]:=Module[{data,scale,j,k,tmp,wdata={}},
removeList[lst_List,rmlst_List]:=Delete[lst,Flatten[Map[Position[lst,#]&,rmlst],1]];
scale=Log[2,Length[signal]];
data=signal;
For[j=scale,j>p,j--,
For[k=0,k≤2^(j-1) -1,k++,
tmp=data[[2^{scale+1-j} k+1]];
data[[2^{scale+1-j} k+1]]=(tmp+data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
data[[2^{scale+1-j} k+2^(scale-j)+1]]=(tmp-data[[2^{scale+1-j} k +2^(scale-j)+1]])/2;
]
];
For[k=0, k≤2^p-1,k++,
AppendTo[wdata,data[[2^(scale-p)  k+1]]];
];
Join[wdata,removeList[data,wdata]]
]

arrangeHaarWaveletTransfList[u,3]

arrangeHaarWaveletTransfList[u,2]

arrangeHaarWaveletTransfList[u,1]

arrangeHaarWaveletTransfList[u,0]

応用例

[ListHaarWaveletTransfLevel[signal, p]のプロットの左側はレベルp の平均化データで描かれていることに注意
後のFourie変換による場合に比べて、平均からのズレを表す右側データの絶対値が「小さい」ことにも注意

Plot[Exp[-x]Sin[10x],{x,0,2Pi},PlotRange->All]

J = 10

signal=Table[Exp[-m]Sin[10m],{m,0,2Pi,2Pi/1023}]//N;

Log[2,Length[signal]]

ListPlot[signal,PlotRange->All]

level 9

ListPlot[HaarScalingCoefficient[signal,9],PlotRange->All]

ListPlot[arrangeHaarWaveletTransfList[signal,9],PlotRange->All]

level 8

ListPlot[HaarScalingCoefficient[signal,8],PlotRange->All]

ListPlot[arrangeHaarWaveletTransfList[signal,8],PlotRange->All]

level 7

ListPlot[HaarScalingCoefficient[signal,7],PlotRange->All]

ListPlot[arrangeHaarWaveletTransfList[signal,7],PlotRange->All]

Fourier変換との比較

スケール J=9 周期L=512

J=9;
L=2^J

Fourier基底を設定

k = 0, ..., L-1

fbase[k_Integer,L_Integer]:=Table[Exp[2 Pi I k *j/L],{j,0,L-1}]

内積の定義

innerproduct[u_List,v_List]:=Module[{i,len=Length[u]},
Sum[u[[j]]Conjugate[v[[j]]],{j,1,len}]
]

離散Fourrier変換

discreteFourierTrans[u_List]:=Table[innerproduct[u,fbase[j,L]],{j,0,L-1}]

Sample

Riemannの稠密な点で微分不可能な連続関数のL点近似

riemann=Table[Sum[Sin[2*Pi*k^2*(j/L-0.2)]/k^2,{k,1,20}],{j,1,L}]//N;

ListPlot[riemann,Filling->Axis]

ListPlot[HaarScalingCoefficient[riemann,8],PlotRange->All]

ListPlot[ListHaarWaveletTransfLevel[riemann,8],PlotRange->All]

次のFourie変換による場合に比べて、Haar Wavelet変換では平均からのズレを表すデータの絶対値が「小さい」 ことにも注意

friemann=discreteFourierTrans[riemann];

ListPlot[Abs[friemann]^2,Filling->Axis]