FFTでウェーブレット変換する

はじめに

フーリエ変換のように周波数の特性を得ることができる連続ウェーブレット変換を試しました。

連続ウェーブレット変換

母関数によっていろいろな種類がありますが、一番よく目にするのはGabor のウェーブレットです。こちらのサイトによると、（ガボールウェーブレット変換(VC++)）

を評価するみたいです。場所とスケールを変えた関数との類似度を計算していて、場所に応じた周波数の特性を求めることができます。フーリエ変換ですと、この場所が周波数高くてここは低いなどの議論はできませんが、ウェーブレット変換ならオッケーです。

実装

前の記事（EigenでFFTする - wildpieの日記）でFFTを使えるようになったので、コンボリューションは周波数領域で行います。周波数領域だとかけ算になるやつです。
たとえばガウシアンフィルタだとこんな感じです（たぶん）。

void CalcGaussKernel(VectorXd& kernel)
{
    const double sigma = 1.0;
    for (int t = 0; t < kernel.size(); t++)
    {
        kernel[t] = exp(-t*t/(2.0 * sigma))  / (sqrt(2.0*M_PI) * sigma);
    }
}

void GaussianFilter(const VectorXd& data, VectorXd& result)
{
    VectorXd kernel(data.size());
    CalcGaussKernel(kernel);

    VectorXcd data_freq = fft.fwd(data);
    VectorXcd kernel_freq = fft.fwd(kernel);
    VectorXcd product = data_freq.array() * kernel_freq.array();
    result = fft.inv(product);
}

これのKernelをGaborに変えます。

void CalcGaborKernel(VectorXcd& kernel, double a)
{
    const double sigma = 1.0;
    for (int i = 0; i < kernel.size(); i++)
    {
        double t = i / a / 8000.0;
        double omega = 2.0 * M_PI;
        kernel[i] = exp(-t*t/(sigma * sigma)) / (2.0 * sqrt(M_PI) * sigma)
                          * exp(std::complex<double>(0.0, omega * t)); // 共役?
    }
}

VectorXcd data_freq = fft.fwd(data);
MatrixXcd map(16, data.size());
for (int m = 1; m <= map.rows(); m++)
{
    double f = 2.0 * m; // 周波数, 0はどうしよう
    double a = 1.0 / f;
    VectorXcd kernel(data.size());
    CalcGaborKernel(kernel, a);

    VectorXcd kernel_freq = fft.fwd(kernel);

    VectorXcd product = data_freq.array() * kernel_freq.array();
    VectorXcd result = fft.inv(product);
    map.row(m - 1) = result / sqrt(a);
}

mapを表示するGUIを作りました。マイク音声をリアルタイムに処理しています
（c#でマイク音声をFFTする - wildpieの日記）。
合っているかあやしいですが、周波数の高い場所を検出できています。