圧縮アルゴリズム

(8) ウェーブレット変換 -1-

この章では、JPEG2000 フォーマットで圧縮アルゴリズムとして採用された「ウェーブレット変換 ( Wavelet Transform )」について取り上げたいと思います。
JPEG2000 では画像圧縮アルゴリズムとして採用されましたが、元々はフーリエ変換の応用形としてデータ解析の利用を目的に考案された手法であり、その内容も多岐に渡っています。全てを説明するのは、量的にも自分の理解度から見ても不可能なため、ここでは、離散ウェーブレット変換とそれを用いた解析手法である多重解像度解析を中心に紹介をしたいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) フーリエ変換とウェーブレット変換

データや関数の特性を分析する手法としてよく知られたものに「フーリエ変換 ( Foueier Transform )」があります (*1-1)。フーリエ変換は、周期性を持った任意の関数を正弦波と余弦波の和で表すことができることを利用してデータを周波数成分に変換するために利用される重要な解析手法です。フーリエ変換とその逆変換は以下のような式で表されます。

F(ω) = ( 1 / √2π )∫{-∞→∞} f(t)・e-iωt dt

f(t) = ( 1 / √2π )∫{-∞→∞} F(ω)・eiωt

上側の式がフーリエ変換で、時間や変位の関数 f(t) が周波数成分の関数 F(ω) に変換されています。また、下側の式は逆フーリエ変換で、F(ω) から f(t) へ戻す場合に利用します。この中には正弦波・余弦波が見当たりませんが、それは e-iωt の中に隠されています。オイラーの定理から

e±iωt = cos ωt ± isin ωt

となるので、

f(t)・e-iωt = f(t)cos ωt - i・f(t)sin ωt

F(ω)・eiωt = F(ω)cos ωt + i・F(ω)sin ωt

と表すこともできます。

この変換法ではデータ全体を周波数成分に変換するため時系列変化や変位などの情報が欠落してしまい、特性が時間や場所などによって変化するようなデータを解析する場合には利用することができません。これを解決する方法として「窓フーリエ変換 ( Window Fourier Transform )」が考案されました。この手法は「短時間フーリエ変換 ( Short-time Fourier Transform )」とも呼ばれ、有限の区域でのみ 0 以外の値を持った「窓関数 ( Window Function )」と呼ばれる関数を分析したいデータと積算してからフーリエ変換することによって、局所的なデータの周波数成分を求めることになります。窓関数を軸上で移動させることで、周波数成分の変化を解析することもできます。このときの変換式は、窓関数の位置を表す変数と周波数を表す変数の、二つの変数を持った関数となります。
窓関数は、ある幅の中でのみ値を持ち、その他の場所では 0 になるような関数であれば、どのようなものでも利用できます。データの一部を切り取ってフーリエ変換を行った場合は、窓関数として「矩形窓 (Rectangular Window)」を利用したことと同じであることは容易に理解できると思います。また、その他にも「ハン窓 (Hann Window ; またはハニング窓 Hanning Window )」「ハミング窓 (Hamming Window)」「ブラックマン窓 (Blackman Window)」などがあります。有限な幅を持っていない関数に対しては、途中で打ち切ってゼロにするか、他の有限な幅を持った関数と掛け合わせて使います。例えば「ガウス窓 (Gaussian Window)」はガウス関数を使うため定義域は実数全体になり端の部分は強制的にゼロにする必要があります。
窓関数の幅はたいてい任意に決めることができます。この幅を小さくすれば軸上の変化を敏感にとらえることができますが、扱うデータ量が少なくなるため周波数成分の精度は悪くなります。逆に幅を大きくすると周波数成分の精度は上がりますが、軸上の変化が追いにくくなってしまいます。このような関係は「不確定性の原理」といい、窓フーリエ変換の大きな欠点になります。

窓フーリエ変換の中で理論上最良の変換方法としては、窓関数にガウス関数を用いた「ガボール変換(Gabor Transform)」がよく知られています。1975 年頃に、フランスの石油探査技師「ジーン・モルレー ( Jean Morlet )」は、石油探査のために地中への振動波の反射をガボール変換によって解析し、油床の場所を特定することを試みていましたが、その解析は非常に困難な作業でした。そのため 1980 年代はじめ頃、三角関数の重ね合わせによって反射信号を表現する代わりに、Wavelet と呼ばれる短い波を拡大縮小・平行移動して得られる波の重ね合わせによって表現することを考案しました。これがウェーブレット変換の始まりです。

三角関数は無限に連続した関数であるのに対し、Wavelet は局所的に値が存在する関数です。フーリエ変換は種々の周波数成分を持った三角関数の重ね合わせによって波を表現しているため、時間や場所などに対して常に一定のパターンを持ったデータの解析には有用ですが、時刻・場所によってパターンの変化するデータを解析するなどの用途には不向きと言えます。しかし、ウェーブレット変換では局所的な波を平行移動したり、拡大縮小したものを使って波を表現するため、有限の区間内にあるデータの特性を解析するには三角関数より適しています。


対象の関数に窓関数 g(t) を積算した結果を使ってフーリエ変換を行ったときの変換式は次のようになります。

F( ω ) = ( 1 / √2π )∫{-∞→∞} g( t )・f( t )・e-iωt dt

窓関数を軸上で移動させることによって、任意の位置における局所的な周波数特性を解析することができます。この位置を変数 b として加えたものが下に示した窓フーリエ変換の式になります。

F( b, ω ) = ( 1 / √2π )∫{-∞→∞} g( t - b )・f( t )・e-iωt dt

ウェーブレット変換では、三角関数の代わりに「マザー・ウェーブレット (Mother Wavelet)」と呼ばれる基本参照波 ψ(t) を拡大縮小・平行移動した参照波

ψa,b( t ) = ( 1 / √a )ψ( ( t - b ) / a )

