ガボールウェーブレット変換

Gabor wavelet transform

ガボールウェーブレット変換のC++コーディングのメモです。

ガボールウェーブレットとウェーブレット変換は、下図に示します。ψ(t)はガボールウェーブレット関数で、f(t)はデジタル信号で、WT(a,b)は変換結果の2次元配列で、ｔは時間で、ｂはシフトで、ａはスケールで、σは減衰係数です。

テストとして、下図のように、1000点の信号データを作ります。0.2秒毎に、周波数を0,8,16,24,0の順で、1秒間の正弦波のデータ作成します。

周波数を０、２，４から３0までの１６個の周波数で、時間を０ｍｓから９９９ｍｓまでの１０００点の時間で、ウェーブレット変換を行って、2次元配列に保存します。

ソースコード：

//定数定義、変数宣言
#define pi2 (8.0 * atan(1.0))   //2π
#define N 1000                  //データ点数                    
#define M 16                    //周波数の個数
#define FS 1000.0               //サンプリングの周波数(HZ)
#define sigma 1                 //係数のσ
#define limit 0.01              //ガボール減衰の最小値
double data[N];                 //信号データ
double result_r[M][N];          //変換結果データ[周は周][時間]実数部
double result_i[M][N];          //変換結果データ[周は周][時間]虚数部
// ψ(t) = 1 / (2 sqrt(π) σ) * exp(-t^2/σ^2) * exp(jωt)
//ガボールウェーブレット変換
void GaborTransform()
{
    //１６個の周波数のループ
    for (int m = 0; m < M; m ++)
    {
        double f = 2.0 * m;    //周波数=0..31HZ
        double a = 1.0 / f;    //スケールのａ＝周期
        //データ点数分までのループ
        for (int p = 0; p < N; p ++)
        {
            double valr = 0.0;
            double vali = 0.0;
            //畳み込み積分のループ
            for (int n = -p ; n < N - p; n ++)
            {
                int p1 = p + n;           //畳み込み乗算の位置
                double t = n / FS / a;    //(t-b)/a
                double alfa = pi2  * t ;  //ωt
                //ガボール減衰係数＝1 / (2 sqrt(π) σ) * exp(-t^2/σ^2)
                double gabor = exp(-t * t / (sigma * sigma)) / (2 * sqrt(pi2) * sigma);
                if (gabor >= limit && 0 <= p1 && p1 < N)
                {
                    //畳み込み乗算：実数、虚数。
                    valr += data[p1] * gabor * cos(alfa);
                    vali += data[p1] * gabor * sin(alfa);
                }
            }
            //結果の保存
            result_r[m][p] = valr / sqrt(a);
            result_i[m][p] = vali / sqrt(a);
        }
    }
}

//VC++のビューのOnPaint()
void CChildView::OnPaint() 
{
    CPaintDC dc(this);
    //0.2秒毎に、周波数を0,8,16,24,0の順で
    //1秒間の正弦波のデータ作成。
    double alfa = 0.0;
    for (int i = 0; i < N; i ++)
    {
        if (i >= 800)
            data[i] = 0.0;
        else
            data[i] = 30.0 * sin(alfa);
        alfa += (i / 200) * 8.0 * pi2 / FS;
        if (i == 0)
            dc.MoveTo(i, 40 - (int)data[i]);
        else
            dc.LineTo(i, 40 - (int)data[i]);
    }
    //ガボールウェーブレット変換
    GaborTransform();

    int K = 50;
    for (int m = 0; m < M; m ++)
    {
        int y = 100 + m * 20;
        for (int n = 0; n < N; n ++)
        {
            //瞬時値を表示
            int val = (int)result_r[m][n] / K;
            //振幅値を表示
            //int val = (int)sqrt(result_r[m][n] * result_r[m][n] + result_i[m][n] * result_i[m][n]) / K;
            if (n == 0)
                dc.MoveTo(n, y - val);
            else
                dc.LineTo(n, y - val);
        }

        CString strF;
        strF.Format(L"F:%d(HZ)", m * 2);
        dc.TextOut(1000, y - 4, strF);
    }
}