ドキュメンテーション

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

フーリエ変換

フーリエ変換の定義

フーリエ変換は、さまざまな大きさ、周波数、位相の複素指数の和として、イメージを表現するものです。フーリエ変換は、強調、解析、復元、圧縮を含む幅広い画像処理アプリケーションで重要な役割を果たします。

f(m,n) が 2 つの離散変数 m および n の関数である場合、f(m,n) の "2 次元フーリエ変換"

F(ω1,ω2)=m=n=f(m,n)ejω1mejω2n.

の関係によって定義されます。

変数 ω1 と ω2 は周波数変数であり、単位はラジアン/サンプルになります。F(ω12) はしばしば、f(m,n) の "周波数領域" 表現と呼ばれます。F(ω12) は複素数値関数で、ω1 と ω2 の両方で周期的であり、周期は 2π となります。周期性のため、通常は πω1,ω2π の範囲のみが表示されます。F(0,0) は f(m,n) のすべての値の和です。このため、F(0,0) はしばしばフーリエ変換の "定数成分" または"DC 成分" と呼ばれます (DC は、直流 (direct current) の略です。電源の電圧が正弦関数的に変化するのとは別に、一定電圧源と考える電気工学での用語です)。

逆変換とは、変換されたイメージから元のイメージを作成する操作です。逆 2 次元フーリエ変換は、

f(m,n)=14π2ω1=ππω2=ππF(ω1,ω2)ejω1mejω2ndω1dω2.

によって求められます。

大まかに言うと、この式は f(m,n) で、異なる周波数の複素指数関数 (正弦波) の無限数の和を表現できることを示します。周波数 (ω12) の大きさと位相は、F(ω12) によって求められます。

フーリエ変換の可視化

例として、四角形の領域で 1、その他の領域で 0 になる関数 f(m,n) を考えます。ブロック線図を簡略化するため、変数 mn が離散していても、f(m,n) を連続関数として表示しています。

四角形関数

次の図は、メッシュ プロットとして、前の図に示されている四角形関数のフーリエ変換、

|F(ω1,ω2)|,

の大きさを示します。大きさのメッシュ プロットは、フーリエ変換の可視化に一般的に使用されます。

四角形関数の大きさのイメージ

プロット中心のピークは F(0,0) で、これは f(m,n) のすべての値の和になります。また、プロットには、F(ω12) 高垂直周波数よりも高水平周波数の方がエネルギーが多いことが示されています。これは、f(m,n) の水平断面が幅の狭いパルスであり、垂直断面が幅の広いパルスである事実を反映しています。幅の狭いパルスは、幅の広いパルスよりも高周波成分をもっています。

フーリエ変換を可視化する他の方法は、次に示すようにイメージを

log|F(ω1,ω2)|

として表示することです。

四角形関数のフーリエ変換の対数

対数を使用すると、F(ω12) が非常に 0 に近い領域のフーリエ変換の詳細が明確になります。

他の単純な形状のフーリエ変換の例を以下に示します。

いくつかの単純な形状のフーリエ変換

離散フーリエ変換

コンピューター上でフーリエ変換を行う場合、通常、離散フーリエ変換 (DFT) という変換を行います。離散変換とは、入力値と出力値が離散サンプルで、コンピューター操作がしやすい変換のことです。この形式の変換を使用する主な理由は 2 つです。

  • DFT の入力と出力は両方が離散で、コンピューターで操作しやすいものである。

  • DFT を計算するためには、高速フーリエ変換 (FFT) として知られるアルゴリズムがある。

DFT は、通常、有限領域 0mM1 および 0nN1 についてのみ非ゼロである離散関数 f(m,n) に対して定義されます。2 次元の MN 列の DFT と MN 列の逆 DFT の関係は、以下の式で表せます。

F(p,q)=m=0M1n=0N1f(m,n)ej2πpm/Mej2πqn/N   p=0, 1, ..., M1q=0, 1, ..., N1

および

f(m,n)=1MNp=0M1q=0N1F(p,q)ej2πpm/Mej2πqn/N   m=0, 1, ..., M1 n=0, 1, ..., N1

F(p,q) の値は、f(m,n) の DFT 係数です。ゼロ周波数係数の F(0,0) は、"DC 成分" と呼ばれます。DC は、Direct Current (直流) を略した電気工学用語です (MATLAB® の行列インデックスは、0 からでなく 1 から始まります。このため、行列の要素 f(1,1) および F(1,1) は、数学上の量としては f(0,0) および F(0,0) に対応します)。

MATLAB 関数 fftfft2fftn は、1 次元 DFT、2 次元 DFT、N 次元 DFT をそれぞれ計算する高速フーリエ変換アルゴリズムを実装します。関数 ifftifft2ifftn は逆 DFT を計算します。

フーリエ変換との関係

DFT 係数 F(p,q) は、フーリエ変換 F(ω12) のサンプルです。

F(p,q)=F(ω1,ω2)|ω1=2πp/Mω2=2πq/Np=0,1,...,M1q=0,1,...,N1

