周波数空間での相関処理 #C#

相関処理とは

相関値とは「どれくらい似ているか」という数値で、ある信号と別の信号がどれくらい似ているかを調べるのが相関処理です。

この図ではleft(青色)が基準となる信号で、right(橙色)がそれをずらした信号です。この2つを相関処理するとcorr(赤色)の信号になります。
この図では複素数の実部のみを表示していますが、実際にはMagnitudeが1になるランダムな複素数が入力されています。35ポイントからのrightはleftと同じく大きさが1なので、相関値も1になっています。一方で60ポイントからのrightはleftの80%の振幅なので、相関値も0.8になっています。基準信号が正規化されていて、入力信号と基準信号が全く同じ形の場合、相関処理を行うと入力信号の大きさが求まります。

上の図の相関処理は、あらかじめ基準となる信号がわかっている場合の処理でした。これは例えばレーダーのように、自分が出した信号の反射を受信して処理するような状況に相当します(他の状況では、重力波望遠鏡のような、「シミュレーションで作った基準信号(本当の形は未知)」に対して相関処理を行う場合もあります)。

逆に、全く知らない信号に対しての相関処理も可能です。

この図ではランダムな信号を110ポイント作り、leftには後ろの100ポイントを、rightには前の100ポイントを与えています。このような信号に対して相関を行っても、正しく10ポイントの位置にピークが出ています(位置がずれているので、相関値は1よりも小さくなります)。この図ではleftに対してrightが遅れている状態ですが、leftに対してrightが進んでいる状態の場合は後ろ側にピークが出現します(一般的なFFTでの負の周波数のような形)。
今回はleftとrightは位置がずれた同じ雑音を使用していますが、一つの雑音に対してさらにleftとrightで個別の雑音が乗った信号でも同じように相関処理を行うことができます。
このような手法は例えばVLBIと呼ばれる分野で使用されています。ブラックホールを撮影したことで話題のアレです。数億光年というとてつもない距離からやってきた電波を離れた2箇所(数千kmオーダー)で受信し、相関処理を行うことで、その2箇所間の正確な距離(ミリメートルオーダー)や運動を求めることができます(これによって大陸の運動や地球の自転を極めて正確に観測できるようになりました)。

時間空間での相関処理

一番シンプルなやり方です。力技とも言います。計算力が大量に必要です。

static Complex[] CorrelationTime(Complex[] left, Complex[] right)
{
    if (right.Length < left.Length) { throw new ArgumentException(); }

    var result = new Complex[right.Length];

    for (int i = 0; i < right.Length; i++)
    {
        var tmp = Complex.Zero;
        for (int j = 0; j < left.Length; j++)
        {
            tmp += left[j].Conjugate() * right[(i + j) % right.Length];
        }
        result[i] = tmp / left.Length;
    }

    return result;
}

このコードではすべての遅延量を計算していますが、ある程度予想ができる場合は必要な範囲だけ計算すれば大幅に計算量を削減できます。

蛇足ですが、昔の人はこれを人力でやったそうです。ゼロイチの二値データ2組を二人で同時に読み上げて、3人目が耳でXORを処理して正の字を書く、というふうにやっていたらしいです(もちろん遅延量(上のコードでは変数i)を変えて何回もやり直すわけです)。昔の人すごい。

周波数空間での相関処理

フーリエ変換で周波数空間に移動し、かけ合わせた上で逆変換を行って時間空間に戻します。

static Complex[] CorrelationFreq(Complex[] left, Complex[] right)
{
    if (right.Length < left.Length) { throw new ArgumentException(); }

    var a = new Complex[right.Length];
    Array.Copy(left, a, left.Length);
    Fourier.Forward(a, FourierOptions.NoScaling);

    for (int i = 0; i < a.Length; i++)
    {
        a[i] = a[i].Conjugate() / (left.Length * right.Length);
    }

    var b = new Complex[right.Length];
    Array.Copy(right, b, right.Length);
    Fourier.Forward(b, FourierOptions.NoScaling);

    for (int i = 0; i < a.Length; i++)
    {
        b[i] = a[i] * b[i];
    }

    Fourier.Inverse(b, FourierOptions.NoScaling);

    return b;
}

このコードでは相関するたびにleftの計算を行っています。leftが毎回異なる場合はこうする他ありません。
leftがいつも同じである場合(例えばレーダー等)は、あらかじめ周波数空間で共役を取ったleftを用意しておけば、計算量を2/3まで減らすことができます。

比較

大雑把に時間空間と周波数空間での計算時間の比較を行ってみました。使用する環境やライブラリによってもかなり変わるはずなので、あくまでも参考程度に。FFTライブラリはマルチスレッドで処理を行うので公平な比較でもないですし。

時間空間での計算は単純なforループなのでポイント数と計算時間が綺麗に比例します。
周波数空間での計算も、やはり傾向としてはポイント数が増えれば計算時間も増えますが、時間空間よりは穏やかな傾きです。また、高速フーリエ変換のアルゴリズムとして、ポイント数が2のべき乗のときは計算が早く、それが下の線です。512ポイントまでは任意のポイント数の場合と比べて1桁(10倍)程度の差があります。1024ポイント以上は差が縮まりますが、それでも任意ポイントの数倍、時間空間と比べれば数百倍程度早く処理できます(512ポイント以下と1024ポイント以上での差は使用しているライブラリの特性に関する挙動だと思います)。