Haar関数系とWavelet分解

2013 Aug 18 Masahiro Mizutani

HaarWavelet[j, k, x]=(Subscript[ψ, j, k])[x]=ψ[(2^j)x-k]

定義

Haarwavelet[x]=ϕ[2x] -ϕ[2x-1]

HaarScale[x_Real]:=If[x<0 || 1 <x, 0,1]

Haarwavelet[x_Real]:=HaarScale[2x]-HaarScale[2x-1]

Plot[Haarwavelet[x],{x,-1,2}]

例

f1[t_]:=-2HaarScale[t+1]+2 HaarScale[t]+3HaarScale[t-1]+3HaarScale[t-2]-HaarScale[t-3]

一般化Haar Wavelet関数の定義

Mathematicaの省略可能な引数機能を使う。実用のために正規化係数 (2^(j/2)) は付けない

HaarScale[x, j, k]=(Subscript[ϕ, j, k])[x]=ϕ[(2^j)x-k] 
Haarwavelet[x, j, k] = ψ[(2^j)x-k] = (Subscript[ψ, j, k])[x]=ϕ[(2^(j+1))x-2k] -ϕ[(2^(J+1))x-2k-1]

一般化Haar sacaling関数を使った定義

HaarScale[x_Real, j_Integer,k_Integer]:=If[x<k/2^j|| (k+1)/2^j <x, 0,1]

Haarwavelet[x_Real, j_Integer, k_Integer]:=HaarScale[x,j+1,2k]-HaarScale[x,j+1,2k+1]

Plot[Haarwavelet[x,0,0],{x,-1,2}]

(Subscript[W, j])=(Subscript[{ψ[(2^j)x-k]}, k∈N])の直交性、(Subscript[およびW, k≠j])との直交性

Plot[Haarwavelet[x,0,2],{x,0.5,4}]

Plot[Haarwavelet[x,1,2],{x,0.5,2}]

Plot[Haarwavelet[x,1,1],{x,-0.5,2}]

Plot[Haarwavelet[x,2,4],{x,-0.5,2}]

spike niseのある関数例

spiked[t_]:=Cos[t^3]Sin[2t]+(Sin[110(t-2)]/(110(t-2)))^2-(Sin[95(t-1)]/(95(t-1)))^2

Plot[spiked[t],{t,0,Pi}]

DFTしてみる

coeff=Table[spiked[t],{t,0,Pi,Pi/255}];

ListPlot[InverseFourier[spikedDFT=Fourier[coeff]],Joined->True]

絶対値で小さいRatio割を零とする関数

compress[y_List,Ratio_Real]:=Module[{n=Length[y],nr,threshold},
nr=Floor[n Ratio];
threshold=Sort[Abs[y],Less][[nr]];
y /. x_ /; Abs[x] < threshold ->0
]

FDTの絶対値の小さい8割を0にしてみる

filteredSpikedDFT=compress[spikedDFT,0.8];

ListPlot[InverseFourier[filteredSpikedDFT],Joined->True]

Haar wavelet分解

f(x) on [a,b] の関数について、((これを2)^j)個の等区間に分割してx=((b-a)/((2^j)-1))kで標本化して、(Subscript[a, k])=f(((b-a)/((2^j)-1))k)k=0,...,(2^j)-1 
として、次の[0,1]上の階段関数(Subscript[f, j])(x)を定義する。

(Subscript[f, j])(x)=(∑_k)(Subscript[a, k])ϕ((2^j)x-k)=(∑_k)(((Subscript[a, 2k])-(Subscript[a, 2k+1]))/2)ψ((2^(j-1))x-k)+(((Subscript[a, 2k])+(Subscript[a, 2k+1]))/2)ϕ((2^(j-1))x-k)=(Subscript[w, j-1])(x)+(Subscript[f, j-1])(x)

この操作を繰り返して
(Subscript[f, j])(x)=(Subscript[w, j-1])(x)+(Subscript[w, j-2])(x)+...+(Subscript[w, j-k])(x)+(Subscript[f, j-k])(x)
における(Subscript[f, k])(x)をHaar k次近似関数。

(Subscript[f, j])(x)=(Subscript[w, j-1])(x)+(Subscript[w, j-2])(x)+...+ (Subscript[w, 0])(x)+(Subscript[f, 0])(x)
まで分解できる。これらは[0,1]上の関数である。(Subscript[f, 0])(x)は階段関数(Subscript[f, j])(x)全体の平均値。

初期分解能 ((b-a)/(2^8))

j=8;

L=2^j-1

関数の定義区間を2^j個の等間隔点でサンプルする

coeff=Table[N[spiked[t]],{t,0,Pi,Pi/(2^j-1)}];

Mathematicaの添字が1から始まることを考慮して

dspiked[x_Real]:=Sum[coeff[[k+1]]HaarScale[x,j,k],{k,0,2^j-1}]

Plot[dspiked[x],{x,0,1}]

具体例

cの奇数番目のリスト

oddCoeff=coeff[[1;;Length[coeff];;2]];

cの偶数番目のリスト

evenCoeff=coeff[[2;;Length[coeff];;2]];

ただしこれらの要素数は半分の2^(j - 1) 個

aveCoeff=(evenCoeff+oddCoeff)/2;

MathematicaのList番号が1からはじまるが故に、0 から添字が始まるときの[even] - [odd] は、Mathematicaでは[odd] - [even] となることに注意する

diffCoeff=(oddCoeff-evenCoeff)/2;

detailed7[x_Real]:=Sum[diffCoeff[[k+1]]Haarwavelet[x,7,k],{k,0,2^7-1}]

averaged7[x_Real]:=Sum[aveCoeff[[k+1]]HaarScale[x,7,k],{k,0,2^7-1}]

