11 周波数領域処理とFFT

前章までは一貫して空間領域（spatial domain）で作業を行ってきた：画像は画素の配列であり、フィルタリングは近傍における重み付き演算であり、すべての問いは「どこのグレースケール値がいくらか」に集中していた。本章では視点を切り替える——周波数領域（frequency domain）が関心を持つのは「どこ」ではなく「どのくらい速いか」である：画像のグレースケール値が空間的にどの程度速く変化するか？緩やかに変動する背景は低周波であり、鮮鋭なエッジと微細で密なテクスチャは高周波である。この視点は単なる数学的な遊びではなく、空間領域では対処できない実際的な問題を解決する。生産ラインの画像が周期的な電磁干渉に汚染され、画像全体に斜めの細かい縞模様が重畳された場合を想像してほしい。空間領域ではこの干渉は画像全体に存在し、画像の内容と絡み合っているため、平均値フィルタでも消し去ることができず、中央値フィルタも同様に無力である。しかし周波数領域では、この画像全体の干渉は2つの孤立した輝点に収縮する——それらを小さな円で切り取れば、干渉はほぼ損失なく消滅する。本章ではこの「キラーアプリケーション」を完全に実演するが、その前にまず周波数領域の言葉を学ぼう。

実験に用いるシーン（図 11.1）は512×512の合成画像である：大きな暗い矩形と大きな円が低周波のブロック状構造を提供し、中央部は周期が正確に8 pxの縦縞帯（既知の周波数の高周波成分）、下方にはドットマトリクスの文字「FFT」（テキスト状の詳細）が配置されている。各成分はスペクトル上に識別可能な特徴を残す。

図 11.1: 512×512の実験シーン：大きな矩形と円（低周波のブロック状構造）、周期8 pxの縦縞帯（既知の高周波）、ドットマトリクス文字「FFT」（詳細成分）。

11.1 2次元フーリエ変換

フーリエ変換の核心的な考え方は、任意の信号を異なる周波数の正弦波の重ね合わせに分解できることである。\(W\times H\)の離散画像\(f[n,m]\)（\(n\)を行、\(m\)を列とする）に対し、2次元離散フーリエ変換（2D Discrete Fourier Transform, DFT）は次のように定義される：

\[ F[u,v] = \sum_{n=0}^{H-1}\sum_{m=0}^{W-1} f[n,m]\, e^{-j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)}, \]

ここで\(u,v\)は水平方向と垂直方向の周波数インデックスである。逆変換（inverse DFT）は周波数成分を再び重ね合わせて画像に戻す：

\[ f[n,m] = \frac{1}{WH}\sum_{u=0}^{W-1}\sum_{v=0}^{H-1} F[u,v]\, e^{+j 2\pi \left(\frac{um}{W} + \frac{vn}{H}\right)}. \]

定義に従って直接DFTを計算するには\(O(N^2)\)回の演算（\(N=WH\)）が必要であり、512×512の画像では約\(7\times 10^{10}\)回に達する——これは許容できない。高速フーリエ変換（Fast Fourier Transform, FFT）は回転因子の対称性を利用して計算量を\(O(N\log N)\)に削減する。これこそが周波数領域手法を実用可能にした前提であり、本章の実験で画像サイズを2の整数べきにした理由でもある。

\(F[u,v]\)は複素数であり、工学的にはほぼ常に2つのより直感的な量に分解して扱う。振幅スペクトル（magnitude spectrum）\(|F[u,v]| = \sqrt{\mathrm{Re}^2 + \mathrm{Im}^2}\)は「周波数\((u,v)\)の正弦波成分がどのくらい存在するか」を示し；位相スペクトル（phase spectrum）\(\arg F[u,v]\)は「その成分がどこに位置するか」を示す——同じ組の正弦波でも位相が異なれば、重ね合わされた画像は完全に異なるものになる。フィルタリング操作は主に振幅に対して行われるが、画像を再構成する際には位相をそのまま保存しなければならない。さもなければ画像の構造は判別不能になる。