を使って変換処理を行います。変換式は下のようになり、これを「連続ウェーブレット変換 ( Continuous Wavelet Transformation; CWT )」と呼びます。

W( b, a ) = ( 1 / √a )∫ f( t )・ψa,b( t ) dt

マザー・ウェーブレットとして利用される基本参照波は通常 ∫{-∞→∞} ψ(t) dt = 0 を満たしています。積分した結果が 0 になるので、信号の平坦な ( 定数に近い ) 位置に対して積算と積分を行った場合、その値はほぼ 0 になります。またウェーブレット ψa,b(t) が関数 f(t) によく似ていれば積分の値は大きくなりますが、そうでなければこの性質から正負の値が打ち消すようになり、積分の結果はやはり 0 に近い値となることが直感的に理解できると思います。
ψ(t) の幅はパラメータ a に対応して a 倍になるため、1 / a は周波数に対応していることになります。つまり a の値が大きい場合は幅が広がるため低周波に、逆に小さい場合は幅が狭くなって高周波に反応します。


窓フーリエ変換と連続ウェーブレット変換のサンプル・プログラムを以下に示します。

/* ----------------------------------------------------------------------------
   Domain : 定義域クラス
   ---------------------------------------------------------------------------- */
struct Domain
{
  double min;  // 最小値
  double max;  // 最大値
  size_t reso; // 解像度

  // コンストラクタ : 定義域と解像度を指定して構築
  Domain( double d0, double d1, size_t r )
    : min( ( d0 < d1 ) ? d0 : d1 ),
      max( ( d0 > d1 ) ? d0 : d1 ),
      reso( r ) {}

  // range : 定義域の幅を返す
  double range() const { return( max - min ); }
};

/* ----------------------------------------------------------------------------
   WindowFunction : 窓関数クラス
   ---------------------------------------------------------------------------- */
template< class Res >
class WindowFunction
{
  vector< Res > data_; // 計算結果
  Domain domain_;      // 計算する幅(値を持つ範囲)

public:

  // コンストラクタ : 値を持つ範囲 domain と関数 f を指定して構築
  template< class F >
  WindowFunction( Domain domain, F f )
    : domain_( domain )
  {
    double inc = range() / reso(); // 増分
    double current = domain_.min;  // 計算対象を最小値で初期化

    for ( size_t i = 0 ; i < reso() ; ++i ) {
      data_.push_back( f( current ) );
      current += inc;
    }
  }

  // operator[] : 計算結果を返す
  // 添字 i を引数としてあらかじめ計算した結果を返す
  Res operator[]( size_t i ) const
  { return( data_.at( i ) ); }

  // range : 値を持つ範囲の幅を返す
  double range() const
  { return( domain_.range() ); }

  // reso : 解像度(範囲内の計算値の個数)を返す
  size_t reso() const
  { return( domain_.reso ); }
};

/* ----------------------------------------------------------------------------
   窓フーリエ変換
   ---------------------------------------------------------------------------- */

/*
  WinFourierTrans : f で計算されるデータを、窓関数 win を使って窓フーリエ変換する

  a : 周波数成分
  b : 位置成分
  win : 窓関数
  f : データ関数
*/
template< class F >
double WinFourierTrans( double a, double b, const WindowFunction< double >& win, F f )
{
  assert( a >= 0 );

  complex< double > res; // 変換結果

  double t = b - win.range() / 2;        // 計算開始位置 b + t
  double inc = win.range() / win.reso(); // t の増分

  for ( size_t i = 0 ; i < win.reso() ; i++ ) {
    // win(b+t)・f(t)・exp(ia(b+t)) = win(t')・f(t'-b)・exp(iat')
    double d = f( t ) * win[i];
    res += complex< double >( d * cos( a * t ), -d * sin( a * t ) );
    t += inc;
  }

  return( sqrt( std::norm( res ) / ( 2 * M_PI ) ) * inc );
}

/*
  WinFourierMain : 各周波数・位置成分毎に窓フーリエ変換を行う

  win : 窓関数
  f : データ関数
  res : 計算結果を保持する配列へのポインタ
  scaleDomain : 周波数成分の定義域と解像度
  transDomain : 位置成分の定義域と解像度
*/
template< class F >
void WinFourierMain( const WindowFunction< double >& win, F f, map< double, map< double, double > >* res, const Domain& scaleDomain, const Domain& transDomain )
{
  double incScale = scaleDomain.range() / scaleDomain.reso; // 周波数成分の増分
  double incTrans = transDomain.range() / transDomain.reso; // 位置成分の増分

  res->clear();
  double a = scaleDomain.min; // 現在の周波数成分を最小値で初期化
  for ( size_t ai = 0 ; ai < scaleDomain.reso ; ++ai ) {
    (*res)[a] = map< double, double >();
    double b = transDomain.min; // 現在の位置成分を最小値で初期化
    for ( size_t bi = 0 ; bi < transDomain.reso ; ++bi ) {
      (*res)[a][b] = WinFourierTrans( a, b, win, f );
      b += incTrans;
    }
    a += incScale;
  }
}

/* ----------------------------------------------------------------------------
   連続ウェーブレット変換
   ---------------------------------------------------------------------------- */

/*
  CalcRes : 積算結果 res を整形して返す

  結果 = 積算結果 x 増分 / √a
  但し、積算結果が複素数の場合はノルムにしてから計算する

  res : 積算結果
  inc : x の増分
  a : スケーリング成分
*/
double CalcRes( double res, double inc, double a )
{
  return( res * inc / sqrt( a ) );
}
double CalcRes( complex< double > res, double inc, double a )
{
  return( sqrt( std::norm( res ) / a ) * inc );
}

