Main Content

このページの内容は最新ではありません。最新版の英語を参照するには、ここをクリックします。

フーリエ変換

フーリエ変換の定義

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

f(m,n) が 2 つの離散空間変数 m および n の関数である場合、f(m,n) の "2 次元フーリエ変換" は次の関係によって定義されます。

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

変数 ω1 および ω2 は周波数変数で、単位はラジアン/サンプルです。F(ω1,ω2)f(m,n) の "周波数領域" 表現とも呼ばれます。F(ω1,ω2) は複素数値関数で、ω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) で、異なる周波数の複素指数関数 (正弦波) の無限数の和を表現できることを示します。周波数 (ω1,ω2) での寄与の大きさと位相は F(ω1,ω2) によって求められます。

フーリエ変換の可視化

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

四角形関数

Plot of rectangular function f(m,n) in the spatial domain

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

|F(ω1,ω2)|,

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

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

Mesh plot of the magnitude of the Fourier transform of the rectangular function f(m,n), plotted as a function of the horizontal frequency, ω1, and vertical frequency, ω2

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

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

log|F(ω1,ω2)|

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

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

2-D pot of the log of the Fourier transform of the rectangular function f(m,n), plotted as a function of the horizontal frequency, ω1, and vertical frequency, ω2

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

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

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

2-D plots of the Fourier transforms of a rectangle tilted forty-five degrees, a circle centered at (0, 0), and an X shape centered at (0, 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(ω1,ω2) のサンプルです。

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")

    Binary image representation of f(m,n)

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

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

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

    2-D plot of the 30-by-30 discrete Fourier transform of the binary rectangular function

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

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

    F = fft2(f,256,256);

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

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

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

    2-D plot of the 30-by-30 discrete Fourier transform of the binary rectangular function with zero padding. The transform with zero padding has a finer frequency resolution.

  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)

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

Mesh plot of the magnitude of the frequency response of a Gaussian filter.

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

フーリエ変換を使用した高速畳み込みの実行

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

メモ: FFT ベースの畳み込み法は、入力量が大きい場合によく使用されます。通常、少ない入力の場合、関数 imfilter を使用する方が高速です。

2 つのシンプルな行列 AB を作成します。A は M 行 N 列で、B は P 行 Q 列の行列です。

A = magic(3);
B = ones(3);

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

A(8,8) = 0;
B(8,8) = 0;

関数 fft2 を使用して AB の 2 次元 DFT を計算します。2 つの DFT を乗算し、この結果の逆 2 次元 DFT を関数 ifft2 を使用して計算します。

C = ifft2(fft2(A).*fft2(B));

結果の非ゼロ部分を抽出し、丸め誤差に原因する虚数部を除去します。

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

    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

FFT ベースの相関を計算してイメージの特徴の位置を検出

この例では、フーリエ変換を使用して、畳み込みと密接に関連する相関を計算する方法を説明します。相関を使用すると、イメージ内の特徴の位置を検出できます。上記に関連して、相関は "テンプレート マッチング" とよく呼ばれます。

サンプル イメージをワークスペースに読み取ります。

bw = imread('text.png');

イメージから文字 "a" を抜き取ることにより、マッチング法のためのテンプレートを作成します。imcrop の対話形式の構文を使用してテンプレートを作成することもできます。

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

テンプレート イメージと元のイメージとの相関を、テンプレート イメージを 180 度回転させ、FFT ベースの畳み込み法を使用することにより計算します。畳み込みは、畳み込みカーネルを 180 度回転すると、相関と等価になります。テンプレートをイメージと一致させるため、関数 fft2 や関数 ifft2 を使用します。結果のイメージ内の明るい部分は、文字が存在している部分です。

C = real(ifft2(fft2(bw) .* fft2(rot90(a,2),256,256)));
figure
imshow(C,[]) % Scale image to appropriate display range.

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

max(C(:))
ans = 68
thresh = 60; % Use a threshold that's a little less than max.
D = C > thresh;
se = strel('disk',5);
E = imdilate(D,se);
figure
imshow(E) % Display pixels with values over the threshold.

参考

| |

関連するトピック