1つの周波数成分を特に取り上げる価値がある：\(u=v=0\)のとき指数項は常に1であるため、\(F[0,0] = \sum f[n,m]\)となり、これはちょうど画像全体のグレースケール値の和、つまり\(WH\)倍の平均グレースケール値である。これが直流成分（DC component）である——チャプター 6 で「カーネルの重みの和は直流ゲインに等しい」と述べたのはまさにここに由来する：直流ゲインが1のフィルタは\(F[0,0]\)を変化させないため、画像の平均輝度も変化させない。

\(|F|\)を直接画像として描画してもほとんど何も見えない：直流成分のエネルギーはほぼすべての高周波成分よりも数桁大きいため、線形のグレースケールマッピングではスペクトル全体が中心の1点を除いてほぼ真っ黒になる。標準的な手法は対数スペクトル\(\log(1+|F|)\)を表示することであり、これは巨大なダイナミックレンジを肉眼で識別可能な区間に圧縮する——本章のすべてのスペクトル図はこの方法で描画されている。

11.2 スペクトルの読み方

図 11.2 は実験シーンの対数振幅スペクトルである（直流成分は画像中心に移動されている）。スペクトルの読み方はすぐに習得できるスキルである——シーンと照らし合わせて各特徴を1つずつ識別しよう。

図 11.2: 実験シーンの対数振幅スペクトル（DC中心）。中心の輝点は低周波エネルギーであり、水平/垂直の輝いた十字は矩形と縞帯の水平・垂直なエッジに由来し、水平軸上の±64にある対称な点対は8 pxの縦縞の特徴である。

中心の輝点：低周波エネルギー。大きな矩形や円のような緩やかに変化するブロック状構造は、エネルギーのほぼすべてを中心付近に集中させる。
水平方向と垂直方向の輝いた十字：矩形の水平/垂直エッジと縞帯の水平境界線はすべて「ある方向に沿った急激な変化」であり、周波数領域では急激な変化はエッジに垂直な方向に線状に広がる。
水平軸上の対称な輝点対：純粋な周期的な縞模様は周波数領域では一対のパルスになる——これが本章で最も重要な対応関係である。

この輝点対の位置は予測してから検証することができる。周期が\(T\)画素の縞模様は、幅\(W\)の画像内でちょうど\(W/T\)回繰り返されるため、そのスペクトルピークは水平周波数インデックス

\[ u = \frac{W}{T} = \frac{512}{8} = 64 \]

の位置、つまり中心の両側\(\pm 64\)の位置に現れる。実験では中心の低周波領域を除外してスペクトル全体で最大振幅を探索したところ、実測されたピーク位置は\((\Delta u, \Delta v) = (-64, 0)\)であり——予測と完全に一致した。この「既知の周期→スペクトルピーク位置」の換算は、生産ラインの周期的な干渉をトラブルシューティングする際に直接使用できる武器である：干渉縞の画素周期を測定すれば、スペクトル上のどこを探せばよいかがわかる。

共役対称性（conjugate symmetry）：実数値画像のスペクトルは\(F[-u,-v] = F^*[u,v]\)を満たすため、振幅スペクトルは中心に関して対称になる——これがスペクトルピークが常に対で現れる理由である。後でノッチフィルタリングを行う際には、各干渉ピークの共役ピークも同時に処理しなければならない。さもなければ逆変換の結果は実数画像ではなくなる。

11.3 周波数領域フィルタリング

周波数領域フィルタリングの理論的な基礎は畳み込み定理（convolution theorem）である：空間領域における畳み込みは周波数領域における点ごとの乗算に等しく、

\[ f * h \;\longleftrightarrow\; F \cdot H. \]

である。これはチャプター 6 のすべての線形フィルタリングを周波数領域で実行できることを意味する：画像にFFTを施し、スペクトルをフィルタの周波数応答\(H\)と点ごとに乗算し、次に逆変換を行う。この2つの経路は数学的に完全に等価であり、工学的なトレードオフは計算量にある。\(K\times K\)のカーネルを用いた空間領域の畳み込みは1画素あたり\(O(K^2)\)回の乗加算を必要とするのに対し、周波数領域経路の計算コストはカーネルサイズに依存しない——小さなカーネルには空間領域を、大きなカーネル（数十画素以上）には周波数領域を使用する。