/*
  WaveletTrans : 連続ウェーブレット変換

  a : スケーリング成分
  b : 位置成分
  win : 窓関数
  f : データ関数
*/
template< class T, class F >
double WaveletTrans( double a, double b, const WindowFunction< T >& win, F f )
{
  assert( a >= 0 );

  if ( a == 0 ) return( 0 );

  T res = 0; // 変換結果

  double t = b - win.range() * a / 2;  // 計算位置 b + at
  double inc = win.range() * a / win.reso(); // t の増分

  for ( size_t i = 0 ; i < win.reso() ; i++ ) {
    // f(b+at)・win(t) = f(t')・win((t'-b)/a)
    res += f( t ) * win[i];
    t += inc;
  }

  return( CalcRes( res, inc, a ) );
}

/*
  WaveletMain : 各スケール・位置 パラメータ毎に連続ウェーブレット変換を行う

  win : 窓関数
  f : データ関数
  res : 計算結果を保持する配列へのポインタ
  scaleDomain : スケーリング成分の定義域と解像度
  transDomain : 位置成分の定義域と解像度
*/
template< class T, class F >
void WaveletMain( const WindowFunction< T >& win, F f, map< double, map< double, double > >* res, const Domain& scaleDomain, const Domain& transDomain )
{
  double incScale = scaleDomain.range() / scaleDomain.reso; // スケーリング成分の増分
  double incTrans = transDomain.range() / transDomain.reso; // 位置成分の増分

  res->clear();
  double a = scaleDomain.min; // 現在のスケーリング成分を最小値で初期化
  for ( size_t ai = 0 ; ai < scaleDomain.reso ; ++ai ) {
    (*res)[a] = map< double, double >();
    double b = transDomain.min; // 現在の位置成分を最小値で初期化
    for ( size_t bi = 0 ; bi < transDomain.reso ; ++bi ) {
      double d = WaveletTrans( a, b, win, f );
      (*res)[a][b] = d;
      b += incTrans;
    }
    a += incScale;
  }
}

窓フーリエ変換と連続ウェーブレット変換に利用する窓関数は、クラス WindowFunction で定義するようにしています。利用する関数を指定して構築をしますが、関数 f の型はテンプレート引数で渡されます。このため、通常の関数と関数オブジェクトのどちらでも使用可能です。その他に定義域 Domain インスタンスを渡し、値を保持する範囲を指定します ( Domain は定義域の最小・最大値と解像度 ( 定義域内のサンプリング数 ) を保持するクラスです )。なお、窓関数の返す値は周波数や位置に依存せず常に一定となるようにしているので、値は構築時に先に計算して配列に保持しておきます。
ところで、関数をコンストラクタの引数として渡すのではなく、関数そのものを純粋仮想関数として定義して、各窓関数の実装は派生クラスとして用意する方法もあります。しかし、構築時にデータをあらかじめ用意するので、その時には関数の内容が実装済みでなければなりません。コンストラクタから派生クラスの仮想関数を呼び出すことはできないので、ここでは引数として渡す方法を採用しています。また、計算結果の型はテンプレート引数 Res で決定します。これは、窓関数に実数・複素数それぞれを返すものが混在しているための処置です。

窓フーリエ変換で使用する関数は WinFoueierTransWinFoueierMain の二つで、WinFoueierTransWinFoueierMain から呼び出されて実行するようになっています。計算処理は WinFoueierTrans が行い、指定された周波数成分と位置成分を使ってデータ・窓関数・三角関数を積算し、その和を求めて積分値としています。データはサンプリングされたものであるため、積分は解析的に実行するのではなく、個々に計算した結果の和 ( 矩形の面積 ) を計算してその和を求める手法をとっていて、これは「区分求積法 ( Quadrature by Parts )」の中の最も単純な処理方法です。なお、積算した結果で和を計算して、最後に増分を掛け合わせることで面積を求めています。
前述の通り、窓関数として渡された win はすでに計算した結果を保持しています。保持したデータは win が持っている定義域 domain_ で決まり、配列の最初の要素が最小値 domain_.min の位置のデータとなります。関数 f と三角関数を計算するときに使う変数 t は、位置成分 b から定義域の幅の半分だけ引いた値で初期化されます。増分 inc を順次加えながら積算・和の計算を win.reso() 回行い、最終的には定義域の最大値まで計算を行うので、結果的には

F( b, a ) = ( inc / √2π )Σk{0→win.reso()-1}( f( b - win.range() / 2 + inc * k )・win[i]・e-ia( b - win.range() / 2 + inc * k ) )

を求めていることになり、窓関数となる配列 win の要素は

win[i] = ψ( win.domain_.min + inc * k )

を意味することになります。窓関数の最初の位置を b から幅の半分だけ負の方向へ移動した位置に合わせ、b が窓関数の中央と一致するように調整していると考えるとわかりやすいと思います。

窓フーリエ変換の場合は窓関数そのもので幅が決まります。a は周波数成分を表し、三角関数だけに作用することになります。また、窓フーリエ変換の場合は窓関数として実数を返す関数だけを利用可能としています。三角関数は複素数を含むため計算結果は複素数ですが、返り値はそのノルム ( 複素数 a + bi に対して ( a2 + b2 )1/2 ) としているので実数値が返されます。

連続ウェーブレット変換は WaveletTransWaveletMain で行い、 WaveletTransWaveletMain から呼び出されて計算処理を行なっています。積分の方法や窓関数 win の使い方は窓フーリエ変換の場合と同じであり、以下のような計算処理が行われています。

W( b, a ) = ( inc / √a )Σk{0→win.reso()-1}( f( b - win.range() * a / 2 + inc * k )・win[i] )

win[i] = ψ( win.domain_.min + inc * k )