Plot[averaged7[x],{x,0,1}]

Plot[detailed7[x],{x,0,1}]

Plot[detailed7[x]+averaged7[x],{x,0,1}]

[a, b] 上の関数f (t) のk次Haar近似とk次Haar wavelet差分

区間[a,b]((を2)^n)に等分割してHaar scaling関数で階段関数(Subscript[f, j])(x)にしておく。k次Haar近似した関数(Subscript[f, j-k])(x)を返す

初期分解能 ((b-a)/(2^8)) のときのspiked関数の階段係数(Haar scaling係数）

coeff=Table[N[spiked[x]],{x,0,Pi,Pi/(2^8-1)}];

k次Haar scaling近似

averagedHaarScaling[x_Real,k_Integer]:=Module[{cf,n,oddc,evenc,avecf,j,i},
(* coef=Table[N[f[x]],{x,a,b,(b-a)/(2^n-1)}]; *)
cf=coeff;
n=Log[2,Length[coeff]];
For[j=0,j<k,j++,
oddc=cf[[1;;2^(n-j);;2]];
evenc=cf[[2;;2^(n-j);;2]];
avecf=(oddc+evenc)/2;
cf=avecf;
];
Sum[cf[[i+1]]HaarScale[x,n-k,i],{i,0,2^(n-k)-1}]
]

k次Haar wavelet差分

k次Haar Scaling分解 (Subscript[f, j])(x)=(Subscript[w, j-1])(x)+(Subscript[w, j-2])(x)+...+ (Subscript[w, j-k])(x)+(Subscript[f, j-k])(x)におけるHaar wavelet成分(Subscript[{(Subscript[w, j-i])(x)}, i=1..k])

HaarWaveletElement[x_Real,k_Integer]:=Module[{cf,n,oddc,evenc,avecf,diffcf,j,i},
(* coef=Table[N[f[x]],{x,a,b,(b-a)/(2^n-1)}]; *)
cf=coeff;
n=Log[2,Length[coeff]];
For[j=0,j<k,j++,
oddc=cf[[1;;2^(n-j);;2]];
evenc=cf[[2;;2^(n-j);;2]];
avecf=(oddc+evenc)/2;
cf=avecf;
];
diffcf=(oddc-evenc)/2;
Sum[diffcf[[i+1]]Haarwavelet[x,n-k,i],{i,0,2^(n-k)-1}]
]

k次Haar近似 = (k +1) 次Haar近似 + (k + 1) 次Haar wavelet差分

(Subscript[f, j-k])(x)=(Subscript[w, j-k-1])(x)+ (Subscript[f, j-k-1])(x)

(Subscript[f, 8])(x)=(Subscript[w, 7])(x)+(Subscript[w, 6])(x)+...+ (Subscript[w, 0])(x)+(Subscript[f, 0])(x)
以下では
(Subscript[w, j])=HaarWaveletElement[x,j-8], (Subscript[f, j])=averagedHaarScaling[x,j-8]

0 次近似(初期分解能 j=2^8のとき)

Plot[averagedHaarScaling[x,0],{x,0,1}]

Plot[HaarWaveletElement[x,1]+averagedHaarScaling[x,1],{x,0,1}]

Plot[HaarWaveletElement[x,1],{x,0,1}]

1 次近似

Plot[averagedHaarScaling[x,1],{x,0,1}]

Plot[HaarWaveletElement[x,2]+averagedHaarScaling[x,2],{x,0,1}]

Plot[HaarWaveletElement[x,2],{x,0,1}]

2 次近似

Plot[averagedHaarScaling[x,2],{x,0,1}]

Plot[HaarWaveletElement[x,3]+averagedHaarScaling[x,3],{x,0,1}]

Plot[HaarWaveletElement[x,3],{x,0,1}]

3次近似

Plot[averagedHaarScaling[x,3],{x,0,1}]

Plot[HaarWaveletElement[x,4]+averagedHaarScaling[x,4],{x,0,1}]

Plot[HaarWaveletElement[x,4],{x,0,1}]

4次近似

Plot[averagedHaarScaling[x,4],{x,0,1}]

Plot[HaarWaveletElement[x,5]+averagedHaarScaling[x,5],{x,0,1}]

Plot[HaarWaveletElement[x,5],{x,0,1}]

5次近似

Plot[averagedHaarScaling[x,5],{x,0,1}]

Plot[HaarWaveletElement[x,6]+averagedHaarScaling[x,6],{x,0,1}]

Plot[HaarWaveletElement[x,6],{x,0,1}]

6次近似

Plot[averagedHaarScaling[x,6],{x,0,1}]

Plot[HaarWaveletElement[x,7]+averagedHaarScaling[x,7],{x,0,1}]

Plot[HaarWaveletElement[x,7],{x,0,1}]

7次近似

Plot[averagedHaarScaling[x,7],{x,0,1}]

Plot[HaarWaveletElement[x,8]+averagedHaarScaling[x,8],{x,0,1}]

Plot[HaarWaveletElement[x,8],{x,0,1}]

8次近似

Plot[averagedHaarScaling[x,8],{x,0,1}]

Plot[{averagedHaarScaling[x,4],averagedHaarScaling[x,5]},{x,0,1}]

reconstruct[x_]=Sum[HaarWaveletElement[x,k],{k,1,8}]+averagedHaarScaling[x,8]

Plot[reconstruct[x],{x,0,1}]

reconstructchoice[x_]=Sum[HaarWaveletElement[x,k],{k,{1,2,3,4,5,6}}]+averagedHaarScaling[x,8]

Plot[reconstructchoice[x],{x,0,1}]