最も単純な周波数領域フィルタは理想ローパスフィルタ（ideal lowpass filter）である：直流を中心とした半径\(D_0\)以内では\(H=1\)、それ以外では\(H=0\)であり、高周波を一括してゼロにする；理想ハイパスフィルタ（ideal highpass）はその正反対である。遮断半径30 px（正規化周波数\(30/256 \approx 0.1172\)）で実験を行った結果を図 11.3 に示す。

ローパスの結果は予測を確認する：縞の周波数64は遮断半径30のはるか外側にあるため、縞帯全体が均一な灰色に平滑化される。しかし矩形、円、文字の周りの同心円状の波紋に注意してほしい——これはギブスリンギング（Gibbs ringing）であり、理想フィルタが払わなければならない代償である。原因は明らかである：周波数領域における矩形の急激な遮断に対応する空間領域の関数はsinc関数であり、sinc関数は無限に続く振動するサイドローブを持つ；畳み込み定理により、周波数領域での乗算は空間領域でこのテールを持つカーネルとの畳み込みに等しいため、すべてのエッジに波紋の列が引きずり出される。このため理想フィルタは工学的にほとんど使用されない——代わりにガウシアンまたはバターワース（Butterworth）フィルタが使用される。これらの周波数応答は滑らかに遷移し、空間領域のカーネルは振動しなくなり、代償は遮断エッジがそれほど「鮮鋭」ではなくなることだけである。

リンギングの教訓は一文でまとめられる：周波数領域で急激に遮断するほど、空間領域での引きずりは長くなる。これはチャプター 6 で述べた「ガウシアンカーネルの周波数応答はサイドローブを持たないため、平滑化の効果はクリーンである」という観察の同じコインの裏側である。

ハイパスの結果も同様にスペクトルの言葉で説明できる：半径30以内の低周波をゼロにすると、ブロック状構造はエッジの輪郭（エッジは高周波である）に還元されるが、縞帯はほぼそのまま残る——その周波数64は遮断周波数30より高く、ハイパスの通過帯域内にある。このためハイパスフィルタリングはエッジ強調と背景抑制の手段としてよく使用される。

11.4 周期ノイズのノッチフィルタリング

本章の中心的な実演に移ろう。シーンに対角方向の正弦波干渉\(40\sin\!\big(2\pi(r+c)/8\big)\)（対角方向の周期は約5.7 px）を重畳し、電磁干渉や機械振動による周期的な縞模様の汚染をシミュレートする。ノイズのある画像と元の画像のRMSEは28.13であり——理論値\(40/\sqrt{2} \approx 28.3\)（正弦波のRMS振幅）と一致している。空間領域フィルタリングはこれに対して無力である：干渉の空間周波数はシーン内の有用な8 pxの縞の周波数と非常に近いため、干渉を抑圧するのに十分な強さの平滑化は有用な縞も同時に消し去ってしまう。

しかし周波数領域では、この干渉はたった2つの点である。図 11.4 に完全な手順を示す。

まず予測する：干渉は行方向と列方向の両方で周期8 pxであるため、スペクトルピークは\((\pm 512/8, \pm 512/8) = (\pm 64, \pm 64)\)に現れるはずである。スペクトルを探索すると\((-64, -64)\)とその共役ピーク\((+64, +64)\)が得られ——再び完全に一致した。次に4段階のノッチフィルタリング（notch filtering）手順を行う：

ノイズのある画像にFFTを施す（DC中心）；
中心の低周波領域を除外し、最大振幅の干渉ピークを探索する；
ピーク位置を中心とした半径5 pxの小さな円内の複素スペクトルをゼロにし、共役ピークに対しても同じ処理を行う；
逆変換で空間領域に戻し、[0, 255]にクランプする。