窓フーリエ変換との大きな違いは、窓関数自体が変数 a によって拡大・縮小されているという点です。プログラム内では関数 f の方を a 倍していますが、これは窓関数を 1 / a していることと同等です。よって、a が大きくなるほど窓関数の幅は小さくなり、より高周波側に対して反応するようになります。なお、f を計算する開始位置は b から窓関数の幅の半分だけ負の方向へシフトしたところで、これは窓フーリエ変換と同じ手順となっています。
また、連続ウェーブレット変換では窓関数に複素数を返す関数を利用することができるようにしてあります。そのため、WaveletTrans の内容が少しだけ複雑になっていますが、前述の通り窓関数は実数と複素数のいずれかを返すため、その返り値の型をテンプレート引数 T としています。しかし、最終的に返す値はどちらも実数なので ( 複素数の場合はノルムを返す )、実数用と複素数用の返り値計算関数 CalcRes を用意して型に応じて切り替わるようにしています。最後の返り値計算以外は全く同じ実装になるため、このような形をとっています。


以下に示すようなデータを使い、サンプル・プログラムによりいくつかの条件で窓フーリエ変換を行ってみます。このデータは、正弦波 sin( πt ) と sin( πt / 2 ) を合成した波に対して、t = π / 2 と t = 7π / 4 付近に大きなピークを付加した擬似的なものです。解析によって、合成された波の周波数とピークの位置を同時に検出したいとします。

テストデータ
図 1-1. テストデータ

まずは窓関数として「ガウス窓 (Gauss Window)」を使います。これは正規分布などでおなじみのガウス関数です。

g( t ) = ( 1 / σ )exp( -t2 / 2σ2 )

正規分布は確率分布なので、全区間での積分値が 1 になるように定数の係数 1 / ( 2π )1/2 が付いていますが、ここでは正規化することは意味がないため省略しています。分散を意味する σ2 が窓の幅を調整するためのパラメータとなり、値が大きく(小さく)なるほど幅も大きく(小さく)なります。ガウス窓を使った窓フーリエ変換は「ガボール変換 (Gabor Transform)」と呼ばれます。

σ = 0.1σ = 1σ = 5
窓フーリエ変換による処理結果 (σ=0.1) 窓フーリエ変換による処理結果 (σ=1) 窓フーリエ変換による処理結果 (σ=5)
図 1-2. 窓フーリエ変換による処理結果 (窓関数 : ガウス窓)

上図において白い部分ほど値が大きく、逆に黒い部分はゼロに近い個所を表しています。なお、値はノルムなので負の数はありません。横軸が位置 b を表し、縦軸が周波数 a で下に向かうほど値は大きくなります。σ = 0.1 は最も幅の小さな窓関数であり、二本のピークに対応した縦方向の白いすじがはっきりと確認できるのに対して二つの異なる周波数の波の合成であることはこの結果から読み取ることができません。σ = 1 では縦横どちらの方向もぼやけた状態になってしまい、かろうじて情報が読み取れるという状態となっています。σ = 5 のときは横方向の白い線が確認できることから二つの周波数成分を持った波があることがわかりますが、ピークに対する位置情報については何の結果も得られません。
このように、窓フーリエ変換では周波数成分と位置成分の両方を同時にとらえることが困難です。窓関数の幅をを大きくすれば周波数成分をとらえやすくなる反面、局在したデータの情報は他のデータによってぼやけてしまいます。逆に幅を小さくすると、局在したデータの情報はとらえやすくなりますが、全体の波の挙動についてはとらえきれなくなります。

ガウス関数を窓関数に使ったフーリエ変換は以下のような式になります。

F( b, ω ) = [ 1 / ( 2π )1/2σ ]∫{-∞→∞} exp( -( t - b )2 / 2σ2 )・f( t )・e-iωt dt

この式は、f(t) にガウス関数と正弦波の積を作用させる形となっています。ガウス関数と正弦波の積は「ガボール・フィルタ (Gabor Filter)」と呼ばれ、帯域フィルタなどにも利用されている関数です。先に紹介したジーン・モルレーは、このガボール変換の代わりに次のような式を使った変換法を考案しました。

ψ( t ) = ( 1 / σ )exp( -t2 / σ2 )[ exp( iπt ) - exp( -π2σ2 / 4 ) ]

これを基本参照波とした変換が連続ウェーブレット変換の始まりとなります。この基本参照波は「モルレー・ウェーブレット ( Morlet Wavelet )」と呼ばれます。また、ガボール変換に登場するガウス関数と正弦波の積は「ガボール・ウェーブレット ( Gabor Wavelet )」と呼ばれますが、本質的には大きな差はありません。以下に、両者で連続ウェーブレット変換を行った結果を示します。窓フーリエ変換の場合とは異なり、二つのピークの位置と周波数成分の帯を読み取ることが可能となります。

Morlet wavelet ; σ = 2Gabor wavelet ; σ = 2
連続ウェーブレット変換による処理結果 (Morlet) 連続ウェーブレット変換による処理結果 (Gabor)
図 1-3. 連続ウェーブレット変換による処理結果

モルレー・ウェーブレットは以下のような形状をしています。基本的にはガウス関数と正弦波の重ねあわせにより作られているので、ガボール・ウェーブレットもほぼ同じ形状です。また、複素関数であり、実数部と虚数部は位相がずれただけで波形が非常によく似ています。ノルムを計算した結果をみるとガウス関数と同じ釣り鐘のような形状をしていることも下図から理解できると思います。

Morlet Wavelet
図 1-4. モルレー・ウェーブレット ( σ = 2 )

ウェーブレットには他にも様々な種類があります。その中から二つほど紹介しておきます。

Haar
wavelet処理結果 (R=0.3)
Haar waveletHaar による処理結果 (R=0.3)
Mexican Hat
wavelet処理結果 (σ=0.25)
MexicanHat waveletMexicanHat による処理結果 (σ=0.25)

上記ウェーブレットは以下のような式になります。これらはモルレー・ウェーブレットやガボール・ウェーブレットとは異なり返り値が実数になります。

