SlideBar

芦野 隆一・Remi Vaillancourt 著：

はやわかり MATLAB（共立出版）より

SlideBar

ウェーブレットツールボックス（Wavelet Toolbox）は MATLAB の数値計算とグラフィックスの機能を拡張し，局所的な解析， 多重スケール解析，あるいは非定常現象の解析に適したウェーブレット解析の数々の手法を包括的に提供するツールボックスである． ウェーブレット解析は非常に新しく進んだ手法であって， フーリエ解析が用いられてきたあらゆる分野で フーリエ解析を補って別の視点から解析を行う． ウェーブレット解析によって， 従来の手法では検知できなかったデータの傾向や崩壊点あるいはデータの導関数の不連続点や自己相似性といった様相を視覚化したり探究することできる．詳しくは拙著「ウェーブレット解析」を参照されたい．

特徴

このツールボックスの特徴は分析，合成，雑音除去，信号や画像の圧縮などのための完全なグラフィカルユーザインタフェースとコマンドラインから使うことができる関数を豊富にもつことである． 特にグラフィカルユーザインタフェースの完成度は高く， 多重スケールで信号解析ができる連続ウェーブレット変換， 1 次元および 2 次元の離散ウェーブレット変換， 信号や画像を解析する多重解像度分解， 1 次元および 2 次元のウェーブレットパケット展開， ウェーブレットパケットにおけるエントロピーにもとづく最適基底や最適レベル，各種の閾値による雑音除去などのウェーブレット解析をすべて対話的に行うことが可能である． したがって，コマンドラインからウェーブレット解析を行う必要はほとんどないが， このホームページでは MATLAB の数値計算機能を理解してもらうため，コマンドラインから行うウェーブレット解析の例を述べよう． 詳しくは拙著「はやわかり MATLAB」を参照されたい．

1 次元サンプルデータの作成

このツールボックスには幾つかのサンプルデータが用意されているが， ここでは解析するサンプルデータを準備しよう．

x = 0:4; y = zeros(1,5); y(3) = 2; t = linspace(0,4,2^9);
sample = interp1(x,y,t);

さらに save sample sample によりディスクに変数 sample を保存しておく． 次に wnoise を使って Quadchirp と呼ばれる振動数がしだいに高くなる雑音を付加したサンプルデータを作ろう．

noise = wnoise(5,9); nsample = sample+0.1*noise;

さらに save nsample nsample によりディスクに変数 nsample を保存しておく． どのようなデータか知りたければ

plot([nsample; sample]')

とすれば次の図を得る．

Samles Plot

コマンドラインからの1 次元連続ウェーブレット解析

help waveinfo により使えるウェーブレット関数を確認する． たとえば，ドベシィのウェーブレットを使うことにし， waveinfo('db') で詳しい説明を読むと db1 or haar, db4, db1 などというウェーブレットを指定する文字列 'wname' が与えられていることがわかる． 1 次元連続ウェーブレット変換の関数はコマンドと関数のリストによれば cwt であるから， help cwt で説明を読むとベクトルのデータ列 S， スケールベクトル SCALES， ウェーブレットを指定する文字列 'wname' に対し， COEFS = CWT(S,SCALES,'wname','plot') により 1 次元連続ウェーブレット変換の係数を計算し， 図示することがわかる． 例として，c = CWT(s,[3 18 12.9 7 1.5],'db2') などがある． サンプルデータ sample がワークスペースになければ load sample でディスクから読み込む．そして

c = cwt(sample,1:32,'db3','plot');

を実行してみると次の図（左）が得られる． カラーマップを colormap(pink) により pink に変えてみると次の図（右）が得られる．

CWT Plot

サンプルデータ sample の第 1 次導関数の不連続点が係数に反映されていることがわかる． 得られた係数は save coef c によりディスクに変数 c を coef というファイル名で保存しておくことができる．

コマンドラインからの1 次元離散ウェーブレット解析

1 次元離散ウェーブレット変換の関数はコマンドと関数のリストによれば

分析・分解

合成・再構成

分解構造の制御

雑音除去・圧縮

dwt

idwt

detcoef

ddencmp

wavedec

waverec

appcoef

wdencmp

wrcoef

upwlev

wden

upcoef

プロンプトで load nsample とタイプし， サンプルデータ nsample をディスクから読み込む． ドベシィの N=3 のウェーブレットを使うことにする． レベル 1 の 1 次元離散ウェーブレット変換係数を計算するには

[cA1,cD1] = dwt(nsample,'db3');

である．この係数からレベル 1 の近似（approximation）を構成するには

ls = length(nsample); A1 = upcoef('a',cA1,'db3',1,ls);

詳細（detail）を構成するには

D1 = upcoef('d',cD1,'db3',1,ls);

である．この結果を図示するには

subplot(1,2,1); plot(A1); subplot(1,2,2); plot(D1);

とする．次の図を得る．

DWT Level 1 Plot

1 次元逆離散ウェーブレット変換は

A0 = idwt(cA1,cD1,'db3',ls);

である． レベル 3 まで 1 次元ウェーブレット分解を計算するには

[C,L] = wavedec(nsample,3,'db3');

である．レベル 3 の近似係数は

cA3 = appcoef(C,L,'db3',3);

レベル 3 からレベル 1 まで詳細係数は

cD3 = detcoef(C,L,3); cD2 = detcoef(C,L,2);
cD1 = detcoef(C,L,1);

である．レベル 3 の近似を構成するには

A3 = wrcoef('a',C,L,'db3',3);

レベル 1 からレベル 3 までの詳細を構成するには

D1 = wrcoef('d',C,L,'db3',1); D2 = wrcoef('d',C,L,'db3',2);
D3 = wrcoef('d',C,L,'db3',3);

である．この結果を図示するには

subplot(2,2,1); plot(A3); subplot(2,2,2); plot(D1);
subplot(2,2,3); plot(D2); subplot(2,2,4); plot(D3);

とする．次の図を得る．

DWT Level 3 Plot

付加した雑音 noise = wnoise(5,9) は Quadchirp と呼ばれる振動数がしだいに高くなる雑音である． レベル 3 の分解によりレベルが高くなるに従ってしだいに低周波部分が 分離されているのがわかる． レベル 3 の分解から元の信号を再構成するには

A0 = waverec(C,L,'db3');

である．元の信号とレベル 3 の分解の近似を

subplot(1,2,1); plot(nsample); v = axis;

subplot(1,2,2); plot(A3); axis(v);

によりプロットする．次の図を得る．

DWT Level 3 Comparing Plot

図の左右を比較すれば， 一般に詳細を捨てることにより信号から高周波の雑音を除去することができて，除去した結果の信号が近似として与えられることがわかる． 関数 ddencmp により閾値を与えてより細かい雑音除去ができる． 正数 t を閾値とするとき， hard thresholding は信号 x が abs(x) <= t を満たす成分は 0 とし， soft thresholding は信号 x が abs(x) <= t を満たす成分は 0 とし，かつ abs(x) > t を満たす成分は sign(x).*(abs(x)-t) とする． 詳しい使い方は help ddencmp によりオンラインヘルプを参照すればよい．

SlideBar
Hits COUNTER (Since Jan. 4, 2001)