結果（図 11.4 (c)）は元の画像とほとんど区別がつかない：RMSEは28.13から0.61に低下し、干渉エネルギーの98%以上が除去され、有用な縦縞（ピークは水平軸上の \(\pm 64\)、干渉ピークからはるかに離れている）は無傷で残っている。これこそが周波数領域手法を代替不能にする理由である：2つの成分が空間領域での「速さ」の点で非常に近くても、方向または周期がわずかでも異なれば、それらはスペクトル上で分離された点になる——空間領域フィルタリングは「速さ」だけで一括して遮断するのに対し、周波数領域のノッチは狙った場所を正確に打ち抜くことができる。残留する0.61の誤差は、切り取られた円内で同時に失われた少量の画像本来のエネルギーに由来する——これはノッチの半径は小さめにするほうが安全であることを示唆している。

11.5 標本化定理の再考

周波数領域の言葉を手に入れたことで、チャプター 1 で述べた「各周期を少なくとも2回標本化する」というエイリアシング（aliasing）の直感をついに正式な記述に昇格させることができる。ナイキスト標本化定理（Nyquist sampling theorem）：信号の最高周波数成分が標本化周波数の半分を下回っていれば、標本値から元の信号を完全に再構成できる。周波数領域ではこれが最も明確になる——標本化により元の信号のスペクトルは標本化周波数間隔で周期的に複製され、信号の帯域幅が標本化周波数の半分を超えると、隣接するスペクトルの複製が移動して重なり合う：高周波成分が低周波に偽装して信号に混入し、事後的に分離することはできない。これはまたチャプター 10 で述べた「画像を縮小する前に平滑化しなければならない」という規則も説明する：間引きによる縮小は標本化周波数を下げることに等しく、プレフィルタリングはまずスペクトルを新しいナイキスト限界以内に刈り込むことを意味する——詳細が失われても、詳細が偽の低周波アーチファクトに変わるよりはましである。

11.6 SciVisionによる実装

本章のすべての実験はSCIMV::SciSvFFTクラスによって実行される；核心的な呼び出しは以下の通りである：

SCIMV::SciSvFFT fft;
SciImage fftImg;
// mode=0：直流成分をスペクトル中心に移動；fftImgは実部/虚部がインターリーブされた2チャンネル32F画像
fft.ApplyFFTGeneric(src, 0, &fftImg, NULL);

// 理想ローパス/ハイパス：frequencyは（半幅W/2に対する）正規化遮断周波数、0.1172 → 実測半径30.0 px
float cutoff = 30.0f / (512 / 2);
SciImage lp;
fft.GenerateLowpass(cutoff, 0, 512, 512, &lp);

// 畳み込み定理の直接的な体現：スペクトルをフィルタと点ごとに乗算
SciImage fftLP, out32;
fft.ApplyFFTConvolution(fftImg, lp, &fftLP);
// 逆変換：必ずSCI_IMAGE_32Fを要求し、自分で[0,255]にクランプすること
fft.ApplyIFFTGeneric(fftLP, 0, &out32, true, SCI_IMAGE_32F);

ApplyFFTGenericが出力する複素スペクトルはImageData()とStep()を通じて直接読み書きできる——\(r\)行\(c\)列の実部と虚部はそれぞれ浮動小数点オフセットr*stride + 2*cとr*stride + 2*c + 1に位置する。セクション 11.4 のノッチフィルタリングはまさにこの方法で、2つの小さな円内の複素数を手動でゼロにすることで実装されている——SDKに既製のノッチインターフェースはないが、必要もない。

このAPIを使用する際に実験的に確認された3つの落とし穴があるので、そのまま記録する：

IFFTの8U出力はmin-max正規化される：直接8ビット出力を要求するとグレースケール値が全体的に伸張され、往復（FFT→IFFT）のRMSEは驚くべき33.62になる——アルゴリズムが壊れているように見えるが、実際には正規化によるものである。SCI_IMAGE_32Fを要求して生の浮動小数点結果を取得し、自分でクランプすると、往復のRMSEは0.0000になる——変換自体は損失がない。
ApplyFFTGenericに内蔵された振幅画像はほぼ真っ黒である：そのimageMagnitude出力は線形正規化を使用しており、直流成分が支配的で他の周波数はすべて見えない。対数スペクトル\(\log(1+|F|)\)は複素数データから自分で計算する必要がある。
RemovePeriodicPatternsByFFTは使用できない：このインターフェースはヘッダファイルで明示的に「未実装」と記載されている；周期ノイズの除去には本章の手順に従って手動でノッチフィルタリングを行うこと。