■ Haar wavelet
ψ(t)=1[ 0 ≤ t < R ]
=0[ t < 0 R ≤ t ]
■ MexicanHat wavelet
ψ(t) = [ 1 / ( 3σ )1/2 ]( 1 - t2 / σ2 )・exp( t2 / 2σ2 )

フーリエ変換には、連続値を扱うもの以外に、離散値を扱う「離散フーリエ変換 ( Discrete Fourier Transform ; DFT )」というものがありました。離散フーリエ変換は以下のような式で表されます。

Fk = (1/N)Σl{0→N-1}( fle-i2πkl/N )

fl = (1/N)Σk{0→N-1}( Fkei2πkl/N )

ウェーブレット変換にも離散値を扱うための「離散ウェーブレット変換 ( Discrete Wavelet Transform ; DWT )」があり、マザー・ウェーブレット ψ(t) を拡大縮小・平行移動した参照波は以下のように表されます。

ψj,k(t) = 2-j/2ψ( 2-jt - k )

j, k はそれぞれ 連続ウェーブレットでのスケーリング・位置成分だった a, b に対応する正数で、a には 2j、b には 2j・k が対応します。

離散ウェーブレット変換の場合、パラメータ j が増加するごとにスケーリング成分が 2 倍ずつ大きくなります。また、位置成分は j と k を変数として持っているので、スケーリング成分が大きくなるほど移動量も増加することになります。この特徴をを利用して、次に説明する多重解像度解析を行うことができます。

*1-1) フーリエ変換については、前回紹介した「JPEG法」の他、「グラフィック・パターンの扱い (5) サンプル補間」の中の「5) 補間関数と畳み込み積分」や「数値演算法 (4) 高速フーリエ変換」などでも紹介しています。

2) 多重解像度解析 ( Multi-Resolution Analysis ; MRA )

多重解像度解析とは、画像や音声などのパターンを異なる周波数成分に分解する作業を何度も行うことによって、その特徴を解析する手法のことです。パターンを分解する際に使用される関数を「スケーリング関数 ( Scaling Function )」と呼び、入力された関数はこのスケーリング関数の和で表されることになります。ここでは、最も単純なスケーリング関数であるハール ( Haar ) のスケーリング関数を例にして説明したいと思います。

Haar のスケーリング関数は以下のように定義される矩形波です ( 前節のハール・ウェーブレットとは形が異なることに注意して下さい )。

φ(t)=1[ 0 ≤ t < 1 ]
=0[ t < 0 ; 1 ≤ t ]

このスケーリング関数を使い、任意の関数 f(t) を以下のように関数 f(0)(t) へ近似することができます。

f(0)( t ) = Σk( skφ( t - k ) )

但し sk = ∫{-∞→∞} f( t )φ( t - k ) dt = ∫{k→k+1} f(t) dt

式で書くとややこしそうに見えますが、要は k から k + 1 までの f( t ) の面積の矩形に置き換わったと考えればわかりやすいと思います。連続しているデータも、測定などはサンプリングで行われるため、そのようなデータは近似された f(0)(t) の形をとります。

Haarのスケーリング関数とそれにより近似された関数
図 2-1. Haar のスケーリング関数とそれにより近似された関数

このスケーリング関数 φ(t) を基本参照波と見なしたとき、これを拡大縮小・平行移動した関数は以下のように表すことができます。

φk(j)( t ) = 2-j/2φ( 2-jt - k )

関数 f(t) は φ(j)k(t)を使って近似することも可能で、この時の近似関数を f(j)(t) としたとき

f(j)( t ) = Σk( sk(j)・φk(j)( t ) )

と表現することができます。このときの sk(j) を「スケーリング係数 ( Scaling Coefficients )」と呼びます。

j の値が 1 増える毎に、φ(t)の幅は 2 倍ずつ大きくなっていきます。よって近似関数は粗くなり、f(j+1)(t) は f(j)(t) に比べて情報が欠落していくことになります。この欠落分を g(j+1)(t) としたとき、次の式が成り立ちます。

f(j)( t ) = f(j+1)( t ) + g(j+1)( t )

このとき、gj+1(t) をレベル j + 1 のウェーブレット成分と呼びます。

Haar のスケーリング関数で近似された関数 f(j+1)(t) は、f(j)(t) 上の 2 つの矩形を平均化してひとつの矩形にしたものを並べた形になります。よって g(j+1)(t) は、以下のような関数で表される基本参照波を用いて表すことができます。この関数 ψ(t) を Haar の「ウェーブレット関数 ( Wavelet Function )」といいます。こちらが、前節の連続ウェーブレット変換で紹介した関数と同じものです。

ψ( t )=1[ 0 ≤ t < 1 / 2 ]
=-1[ 1 / 2 ≤ t < 1 ]
=0[ 0 < t ; t ≤ 1 ]
 
Haarのウェーブレット関数及び近似関数と誤差の関係図
図 2-2. Haar のウェーブレット関数及び近似関数と誤差の関係図

Haar のウェーブレット関数も、スケーリング関数と同様に拡大縮小・平行移動した関数として以下のように表すことができます。

ψk(j)( t ) = 2-j/2ψ( 2-jt - k )

また、ウェーブレット成分は下式のように表すことができます。

g(j)( t ) = Σkωk(j)ψk(j)( t )

このとき、ωk(j) をレベル j の「ウェーブレット係数 ( Wavelet Coefficients )」と呼びます。

f(0)(t) は f(1)(t) と g(1)(t)の和で表され、f(1)(t) はさらに f(2)(t) と g(2)(t) の和で表されます。これを繰り返すと、関数 f(0)(t) を以下のように表すことができます。

f0(t)=g(1)(t) + g(2)(t) + ... + g(j)(t) + f(j)(t)
=Σi{1→j}( g(i)(t) ) + f(j)(t)

このように、関数を幅の異なる ( つまり解像度の異なる ) ウェーブレット成分の和で表現して解析する方法を「多重解像度分析 ( Multi-Resolution Analysis ; MRA )」といいます。