離散フーリエ変換の可視化

  1. フーリエ変換の定義」の例の関数 f(m,n) と類似した行列 f を作成します。f(m,n) は、四角形の内部が 1、外部が 0 であることを思い出してください。そのため、バイナリ イメージを使用して、f(m,n) を表示します。

    f = zeros(30,30);
    f(5:24,13:17) = 1;
    imshow(f,'InitialMagnification','fit')

  2. 次のコマンドを使用して、f の 30 行 30 列の DFT を計算し、可視化しましょう。

    F = fft2(f);
    F2 = log(abs(F));
    imshow(F2,[-1 5],'InitialMagnification','fit');
    colormap(jet); colorbar

    パディングなしの離散フーリエ変換

    このプロットは「フーリエ変換の可視化」に表示されているフーリエ変換とは異なります。まず、フーリエ変換のサンプルが、非常に粗くなっています。2 番目に、ゼロ周波数係数が、従来設定されている図の中心の位置でなく、左上隅に位置しています。

  3. フーリエ変換の細かいサンプリングを得るには、DFT を計算するときに f をゼロ パディングします。ゼロ パディングと DFT の計算は、次のコマンドを使用すれば、1 つの手順で実行できます。

    F = fft2(f,256,256);

    このコマンドは、DFT を計算する前に f をゼロ パディングして、サイズを 256 行 256 列にします。

    imshow(log(abs(F)),[-1 5]); colormap(jet); colorbar

    パディング付き離散フーリエ変換

  4. しかし、ゼロ周波数係数は、中央ではなく左上隅に表現されたままです。関数 fftshift を使用すると、この問題を回避できます。この関数は、F の 4 象限を入れ換えて、ゼロ周波数係数を中央にします。

    F = fft2(f,256,256);F2 = fftshift(F);
    imshow(log(abs(F2)),[-1 5]); colormap(jet); colorbar

    結果としてできるプロットは、「フーリエ変換の可視化」のプロットのようになります。

フーリエ変換の応用

この節は、フーリエ変換の画像処理に関連したいくつかの応用例を紹介します。

線形フィルターの周波数応答

線形フィルターのインパルス応答のフーリエ変換は、フィルターの周波数応答を与えます。関数 freqz2 は、フィルターの周波数応答を計算し、表示します。ガウスたたみ込みカーネルの周波数応答は、このフィルターが低域を通過させ、高域を減衰させることを示しています。

h = fspecial('gaussian');
freqz2(h)

ガウス フィルターの周波数応答

線形フィルター処理、フィルター設計および周波数応答についての詳細は、「周波数領域での線形フィルターの設計」を参照してください。

高速たたみ込み

フーリエ変換の重要な性質は、2 つのフーリエ変換の乗算が関連した空間関数のたたみ込みに対応していることです。この性質は、高速フーリエ変換と共に高速たたみ込みアルゴリズムに対する基礎となっています。

    メモ:   FFT ベースのたたみ込み法は、入力量が大きい場合によく使用されます。少ない入力の場合、imfilter を使用することで、一般に高速化できます。

このことを示すため、次の例で、AB のたたみ込みを行います。ここで A は M 行 N 列で、B は P 行 Q 列の行列です。

  1. 2 つの行列を作成します。

    A = magic(3);
    B = ones(3);
  2. AB が少なくとも (M+P-1) 行 (N+Q-1) 列になるようにゼロ パディングします (fft2 は 2 のべき乗のサイズで高速処理を行うので、AB をゼロ パディングして、2 のべき乗のサイズにすることがよくあります)。この例は、8 行 8 列になるように行列をパディングします。

    A(8,8) = 0;
    B(8,8) = 0;
  3. fft2 を使用して AB の 2 次元 DFT を計算し、2 つの DFT を一緒に乗算してから、ifft2 を使用して結果の逆 2 次元 DFT を計算します。

    C = ifft2(fft2(A).*fft2(B));
  4. 結果の非ゼロ部分を抽出し、丸め誤差に原因する虚数部を除去します。

    C = C(1:5,1:5);
    C = real(C)

    この例では以下の結果が作成されます。

    C =
    
        8.0000    9.0000   15.0000    7.0000    6.0000
       11.0000   17.0000   30.0000   19.0000   13.0000
       15.0000   30.0000   45.0000   30.0000   15.0000
        7.0000   21.0000   30.0000   23.0000    9.0000
        4.0000   13.0000   15.0000   11.0000    2.0000

イメージ内の特徴の検出

フーリエ変換は、たたみ込みと非常に近い関係のある相関を計算するときにも使用されます。相関は、イメージ内の特徴の位置を求めるときにも使用されます。この操作の中で、相関は、"テンプレート マッチング" と呼ばれます。

次の例は、文字列を含むイメージの中に文字 "a" の配置されている位置を相関を使って検出する方法を示しています。

  1. サンプル イメージを読み取る。

    bw = imread('text.png');
  2. イメージから文字 "a" を抜き取ることにより、マッチング法のためのテンプレートを作成します。

    a = bw(32:45,88:98);

    imcrop の対話型バージョンを使ってテンプレート イメージを作成することもできます。

    次の Figure は、元のイメージとテンプレートを示します。

    imshow(bw);
    figure, imshow(a);

    イメージ (左) と相関に使用するテンプレート (右)

  3. テンプレート イメージと元のイメージとの相関を、テンプレート イメージを 180 o 回転させ、「高速たたみ込み」で説明した FFT ベースのたたみ込み法を使用することにより計算します。

    たたみ込みは、たたみ込みカーネルを 180o 回転すると、相関と等価になります。テンプレートをイメージと一致させるため、関数 fft2 や関数 ifft2 を使用します。

    C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256)));

    次のイメージは、相関結果を示しています。イメージ内の明るい部分は、文字が存在している部分です。

    figure, imshow(C,[]) % Scale image to appropriate display range.

    相関を計算したイメージ

  4. イメージ内のテンプレートの位置を表示するには、最大ピクセル値を検出し、この値より小さいしきい値を設定します。これらのピークの位置は、しきい値の相関イメージの中で、白色スポットで定義されています (この図では位置を発見しやすいように、しきい値のイメージを膨張させて、点のサイズを拡大しています)。

    max(C(:))
    ans =
    	68.0000
    thresh = 60; % Use a threshold that's a little less than max.
    figure, imshow(C > thresh) % Display pixels over threshold.

    テンプレートの位置を示す相関を計算し、しきい値を適用したイメージ

この情報は役に立ちましたか?