本章のすべての画像を生成する完全なプロジェクトはcode/fourier/に配置されている；スペクトルピークの予測値と実測値の対照はコンソールに直接出力される。

産業事例：織物表面の織り目干渉

ある織物表面欠陥検査プロジェクトでは、織物自体の経緯の織り目が強力な周期構造であり、そのグレースケール変動は糸切れや油汚れなどの欠陥信号をはるかに超えていたため、元の画像上では欠陥はほとんど見えなかった。空間領域フィルタリングはジレンマに陥った：小さなカーネルでは織り目をきれいに除去できず、大きなカーネルでは欠陥が織り目と一緒にぼやけてしまった。周波数領域手法に切り替えることで問題は解消した——織り目の基本周波数とその高調波はスペクトル上に位置の安定した輝点対の集合を形成し；各点対に小さな半径のノッチを適用して空間領域に逆変換すると、織り目は全体的に除去され、残差画像上で欠陥が明瞭に見えるようになった。現場調整で得られた重要な教訓はノッチの半径は小さめにするほうが安全であることである：欠陥エネルギーも織り目のピーク付近に分布しており、半径を大きくすると織り目はより完全に除去できるが、欠陥エネルギーを飲み込んで検出のコントラストを低下させてしまう——ノッチのサイズはスペクトルピークをちょうど覆う程度にするのがよい。

11.7 まとめ

周波数領域は「どのくらい速いか」に答え、空間領域は「どこ」に答える：振幅スペクトルは各周波数成分の存在量を示し、位相スペクトルはそれらの位置を示す；直流成分\(F[0,0]\)は全グレースケール値の和であり、対数表示\(\log(1+|F|)\)はスペクトルを見るための標準的な方法である。
スペクトルは読むことができ、さらに予測することもできる：幅\(W\)の画像内の周期\(T\)の縞は\(u=W/T\)に共役なピーク対を生じる——8 pxの縞は±64と予測され、実測は完全に一致した。
畳み込み定理\(f*h \leftrightarrow F\cdot H\)は2つの領域をつなぐ：小さなカーネルには空間領域の畳み込みを、大きなカーネルには周波数領域経路を使用する；周波数領域での急激な遮断はギブスリンギング（sincのテール）を引き起こすため、工学的には滑らかに遷移するガウシアン/バターワースフィルタで理想フィルタを代替する。
ノッチフィルタリングは周期干渉に対する特効薬である：画像全体の周期的な汚染はスペクトル上では数個の点に過ぎず、小さな円をゼロにして逆変換すればよい——実験ではRMSEは28.13から0.61に低下し、これはどの空間領域フィルタも達成できない。
周波数領域で記述されるナイキスト定理：信号の最高周波数は標本化周波数の半分を下回らなければならず、さもなければ移動したスペクトルの複製が重なり合ってエイリアシングが生じる——これが間引きによる縮小の前にプレフィルタリングが必須である根本的な理由である。

FFTを実用可能にした高速アルゴリズムはCooleyとTukeyの1965年の古典的な論文(Cooley と Tukey 1965)に由来し、DFTの計算量を\(O(N^2)\)から\(O(N\log N)\)に削減した；デジタル画像処理におけるフーリエ変換の体系的な解説（周波数領域フィルタリングとノッチフィルタを含む）はGonzalezとWoodsの教科書(Gonzalez と Woods 2018)にある。産業画像処理における周波数領域手法の体系的な解説（最適フィルタと周波数領域特徴を含む）については、Stegerらの著作(Steger, Ulrich, と Wiedemann 2018)をさらに参照されたい。