式だけではわかりづらいので、例としてHaarのスケーリング関数を使って関数 y = x ( 0 ≦ x ≦ 4 ) を表現してみたいと思います。

表 2-1. 多重解像度解析の例
f(0)( t ) = (1/2)φ( t ) + (3/2)φ( t - 1 ) + (5/2)φ( t - 2 ) + (7/2)φ( t - 3 ) f(x)のスケーリング関数による表現
φk(1)( t ) = (1/√2)φ( t / 2 - k ) より
 f(1)( t )=φ( t / 2 ) + 3φ( t / 2 - 1 )
=√2φ0(1)( t ) + 3√2φ1(1)( t )
s0(1) = √2, s1(1) = 3√2
ψk(1)( t ) = (1/√2)ψ( t / 2 - k )より
 g(1)( t )=-0.5ψ( t / 2 ) - 0.5ψ( t / 2 - 1 )
=-(√2/2)ψ0(1)( t ) - (√2/2)ψ1(1)( t )
ω0(1) = -√2 / 2, ω1(1) = - √2 / 2
f(0)( t )=g(1)( t ) + f(1)( t )
=-(√2/2)ψ0(1)( t ) - (√2/2)ψ1(1)( t ) + √2φ0(1)( t ) + 3√2φ1(1)( t )
=-(1/2)ψ( t / 2 ) - (1/2)ψ( t / 2 - 1 ) + φ( t / 2 ) + 3φ( t / 2 - 1 )
レベル1の分解
φk(2)( t ) = (1/2)φ( t / 4 - k ) より
 f(2)( t )=2φ( t / 4 )      
=0(2)( t )
s0(2) = 4
ψk(2)( t ) = (1/2)ψ( t / 4 - k ) より
g(2)( t )=-ψ( t / 4 )      
=-2ψ0(2)( t )
ω0(2) = -2
f(1)( t )=g(2)( t ) + f(2)( t )
=-2ψ0(2)( t ) + 4φ0(2)( t )
=-ψ( t / 4 ) + 2φ( t / 4 )
f(0)( t )=g(1)( t ) + g(2)( t ) + f(2)( t )
=-(√2/2)ψ0(1)( t ) - (√2/2)ψ1(1)( t ) -2ψ0(2)( t ) + 4φ0(2)( t )
=-(1/2)ψ( t / 2 ) - (1/2)ψ( t / 2 - 1 ) - ψ( t / 4 ) + 2φ( t / 4 )
レベル2の分解

レベル j からレベル j + 1 への変換をする場合、レベル j の近似データに対して隣りあったデータの平均値を取り、それを 2(j+1)/2 倍するとスケーリング係数が得られます。ウェーブレット係数は、平均値と元のデータの誤差 ( 二つのうちどちらを使っても、符号が反転するだけで絶対値は等しくなります ) を 2(j+1)/2 倍することで求めることができます。あるレベルまで計算を行う場合、レベル 0 から順に同じ処理を繰り返すことになります。
一回の処理によって、n 個のデータが n / 2 個のスケーリング係数と n / 2 個のウェーブレット係数に分割されるため、データ量としては変化しないことになります。しかし、JPEG の章でも紹介したように、自然画像は近隣の画素間の変化量が小さい傾向があるため、上記のような近似を行ったときの誤差 g(j)(t) はそれほど大きな値を持たないことが予想できます。近似を行った結果、情報量は半分になり、誤差成分は小さな値となるので圧縮が可能になります。逆に近隣の画素間に大きな変化があった場合、その変化は g(j)(t) に強く表れるため、多重解像度解析によりエッジ検出などを行うことも可能です。
なお、画像は二次元平面上にあるため、ウェーブレット変換を使う場合は通常、水平方向と垂直方向の二回に分けて行う必要があります。また、逆変換を行う場合、平均値を表す s(j) とその誤差を表す ω(j) との和と差をそれぞれ計算する事で行う事ができます。


Haar のスケーリング関数とウェーブレット関数を使った変換処理のサンプル・プログラムを以下に示します。

/*
  WaveletCoeff : スケーリング・ウェーブレット係数成分
*/
struct WaveletCoeff
{
  vector< double > s;  // スケーリング係数計算結果
  vector< double > wv; // ウェーブレット展開係数計算結果(垂直方向成分)
  vector< double > wh; // ウェーブレット展開係数計算結果(水平方向成分)
  vector< double > wd; // ウェーブレット展開係数計算結果(対角方向成分)

  // 配列のサイズを初期化(リサイズ)する
  void resize( Coord< int > size )
  {
    s.resize( size.x * size.y );
    wv.resize( size.x * size.y );
    wh.resize( size.x * size.y );
    wd.resize( size.x * size.y );
  }
};

/*
  Wavelet_Haar_1D : 一次元方向のWavelet変換

  tp : 入力信号 ( レベルがひとつ低いスケーリング係数 ) へのポインタ
  sp : スケーリング係数計算結果へのポインタ
  wp : ウェーブレット展開係数計算結果へのポインタ
  size : 信号のデータ長
  inc : indexの増分(X/Y方向を示す)
*/
void Wavelet_Haar_1D( const double* tp, double* sp, double* wp, size_t size, size_t inc )
{
  for ( size_t i = 0 ; i < size / 2 ; ++i ) {
    /*
      s = ( t0 + t1 ) / √2
      w = √2t0 - s = ( t0 - t1 ) / √2
    */
    *sp = ( *tp + tp[inc] ) / sqrt( 2 );
    *wp = *tp * sqrt( 2 ) - *sp;
    /*
      s = ( t0 + t1 ) / 2
      w = t0 - s = ( t0 - t1 ) / 2
    */
    //*sp = ( *tp + tp[inc] ) / 2;
    //*wp = *tp - *sp;
    tp += inc * 2;
    sp += inc;
    wp += inc;
  }
}

/*
  Wavelet_Haar_2D : 二次元方向のWavelet変換

  yuv : YUV成分
  wc : スケーリング・ウェーブレット係数成分へのポインタ
  width : 画像の幅
  height : 画像の高さ
*/
void Wavelet_Haar_2D( const vector< double >& yuv, WaveletCoeff* wc, Coord< int > size )
{
  size_t halfWidth = size.x / 2;
  vector< double > sx( halfWidth * size.y );
  vector< double > wx( halfWidth * size.y );

  wc->resize( Coord< int >( halfWidth, size.y / 2 ) );

  // YUV 成分を水平方向に変換 → 水平方向のスケーリングとウェーブレット
  for ( size_t i = 0 ; i < static_cast< size_t >( size.y ) ; i++ )
    Wavelet_Haar_1D( &yuv[i * size.x], &sx[i * halfWidth], &wx[i * halfWidth], size.x, 1 );

  // 水平方向のスケーリング・ウェーブレット成分をそれぞれ垂直方向に変換
  // 水平方向のスケーリング → 両方向のスケーリングと縦方向のウェーブレット
  // 水平方向のウェーブレット → 水平方向のウェーブレットと対角方向のウェーブレット
  for ( size_t i = 0 ; i < halfWidth ; i++ ) {
    Wavelet_Haar_1D( &sx[i], &( wc->s )[i], &( wc->wv )[i], size.y, halfWidth );
    Wavelet_Haar_1D( &wx[i], &( wc->wh )[i], &( wc->wd )[i], size.y, halfWidth );
  }
}

/*
  IWavelet_Haar_1D : 一次元方向のWavelet逆変換

  sp : スケーリング係数計算結果へのポインタ
  wp : ウェーブレット展開係数計算結果へのポインタ
  tp : 逆変換結果 ( レベルがひとつ低いスケーリング係数 ) へのポインタ
  size : 信号のデータ長
  inc : indexの増分(X/Y方向を示す)
*/
void IWavelet_Haar_1D( const double* sp, const double* wp, double* tp, size_t size, size_t inc )
{
  for ( size_t i = 0 ; i < size ; i++ ) {
    /* t0 = ( s + w ) / √2 */
    *tp = ( *sp + *wp ) / sqrt( 2 );
    /* t0 = s + w */
    //*tp = *sp + *wp;
    tp += inc;
    /* t1 = ( s - w ) / √2 */
    *tp = ( *sp - *wp ) / sqrt( 2 );
    /* t1 = s - w */
    //*tp = *sp - *wp;
    tp += inc;
    sp += inc;
    wp += inc;
  }
}

/*
  IWavelet_Haar_2D : 二次元方向のWavelet逆変換

  WaveletCoeff* wc : 入力信号
  double* yuv : 逆変換結果(レベルがひとつ低いスケーリング係数)
  unsigned int width : 画像の幅
  unsigned int height : 画像の高さ
*/
void IWavelet_Haar_2D( const WaveletCoeff& wc, vector< double >* yuv, Coord< int > size )
{
  Coord< size_t > half( size.x / 2, size.y / 2 );
  vector< double > sx( half.x * size.y );
  vector< double > wx( half.x * size.y );

  yuv->resize( size.x * size.y );

  // 両方向のスケーリング + 縦方向のウェーブレット → 水平方向のスケーリング
  // 水平方向のウェーブレット + 対角方向のウェーブレット → 水平方向のウェーブレット
  for ( size_t i = 0 ; i < half.x ; i++ ) {
    IWavelet_Haar_1D( &( wc.s )[i], &( wc.wv )[i], &sx[i], half.y, half.x );
    IWavelet_Haar_1D( &( wc.wh)[i], &( wc.wd )[i], &wx[i], half.y, half.x );
  }

  // 水平方向のスケーリングとウェーブレット → YUV 成分
  for ( size_t i = 0 ; i < static_cast< size_t >( size.y ) ; i++ )
    IWavelet_Haar_1D( &sx[i * half.x], &wx[i * half.x], &(*yuv)[i * size.x], half.x, 1 );
}

入力信号は double 型の変数で表されています。画素は通常 RGB ですが、変換時に RGB から YUV への変換を行うことを前提としているためこのようになっています。JPEG の場合と同様に、YUV 各成分ごとに、変換後の圧縮処理に対する圧縮率を細かく制御することが期待できますが、ウェーブレット変換処理だけに限れば RGB 成分をそのまま変換することも有用で、RGB 各成分のエッジ抽出などに利用可能です。

四つの関数のうち、実際に変換処理と逆変換処理を行っているのはそれぞれ Wavelet_Haar_1DIWavelet_Haar_1D で、残りの Wavelet_Haar_2DIWavelet_Haar_2D は二次元方向へ拡張するために Wavelet_Haar_1DIWavelet_Haar_1D を x, y 方向に対して二回呼び出しているだけの関数です。

変換処理は、次の二式で行われています。

*sp = ( *tp + tp[inc] ) / sqrt( 2 );

*wp = *tp * sqrt( 2 ) - *sp;

また、逆変換処理は

*yuv = ( *s1 + *w1 ) / sqrt( 2 );
 :
*yuv = ( *s1 - *w1 ) / sqrt( 2 );

が行います。一回の処理でレベルが一つ上がり、これを何度も呼び出すことによって、任意のレベルまで処理を行うことができます。処理内容は、単純に二つのピクセルの平均とその誤差を求めるだけなので、JPEG 法での処理に比べて処理自体は非常にシンプルで、高速な処理も期待することができます。なお、変換処理と逆変換処理の両方に係数を等分するため √2 で割るという処理が含まれていますが、本来の意味をそのままコーディングする場合はコメントアウトされた方を利用することで実現ができます ( こちらの方がよりシンプルで、実際の処理の意味もとらえやすくなります ) (*2-1)。

サンプル・プログラムを使い、画像を YUV 変換した上でウェーブレット変換処理を一回だけ実行した結果を以下に示します。テストに使用したのは、画像処理テスト用としてよく利用されることで有名な "Lenna" です。なお、実際には 512 x 512 のサイズをテストに使用しており、下に示した画像は縮尺されています。

図 2-3. テスト画像 (Lenna)
Lenna
 
図 2-4. 変換結果
ウェーブレット変換結果

変換結果の中で、左上側はスケーリング係数 ( 低周波成分 )、その他の三つはウェーブレット展開係数 ( 高周波成分 )で、右上側が水平成分、左下側が垂直成分、そして右下側が対角成分をそれぞれ表しています。次のレベルの変換を行う場合は、左上側のスケーリング係数に対して同じ処理を繰り返すことになります。
少し補足しておくと、スケーリング係数は得られたデータを 1 / 2 にしてから RGB 変換しています。スケーリング係数を計算する時に二つのデータの和の 1 / 2 ではなく 1 / √2 倍しているので水平・垂直方向に処理した時に実際の平均値よりも 2 倍大きな値になっています。そのため、描画前に 1 / 2 倍して補正を行う必要があります。逆にウェーブレット展開係数は値自体が小さいので、計算値を 50 倍してから RGB へ変換しています。

ウェーブレット係数は水平・垂直・対角それぞれの方向に対して画像のエッジ部分が強くなっていることがこの結果からわかります。スケーリング係数にさらに同様の処理をすれば、少し離れたピクセルどうしでエッジ抽出をすることになり、これを繰り返すことで各周波数成分に対するエッジが検出できる様子が直感的にも理解できると思います。また、実際のウェーブレット係数の値はほとんどが小さく、この部分が圧縮可能な個所になります。
なお、このサンプル・プログラムでは画像のスケールを単純に半分 ( 小数点部分を切り捨て ) にしているので、画像サイズが奇数であるようなときは端の部分が切り捨てられていきます。画像のサイズが 2 のべき乗で表されるのなら何度処理しても問題はありませんが、通常はそのようなことはないので本来ならサイズが奇数の場合の処理を加える必要があります。やり方としては、仮想のピクセルがあると仮定して処理するか、スケーリング係数を対象のピクセルの値そのままに、ウェーブレット係数はゼロにする方法が考えられます。後者は、仮想のピクセルが端のピクセルと全く同じ色であると仮定した場合に相当します。


*2-1) スケーリング・ウェーブレット関数の定義を考慮せずに処理の中身を見れば、処理の内容は隣り合ったデータの平均と、各データと平均との誤差を求めているだけの単純なものです。和を √2 で割る代わりに、平均を計算するため 2 で割っても問題はなく、サンプル・プログラムの中でも示した通り処理は可能です。切り替えた場合は、逆変換もスケーリング係数 ( 平均 ) にウェーブレット係数 ( 誤差 ) を加算・減算するだけという非常にシンプルな処理になります。
√2 で割る処理にしてあるのは、定義したスケーリング・ウェーブレット関数に合わせるという意味もありますが、こうすることでウェーブレット変換・逆変換のどちらにも 1 / √2 という同じ係数が付くことになります。実はフーリエ変換でも同様に二種類の定義があり、本章で紹介した計算式以外にも

F(ω) = ∫{-∞→∞} f(t)・e-iωt dt

f(t) = ( 1 / 2π )∫{-∞→∞} F(ω)・eiωt

という表し方があります。実際、参考にした書籍の中ではこの式を採用したものもあり、本章で紹介した表し方は数学者に多いそうです。さらに、ω は角周波数で単位時間あたりの角度を意味するので、代わりに単位時間あたりの振動数 f を利用すれば

ω = 2πf

より

F(f) = ∫{-∞→∞} f(t)・e-i2πft dt

f(t) = ∫{-∞→∞} F(f)・ei2πft df

となって係数は表れなくなります。こちらは電気工学関連の書籍に多いようです。

本章で表現したフーリエ変換に合わせてウェーブレット変換の式を定義しているだけなので、どちらがいいかは気にする必要はありません。しかし、実装に関してはコメントアウトした側の方が単純で処理後のデータも扱いやすいので、実際に利用する場合は切り替えた方がよさそうです。


Haar のスケーリング関数とウェーブレット関数を使って、多重解像度解析により二種類の成分に分解する手法について説明してきました。この手法は考え方も実装方法も非常に単純であり、処理も高速になるのでかなり有用な方法になりますが、スケーリング関数としてはこの他にもいくつかの種類があります。次の章でもウェーブレット変換について取り上げ、Haar の関数に代わるものとして Daubechies の関数などを紹介したいと思います。また、変換後の量子化や圧縮処理についても説明をする予定です。


<参考文献>
  1. C MAGAZINE 2003年1月号 「アルゴリズムラボ 〜 ウェーブレット変換による圧縮処理(1) ウェーブレット変換の基礎」
  2. 「これなら分かる応用数学教室」 (金谷健一 著, 共立出版)
  3. 「ウェーブレット解析の基礎理論」 (新井康平 著, 森北出版)
  4. 「ウェーブレット10講」 (Ingrid Daubechies 著, シュプリンガー・ジャパン)
  5. PIT Lab短時間フーリエ変換と連続ウェーブレット変換」 「離散ウェーブレット変換について
  6. 大阪教育大学 Ashino's Homepageウェーブレットの歴史
  7. サイバネット システム 株式会社 「Wavelet Toolboxの紹介と利用の開始 (リンク切れ)」
  8. Wikipedia

◆◇◆更新履歴◆◇◆

◎ 図 2-4. 変換結果を差し替えました。ウェーブレット変換結果の内容がほとんどわからない状態でした (2016-05-24)


[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system