---

【初心者向け】Python NumPy/SciPyを使ったFFT（高速フーリエ変換）入門

はじめに

データ分析や信号処理の世界に足を踏み入れると、「FFT」という言葉をよく耳にするでしょう。FFTは「高速フーリエ変換」の略で、時間とともに変化する信号を、どのような周波数成分が含まれているかに分解するための強力なツールです。音声、画像、振動、金融市場のデータなど、様々な分野で応用されています。

しかし、「フーリエ変換」や「FFT」と聞くと、難解な数式を想像して尻込みしてしまう方もいるかもしれません。安心してください。Pythonを使えば、複雑な数式を知らなくても、FFTを比較的簡単に実行し、その結果を分析することができます。特に、科学技術計算で広く使われているNumPyとSciPyライブラリは、FFTを扱うための強力な機能を提供しています。

この記事は、Pythonを使ったFFTが初めてという方を対象としています。FFTの基本的な考え方から、NumPyおよびSciPyを使った具体的な実装方法、そして得られた結果の読み取り方までを、豊富なコード例とともに丁寧に解説します。この記事を読み終える頃には、あなたもPythonを使って信号の周波数成分を分析できるようになっているはずです。

この記事で学ぶこと:

• フーリエ変換とは何か、なぜFFTが必要なのか
• PythonでFFTを行うための準備（NumPy, SciPy, Matplotlibのインストール）
• NumPyを使ったFFTの基本的な実行方法 (numpy.fft モジュール)
• FFT結果の解釈（周波数スペクトル、振幅、パワー）
• SciPyを使ったFFTの実行方法 (scipy.fft モジュール)
• データ分析におけるFFTの重要な注意点（サンプリングレート、ナイキスト周波数、データ長）
• 実践的なFFTの使い方（窓関数など）

さあ、Pythonを使って信号の「声」を聞く旅に出かけましょう！

フーリエ変換とは何か？

まずは、FFTの元となる「フーリエ変換」の基本的な考え方から始めましょう。

私たちが普段扱っているデータや信号の多くは、「時間とともに変化する値」として記録されています。これを「時間領域のデータ」と呼びます。例えば、マイクで録音した音声データは、時間の経過とともに空気の振動（音圧）がどのように変化したかを記録したものです。心電図の波形も、時間の経過とともに心臓の電気信号がどのように変化したかを示しています。

フーリエ変換は、この時間領域のデータを「周波数領域のデータ」に変換する数学的な手法です。周波数領域のデータとは、「元の信号にどのような周波数の成分が、それぞれどれくらいの強さで含まれているか」を示したものです。

例えるなら、時間領域のデータが「オーケストラの演奏全体の音」だとすると、フーリエ変換は「その演奏を構成している個々の楽器（周波数）の音とその大きさ（強さ）」に分解する作業です。フルートの音（高周波数）、チェロの音（低周波数）、それぞれの音がどのくらいの音量で鳴っているかを分析するイメージです。

この変換の根底にあるのは、「どんなに複雑な波形も、様々な周波数を持つ単純なサイン波（sin波）とコサイン波（cos波）を重ね合わせることで表現できる」という驚くべき事実です。フーリエ変換は、元の信号を構成するサイン波とコサイン波の「リスト」（それぞれの周波数と、その波の振幅および位相）を取り出す作業だと言えます。

なぜ周波数領域への変換が重要なのでしょうか？それは、特定の周波数成分に注目することで、信号の特性をより深く理解したり、ノイズを除去したり、信号を圧縮したりといった処理が可能になるからです。例えば、音声信号から特定の周波数帯のノイズを除去したり、機械の振動データから異常を示す特定の周波数成分を検出したりするのに役立ちます。

コンピュータで信号を扱う場合、信号は連続的な波形ではなく、一定の間隔でサンプリング（標本化）された離散的な値の集まりとして表現されます。この離散的なデータに対して行われるのが「離散フーリエ変換 (DFT: Discrete Fourier Transform)」です。

DFTは、理論上はどのような離散信号に対しても適用できますが、データ点数が多い場合に計算量が膨大になるという欠点があります。データ点数をNとすると、DFTの計算量は概ねN²に比例します。Nが大きくなると、計算時間が現実的ではなくなってしまいます。

FFT（高速フーリエ変換）とは何か？

ここで登場するのが「FFT（高速フーリエ変換）」です。

FFTは、DFTを計算するための効率的なアルゴリズムです。DFTと同じ結果を得られますが、計算量が大幅に削減されます。データ点数Nが2のべき乗の場合に最も効率的で、計算量は概ね N log₂(N) に比例します。

例えば、N=1024（約1000）の場合を考えてみましょう。
* DFTの計算量: N² ≈ (10³)² = 10⁶
* FFTの計算量: N log₂(N) = 1024 * log₂(1024) = 1024 * 10 ≈ 10⁴

計算量が100分の1になることがわかります。データ点数がさらに増えれば、その差はさらに顕著になります。N=1,000,000 (10⁶) の場合、DFTは (10⁶)² = 10¹²、FFTは 10⁶ * log₂(10⁶) ≈ 10⁶ * 20 ≈ 2 * 10⁷ となり、その差は歴然です。

この計算効率のおかげで、FFTはリアルタイムの信号処理や大規模なデータ分析において不可欠なツールとなりました。私たちが普段使っているデジタル機器の多く（スマートフォン、音楽プレイヤー、画像編集ソフトなど）の内部でも、FFTが様々な処理に利用されています。

NumPyやSciPyで「FFTを実行する」と言う場合、実際にはこの高速なFFTアルゴリズムを使って離散フーリエ変換を行っている、と考えて差し支えありません。

PythonでFFTを使うための準備

PythonでFFTを使うためには、以下のライブラリが必要です。

• NumPy: 数値計算、特に配列（多次元配列）を効率的に扱うための基本ライブラリです。FFT機能も提供しています。
• SciPy: 科学技術計算のための様々な機能を提供するライブラリ集です。FFTに関しても、NumPyよりもさらに特化した機能や、関連する信号処理機能（窓関数など）を提供しています。
• Matplotlib: グラフ描画ライブラリです。FFTの結果を視覚化するために使います。

これらのライブラリがインストールされていない場合は、pipコマンドを使ってインストールしてください。

bash
pip install numpy scipy matplotlib

インストールが完了したら、PythonスクリプトやJupyter Notebookで以下のコードを実行し、エラーが出ないことを確認してください。

“`python
import numpy as np
import scipy.fft as fft
import matplotlib.pyplot as plt

print(f”NumPy version: {np.version}”)
print(f”SciPy version: {scipy.version}”)
print(f”Matplotlib version: {matplotlib.version}”)
“`

エラーが出ずに各ライブラリのバージョンが表示されれば準備完了です。

NumPyを使ったFFT

まず、NumPyを使ってFFTを実行してみましょう。NumPyのfftモジュールには、FFTを行うための基本的な関数が揃っています。

主な関数：

• numpy.fft.fft(a): 1次元配列aに対してFFTを実行します。結果は複素数の配列になります。
• numpy.fft.ifft(a): fftの逆変換（逆FFT）を実行します。周波数領域のデータから時間領域のデータに戻します。
• numpy.fft.fftfreq(n, d=1.0): FFTの結果に対応する周波数ビン（各データ点がどの周波数に対応するか）の配列を生成します。nはデータ点数、dはサンプリング間隔（1/サンプリングレート）です。
• numpy.fft.fftshift(x, axes=None): FFTの結果を、0Hzを中心とするようにシフトします。プロットする際に見やすくなります。
• numpy.fft.ifftshift(x, axes=None): fftshiftの逆を行います。

では、実際にサイン波のデータに対してFFTを実行してみましょう。

実践例1: 単一のサイン波

周波数5Hzのサイン波を生成し、そのFFTを計算します。

“`python
import numpy as np
import matplotlib.pyplot as plt

サンプリングレート (Hz)

Fs = 100

サンプリング間隔 (秒)

T = 1.0 / Fs

信号の長さ (秒)

duration = 2.0

サンプル点の数

N = int(Fs * duration)

時間軸データ

t = np.linspace(0.0, duration, N, endpoint=False)

元の信号 (周波数5Hzのサイン波)

freq1 = 5
y = np.sin(freq1 * 2.0 * np.pi * t)

元の信号をプロット

plt.figure(figsize=(10, 4))
plt.plot(t, y)
plt.title(“Original Signal (5Hz Sine Wave)”)
plt.xlabel(“Time [s]”)
plt.ylabel(“Amplitude”)
plt.grid(True)
plt.show()

——– FFTを実行 ——–

FFTを実行

yf = np.fft.fft(y)

周波数軸を生成

N: サンプル点の数

T: サンプリング間隔 (1/Fs)

xf = np.fft.fftfreq(N, T)

FFT結果をプロット

plt.figure(figsize=(10, 4))

振幅スペクトル (絶対値)

plt.plot(xf, np.abs(yf))
plt.title(“FFT Result (Amplitude Spectrum)”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Amplitude |Y(f)|”)
plt.grid(True)
plt.show()

FFT結果をfftshiftでシフトしてプロット

plt.figure(figsize=(10, 4))
plt.plot(np.fft.fftshift(xf), np.fft.fftshift(np.abs(yf)))
plt.title(“FFT Result (Amplitude Spectrum, Shifted)”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Amplitude |Y(f)|”)
plt.grid(True)
plt.show()
“`

コードの解説:

1. パラメータ設定:
   • Fs: サンプリングレートです。1秒間にいくつのデータを取得するかを示します。ここでは100Hz（1秒間に100点）とします。
   • T: サンプリング間隔です。1/Fsで計算できます。ここでは0.01秒です。
   • duration: 信号の全体の長さです。ここでは2秒とします。
   • N: サンプル点の総数です。Fs * durationで計算できます。ここでは200点です。
2. 時間軸データの生成: np.linspace(0.0, duration, N, endpoint=False)を使って、0秒からduration秒までの時間をN等分した配列を作成します。endpoint=Falseとすることで、最後の点は含めず、0, T, 2T, …, (N-1)T となります。これはFFTが周期信号を仮定しているため都合が良い場合が多いです。
3. 元の信号生成: np.sin(freq1 * 2.0 * np.pi * t)で、周波数freq1（ここでは5Hz）のサイン波を生成します。2.0 * np.piをかけるのは、np.sinがラジアンを引数にとるためです。周波数fの角周波数は 2πf となります。
4. FFTの実行: np.fft.fft(y)で、生成したサイン波データyのFFTを実行します。結果yfはyと同じ長さの配列ですが、各要素は複素数です。
5. 周波数軸の生成: np.fft.fftfreq(N, T)を使って、FFT結果yfの各要素がどの周波数に対応するかを示す配列xfを生成します。Nはデータ点数、Tはサンプリング間隔です。
6. 結果のプロット: matplotlib.pyplot.plot(xf, np.abs(yf))で、FFT結果の振幅スペクトルをプロットします。FFT結果yfは複素数なので、その絶対値np.abs(yf)を取ることで振幅（パワースペクトルの平方根に相当）が得られます。
7. fftshiftによるシフト: np.fft.fftshift(xf)とnp.fft.fftshift(np.abs(yf))を使って、周波数軸とFFT結果をシフトしてからプロットします。

プロットの解釈:

• 最初のプロットは、私たちが生成した時間領域のサイン波です。
• 2番目のプロットは、FFTを直接適用した結果です。周波数軸はfftfreqで生成されたもので、0Hzが中央ではなく左端（正の周波数）と右端（負の周波数）にあります。データ点数200の場合、-50Hzから+49.5Hzまでの周波数ビンが生成されます（詳細は後述）。
• 3番目のプロットは、fftshiftを使って結果をシフトしたものです。これで0Hzが中央に来て、負の周波数が左側、正の周波数が右側に並びます。

3番目のプロットを見ると、±5Hzのところに大きなピークがあることがわかります。これは、元の信号が5Hzのサイン波だったことを正確に捉えています。サイン波は正の周波数成分と負の周波数成分の両方を持つため、±5Hzの両方にピークが現れます。

ピークの高さは約100になっています。これはデータ点数N=200の半分です。FFT結果の振幅は、通常、元の信号の振幅にデータ点数Nをかけたものになります。サイン波は実数信号なので、周波数成分は正負で半分ずつに分かれます。したがって、特定の周波数f (f > 0) の振幅を知りたい場合は、FFT結果の周波数fに対応する値の絶対値を2/N倍するか（DC成分やナイキスト周波数成分は例外）、ピーク値を読んでから後で調整する必要があります。単純に絶対値を取るだけでは、各周波数成分の「強さの相対的な関係」はわかりますが、元の信号の物理的な振幅とはスケールが異なります。後ほど正規化についても詳しく説明します。

実践例2: 複数のサイン波の合成

次に、複数の周波数を持つ信号に対してFFTを適用してみましょう。周波数5Hzと15Hzのサイン波を合成した信号を分析します。

“`python
import numpy as np
import matplotlib.pyplot as plt

パラメータ設定 (実践例1と同じ)

Fs = 100
T = 1.0 / Fs
duration = 2.0
N = int(Fs * duration)
t = np.linspace(0.0, duration, N, endpoint=False)

元の信号 (周波数5Hzと15Hzのサイン波の合成)

freq1 = 5
freq2 = 15
amplitude1 = 1.0
amplitude2 = 0.5 # 2つ目のサイン波は振幅を小さくしてみる
y_composite = amplitude1 * np.sin(freq1 * 2.0 * np.pi * t) + \
amplitude2 * np.sin(freq2 * 2.0 * np.pi * t)

元の信号をプロット

plt.figure(figsize=(10, 4))
plt.plot(t, y_composite)
plt.title(“Original Composite Signal (5Hz + 15Hz Sine Waves)”)
plt.xlabel(“Time [s]”)
plt.ylabel(“Amplitude”)
plt.grid(True)
plt.show()

——– FFTを実行 ——–

FFTを実行

yf_composite = np.fft.fft(y_composite)

周波数軸を生成

xf_composite = np.fft.fftfreq(N, T)

FFT結果をシフト

xf_composite_shifted = np.fft.fftshift(xf_composite)
yf_composite_shifted = np.fft.fftshift(yf_composite)

FFT結果をプロット (シフト済み)

plt.figure(figsize=(10, 4))
plt.plot(xf_composite_shifted, np.abs(yf_composite_shifted))
plt.title(“FFT Result (Amplitude Spectrum, Shifted)”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Amplitude |Y(f)|”)
plt.grid(True)
plt.show()
“`

プロットの解釈:

時間領域のプロットを見ると、複雑な波形になっています。これは5Hzと15Hzの波が重なり合った結果です。

FFT結果のプロットを見ると、±5Hzと±15Hzのところにピークが現れていることがわかります。±5Hzのピークは、元の信号に含まれていた5Hz成分に対応し、その高さは約100です（元の振幅1.0 × N/2 = 1.0 × 200/2 = 100）。±15Hzのピークは、15Hz成分に対応し、その高さは約50です（元の振幅0.5 × N/2 = 0.5 × 200/2 = 50）。

このように、FFTを使うことで、複雑な信号がどのような周波数成分から構成されているかを簡単に特定できます。

fftshift について補足

np.fft.fftfreq(N, T)で生成される周波数軸は、通常、0Hzから始まり、次に正の周波数成分が並び、最後に負の周波数成分が並びます。例えば、N=200, Fs=100Hzの場合、周波数軸は以下のようになります。

[ 0.0, 0.5, 1.0, …, 49.5, -50.0, -49.5, …, -0.5 ]

これはコンピュータ内部での処理の都合による並び順ですが、グラフとして表示する場合、負の周波数を左側、正の周波数を右側に配置し、0Hzを中央にするのが一般的で見やすいです。np.fft.fftshift()関数は、この並び替えを行ってくれます。

[ -50.0, -49.5, …, -0.5, 0.0, 0.5, 1.0, …, 49.5 ]

このように、fftshiftはFFTの結果自体を変更するのではなく、結果の並び順とそれに合わせて周波数軸の並び順を変えるだけです。プロットする際に利用することが多い関数です。

SciPyを使ったFFT

SciPyのfftモジュール（または古いバージョンではfftpackモジュール）もNumPyと同様のFFT機能を提供しています。多くの基本的な関数名や使い方はNumPyと共通しています。

なぜNumPyとSciPyの両方にFFT機能があるのでしょうか？
* NumPyはより基本的な数値計算をカバーしており、FFTはその一部として含まれています。
* SciPyはより高度で専門的な科学技術計算機能を提供しており、FFTはその信号処理（scipy.signalなど）やその他のモジュールと連携して使われることが多いです。SciPyのFFT実装は、特定のプラットフォームでNumPyより高性能な場合があったり、追加機能を提供していたりすることがあります（例：窓関数など）。

NumPyのnp.fftモジュールは多くの用途で十分高速ですが、大規模な計算や特定の最適化が必要な場合は、SciPyのscipy.fftを検討する価値があります。最近のバージョンでは、SciPyのscipy.fftがNumPyのnp.fftよりも推奨される傾向にあります。

SciPyの主な関数：

• scipy.fft.fft(x): NumPyと同様にFFTを実行します。
• scipy.fft.ifft(x): NumPyと同様に逆FFTを実行します。
• scipy.fft.fftfreq(n, d=1.0): NumPyと同様に周波数軸を生成します。
• scipy.fft.fftshift(x, axes=None): NumPyと同様に結果をシフトします。

使い方はNumPyとほとんど同じです。

実践例3: SciPyを使ったFFT (NumPyとの比較)

実践例2で使用した合成波データに対して、SciPyを使ってFFTを実行してみましょう。

“`python
import numpy as np
import scipy.fft as fft # scipy.fft をインポート
import matplotlib.pyplot as plt

パラメータ設定 (実践例2と同じ)

Fs = 100
T = 1.0 / Fs
duration = 2.0
N = int(Fs * duration)
t = np.linspace(0.0, duration, N, endpoint=False)

元の信号 (周波数5Hzと15Hzのサイン波の合成)

freq1 = 5
freq2 = 15
amplitude1 = 1.0
amplitude2 = 0.5
y_composite = amplitude1 * np.sin(freq1 * 2.0 * np.pi * t) + \
amplitude2 * np.sin(freq2 * 2.0 * np.pi * t)

——– SciPyでFFTを実行 ——–

FFTを実行 (scipy.fft.fft を使用)

yf_scipy = fft.fft(y_composite)

周波数軸を生成 (scipy.fft.fftfreq を使用)

xf_scipy = fft.fftfreq(N, T)

FFT結果をシフト (scipy.fft.fftshift を使用)

xf_scipy_shifted = fft.fftshift(xf_scipy)
yf_scipy_shifted = fft.fftshift(yf_scipy)

FFT結果をプロット (シフト済み)

plt.figure(figsize=(10, 4))
plt.plot(xf_scipy_shifted, np.abs(yf_scipy_shifted))
plt.title(“FFT Result (Amplitude Spectrum, Shifted) using SciPy”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Amplitude |Y(f)|”)
plt.grid(True)
plt.show()

NumPyの結果と比較 (オプション)

実践例2のyf_composite_shifted, xf_composite_shiftedと np.abs(yf_scipy_shifted), xf_scipy_shifted はほぼ同じになるはずです。

np.allclose(np.abs(yf_composite_shifted), np.abs(yf_scipy_shifted)) # True になるはず

“`

ご覧のように、NumPyを使った場合と全く同じ結果が得られます。基本的なFFTの計算においては、NumPyとSciPyのfftモジュールは同様に利用できます。

窓関数 (Window Functions) について

SciPyのfftモジュール自体ではなく、関連するscipy.signalモジュールには、FFTを行う際に非常に役立つ「窓関数」の機能があります。

窓関数とは？ なぜ必要？

FFTは、対象となる信号が「無限に続いている周期信号」であるかのように扱います。しかし、実際のFFTの入力データは、観測された信号の一部を切り取った「有限長」のデータです。

この有限長のデータをそのままFFTにかけると、データの開始点と終了点がつながっているかのように処理されます。もし開始点と終了点の値が大きく異なると、そこで不連続なジャンプが発生しているとみなされ、本来存在しない高い周波数成分（リーケージと呼ばれる現象）がFFT結果に現れてしまいます。これは、意図しないノイズのように見えたり、本来の周波数成分のピークが広がって見えたりする原因となります。

窓関数は、このリーケージを軽減するためのテクニックです。FFTをかける前に、信号データに窓関数を「かける」ことで、データの開始点と終了点の値を滑らかにゼロに近づけます。これにより、信号が周期的な波形であるというFFTの仮定により近づき、リーケージを抑えることができます。

窓関数には様々な種類があります（例: Hanning窓、Hamming窓、Blackman窓など）。どの窓関数を使うかは、分析したい信号の性質や目的に応じて選びますが、初心者としてはまずHanning窓などから試してみるのが良いでしょう。

SciPyのscipy.signal.windowsモジュールは、これらの窓関数を簡単に生成する機能を提供しています。

• scipy.signal.windows.hann(M): データ長Mのハニング窓を生成します。
• scipy.signal.windows.hamming(M): データ長Mのハミング窓を生成します。
• 他にも多くの窓関数が利用可能です。

窓関数を使ったFFTの手順は以下のようになります。

1. 元の信号データを準備する。
2. 元の信号データと同じ長さの窓関数を生成する。
3. 元の信号データと窓関数を、要素ごとに掛け合わせる。
4. 窓関数をかけた信号データに対してFFTを実行する。

実践例として、先ほどの合成波に窓関数を適用してFFTを実行してみましょう。

“`python
import numpy as np
import scipy.fft as fft
import scipy.signal.windows as windows # windows モジュールをインポート
import matplotlib.pyplot as plt

パラメータ設定 (実践例2と同じ)

Fs = 100
T = 1.0 / Fs
duration = 2.0
N = int(Fs * duration)
t = np.linspace(0.0, duration, N, endpoint=False)

元の信号 (周波数5Hzと15Hzのサイン波の合成)

freq1 = 5
freq2 = 15
amplitude1 = 1.0
amplitude2 = 0.5
y_composite = amplitude1 * np.sin(freq1 * 2.0 * np.pi * t) + \
amplitude2 * np.sin(freq2 * 2.0 * np.pi * t)

——– 窓関数を適用してFFTを実行 ——–

ハニング窓を生成 (データ長と同じサイズ)

hann_window = windows.hann(N)

元の信号に窓関数をかける

y_windowed = y_composite * hann_window

窓関数をかけた信号をプロット (元の信号と比較)

plt.figure(figsize=(10, 4))
plt.plot(t, y_composite, label=”Original Signal”)
plt.plot(t, y_windowed, label=”Windowed Signal”)
plt.title(“Signal before and after Windowing”)
plt.xlabel(“Time [s]”)
plt.ylabel(“Amplitude”)
plt.legend()
plt.grid(True)
plt.show()

FFTを実行 (窓関数適用後のデータに対して)

yf_windowed = fft.fft(y_windowed)

周波数軸を生成

xf_windowed = fft.fftfreq(N, T)

FFT結果をシフト

xf_windowed_shifted = fft.fftshift(xf_windowed)
yf_windowed_shifted = fft.fftshift(yf_windowed)

FFT結果をプロット (シフト済み, 窓関数適用後)

plt.figure(figsize=(10, 4))
plt.plot(xf_windowed_shifted, np.abs(yf_windowed_shifted))
plt.title(“FFT Result (Amplitude Spectrum, Shifted) after Windowing”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Amplitude |Y(f)|”)
plt.grid(True)
plt.show()
“`

プロットの解釈:

• 窓関数適用前後の信号を比べたプロットを見ると、窓関数適用後の信号は両端に向かって振幅が滑らかにゼロに近づいているのがわかります。
• 窓関数適用後のFFT結果を見ると、±5Hzと±15Hzにピークがあるのは同じですが、ピークの周りの裾野（他の周波数成分への広がり）が小さくなっていることが期待されます。特に、単一の周波数成分を分析する場合など、意図しないリーケージ成分を抑えたい場合に窓関数は有効です。

窓関数を適用した場合、ピークの高さが元の信号の振幅と直接対応しなくなることがあります（窓関数の種類によってスケールが変わるため）。正確な振幅を知りたい場合は、窓関数による影響を考慮して正規化を行う必要がありますが、周波数成分の位置を特定する目的であれば、窓関数は非常に有効です。

FFT結果の解釈と注意点

FFTの結果を正しく理解し、誤った結論を導かないために、いくつか重要なポイントがあります。

振幅スペクトルとパワースペクトル、正規化

FFTの結果yfは複素数です。この複素数の各要素は、その周波数成分の振幅と位相の情報を持っています。

• 振幅スペクトル: 各周波数成分の「強さ」を示します。これは、複素数の絶対値np.abs(yf)で得られます。一般的にFFTの結果をプロットする際は、この振幅スペクトルを見ることが多いです。
• パワースペクトル: 各周波数成分の「パワー」を示します。これは、振幅の2乗np.abs(yf)**2で得られます。エネルギーやパワーに注目したい場合に利用されます。物理的な意味合いを持つことが多いです。

正規化: FFTの結果の振幅やパワーの絶対値は、入力データのスケールやデータ点数Nに依存します。元の信号の物理的な振幅やパワーを知りたい場合は、正規化が必要です。

一般的な正規化の方法はいくつかあります。

• 振幅スペクトル: np.abs(yf) / N または実数信号の場合 np.abs(yf) / (N/2) (DC成分とナイキスト周波数成分は np.abs(yf) / N)
• パワースペクトル: np.abs(yf)**2 / N**2 または実数信号の場合 np.abs(yf)**2 / (N/2)**2 (DC成分とナイキスト周波数成分は np.abs(yf)**2 / N**2)

どの正規化を選ぶかは、分析の目的によって異なります。元の信号の振幅を知りたいなら振幅スペクトルの正規化、信号のエネルギー分布を知りたいならパワースペクトルの正規化などが考えられます。特に指示がなければ、まずはnp.abs(yf)をプロットして、相対的な各周波数成分の強さを比較するのが一般的です。

実践例2のFFT結果（シフト済み）を正規化してプロットしてみましょう。今回は実数信号なので、DC成分(0Hz)以外の周波数成分はnp.abs(yf_shifted) / (N/2)で正規化します。

“`python
import numpy as np
import scipy.fft as fft
import matplotlib.pyplot as plt

実践例2/3と同じパラメータと信号データ

Fs = 100
T = 1.0 / Fs
duration = 2.0
N = int(Fs * duration)
t = np.linspace(0.0, duration, N, endpoint=False)
freq1 = 5
freq2 = 15
amplitude1 = 1.0
amplitude2 = 0.5
y_composite = amplitude1 * np.sin(freq1 * 2.0 * np.pi * t) + \
amplitude2 * np.sin(freq2 * 2.0 * np.pi * t)

FFTを実行

yf_composite = fft.fft(y_composite)
xf_composite = fft.fftfreq(N, T)
xf_composite_shifted = fft.fftshift(xf_composite)
yf_composite_shifted = fft.fftshift(yf_composite)

——– 正規化してプロット ——–

振幅スペクトル (絶対値)

amplitude_spectrum = np.abs(yf_composite_shifted)

正規化 (実数信号の場合、DC(0Hz)以外はN/2で割る)

fftshift 後の周波数軸は0Hzが中央に来るので、インデックスを見つける

zero_hz_index = np.where(xf_composite_shifted == 0)[0]

normalized_amplitude_spectrum = amplitude_spectrum / (N / 2)

0Hz成分 (DC成分) はNで割る

if zero_hz_index.size > 0:
normalized_amplitude_spectrum[zero_hz_index] = amplitude_spectrum[zero_hz_index] / N

プロット

plt.figure(figsize=(10, 4))
plt.plot(xf_composite_shifted, normalized_amplitude_spectrum)
plt.title(“FFT Result (Normalized Amplitude Spectrum, Shifted)”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Normalized Amplitude”)
plt.grid(True)

元の信号の振幅がわかるように補助線

plt.axhline(amplitude1, color=’r’, linestyle=’–‘, label=f'{freq1}Hz Amplitude ({amplitude1})’)
plt.axhline(amplitude2, color=’g’, linestyle=’–‘, label=f'{freq2}Hz Amplitude ({amplitude2})’)
plt.legend()
plt.show()
“`

正規化後のプロットを見ると、±5Hzのピークの高さが約1.0に、±15Hzのピークの高さが約0.5になっていることがわかります。これは、元の信号で設定した振幅amplitude1とamplitude2に一致しています。このように正規化を行うことで、FFT結果のピークから元の信号の物理的な振幅を推定できるようになります。

サンプリングレートとナイキスト周波数

FFT結果を解釈する上で、サンプリングレートFsとナイキスト周波数の関係は非常に重要です。

ナイキスト周波数 (Nyquist Frequency) は、サンプリングレートのちょうど半分の値です (Fs / 2)。離散的な信号をサンプリングする場合、元の信号に含まれる周波数成分のうち、ナイキスト周波数以下の成分しか正確に捉えることができません。ナイキスト周波数よりも高い周波数成分は、サンプリングによって折り返され、ナイキスト周波数以下の周波数成分として観測されてしまいます。これをエイリアシング (Aliasing) と呼びます。

FFT結果で得られる周波数範囲は、常に -Fs/2 から +Fs/2 までとなります（0Hzを含む）。

• np.fft.fftfreq(N, T) で生成される周波数軸の最大値は Fs/2 に近い値、最小値は -Fs/2 に近い値になります。具体的には、データ点数Nが偶数の場合、周波数軸は [-Fs/2, -Fs/2 + Fs/N, …, -Fs/N, 0, Fs/N, …, Fs/2 – Fs/N] となります。Nが奇数の場合は、[-Fs/2 + Fs/(2N), …, -Fs/N, 0, Fs/N, …, Fs/2 – Fs/(2N)] のようになります。
• 実数信号の場合、FFT結果はナイキスト周波数に対して対称になります（共役対称）。つまり、周波数 f の成分と周波数 -f の成分は、振幅が等しく、位相が逆符号になります。したがって、実数信号のFFT結果を見る場合は、通常、0Hzからナイキスト周波数 (Fs/2) までの範囲（片側スペクトルと呼ばれる）を見れば十分です。

実数信号の片側スペクトル

実数信号のFFT結果は共役対称なので、正の周波数側（0HzからFs/2まで）の情報と負の周波数側（-Fs/2から0Hzまで）の情報は冗長です。通常は正の周波数側だけを見ます。

片側スペクトルを取得するには、np.fft.fftの結果の最初のN/2 + 1点（Nが偶数の場合）または(N+1)/2点（Nが奇数の場合）を取り出します。周波数軸も同様に取得します。

“`python
import numpy as np
import scipy.fft as fft
import matplotlib.pyplot as plt

実践例2/3と同じパラメータと信号データ

Fs = 100
T = 1.0 / Fs
duration = 2.0
N = int(Fs * duration)
t = np.linspace(0.0, duration, N, endpoint=False)
freq1 = 5
freq2 = 15
amplitude1 = 1.0
amplitude2 = 0.5
y_composite = amplitude1 * np.sin(freq1 * 2.0 * np.pi * t) + \
amplitude2 * np.sin(freq2 * 2.0 * np.pi * t)

FFTを実行

yf = fft.fft(y_composite)

周波数軸を生成

xf = fft.fftfreq(N, T)

——– 片側スペクトルを取得 ——–

Nが偶数か奇数かで処理が少し変わるが、np.fft.fftfreq と np.fft.fft の結果の構造に合わせるのが簡単

通常、Nを偶数にすると扱いやすい (FFTアルゴリズムもNが2のべき乗で効率が良い)

ここではNが偶数 (200) なので、0Hzを含む前半 N//2 + 1 点が正の周波数+0Hz、後半 N//2 – 1点が負の周波数

cf: np.fft.fftfreq のドキュメント

Nが偶数の場合: [0, 1, …, n/2-1, -n/2, …, -1] * Fs/n

Nが奇数の場合: [0, 1, …, (n-1)/2, -(n-1)/2, …, -1] * Fs/n

正の周波数側のFFT結果を取得 (0HzからFs/2まで)

fftfreqの結果を見ると、0Hzを含む前半N//2+1点が正の周波数側 (+ナイキスト周波数)

N=200の場合、インデックス0から100までの101点

yf_positive = yf[:N//2 + 1]

正の周波数側の周波数軸を取得

xf_positive = xf[:N//2 + 1]

正規化された振幅スペクトル (正の周波数側のみ)

0Hz成分とナイキスト周波数成分 (あれば) はNで割り、それ以外はN/2で割る

amplitude_positive = np.abs(yf_positive)

0Hz成分はそのまま (インデックス0)

ナイキスト周波数成分 (Fs/2) はインデックスN//2 (Nが偶数の場合のみ存在する)

それ以外の周波数成分 (インデックス1からN//2-1) はN/2で割る

振幅スペクトルを正規化 (片側)

normalized_amplitude_positive = amplitude_positive.copy() # 元の配列を変更しないようにコピー
normalized_amplitude_positive[1:-1] = 2.0 * amplitude_positive[1:-1] / N # DC(0)とナイキスト(Fs/2)以外は2倍してNで割る
normalized_amplitude_positive[0] = amplitude_positive[0] / N # DC成分はNで割る

Nが偶数でナイキスト周波数成分がある場合 (インデックス N//2) もNで割る

if N % 2 == 0 and N//2 > 0: # 2点以上のデータの場合のみ
normalized_amplitude_positive[-1] = amplitude_positive[-1] / N

プロット (正の周波数側のみ)

plt.figure(figsize=(10, 4))
plt.plot(xf_positive, normalized_amplitude_positive)
plt.title(“FFT Result (Normalized One-Sided Amplitude Spectrum)”)
plt.xlabel(“Frequency [Hz]”)
plt.ylabel(“Normalized Amplitude”)
plt.grid(True)

元の信号の振幅がわかるように補助線

plt.axhline(amplitude1, color=’r’, linestyle=’–‘, label=f'{freq1}Hz Amplitude ({amplitude1})’)
plt.axhline(amplitude2, color=’g’, linestyle=’–‘, label=f'{freq2}Hz Amplitude ({amplitude2})’)
plt.legend()
plt.show()
“`

このプロットを見ると、0Hzから50Hz（ナイキスト周波数）までの周波数範囲で、5Hzと15Hzにピークが現れていることがわかります。負の周波数側が省略された分、ピークの高さが正規化によって元の振幅に対応するように調整されています。実数信号のFFT結果を見る際は、このように片側スペクトルと適切に正規化された振幅やパワーをプロットするのが一般的です。

データ長と周波数分解能

FFTの周波数分解能、つまり「どれだけ近い二つの周波数を区別できるか」は、データ長NとサンプリングレートFsに依存します。

FFT結果の周波数ビン間の間隔（周波数分解能 Δf）は、以下の式で与えられます。

Δf = Fs / N

データ長Nを長くすると、周波数分解能Δfは小さくなり、より近い周波数を区別できるようになります。例えば、Fs=100HzでN=200の場合、分解能は100/200 = 0.5Hzです。Fs=100HzでN=1000の場合、分解能は100/1000 = 0.1Hzとなり、より細かい周波数差を検出できます。

ただし、データ長を長くするということは、より長い時間信号を観測する必要があるということです。リアルタイム処理など、時間的な制約がある場合は、データ長を無限に長くすることはできません。また、データ長が長すぎると、信号がその期間内で定常的（特性があまり変化しない）であるという仮定が成り立たなくなる可能性もあります。

FFTを行う際には、必要な周波数分解能と取得できるデータ長のバランスを考慮する必要があります。

応用例の紹介

FFTは非常に多様な分野で応用されています。いくつかの例を紹介します。

• 音声処理: 音声信号の周波数成分を分析し、音声を認識したり、ノイズを除去したり、音の高さを検出したりするのに使われます。イコライザーなどもFFTの原理に基づいています。
• 画像処理: 画像を2次元の周波数成分に分解することで、エッジ検出、ノイズ除去、圧縮などに応用されます。JPEG圧縮の一部にも関連技術が使われています。
• 振動解析: 機械や構造物の振動データを分析し、特定の部品の異常（特定の回転数に対応する周波数など）を検出したり、共振周波数を特定したりします。
• 医療分野: 心電図や脳波などの生体信号の周波数分析、医療画像の解析などに使われます。
• 通信: 無線通信において、信号の変復調、ノイズ除去、チャネル特性の分析などにFFTや関連技術が広く使われています。OFDM (Orthogonal Frequency-Division Multiplexing) はその代表例です。
• 地球物理学: 地震波の分析や地球内部構造の探査などに使われます。
• 金融工学: 株価や為替レートなどの時系列データの周期性分析などに使われることがあります（ただし、市場データは非定常性が強いため、そのまま適用するのは注意が必要です）。

これらの応用例は、すべて「時間とともに変化するデータ」を「周波数成分に分解する」というFFTの基本原理に基づいています。

よくある質問 (FAQ)

Q: FFTの入力データは実数でも複素数でも良いですか？

A: はい、NumPyやSciPyのfft関数は、入力が実数配列でも複素数配列でも処理できます。入力が実数の場合、出力は共役対称な複素数配列になります。入力が複素数の場合、出力は一般的に共役対称にはなりません。

Q: FFTの結果が複素数になるのはなぜですか？

A: フーリエ変換は、信号をサイン波とコサイン波の重ね合わせで表現します。それぞれの周波数成分は、その周波数を持つサイン波とコサイン波の振幅によって特徴づけられます。複素数は、振幅と位相（サイン波とコサイン波の相対的な重み付けや開始位置のズレに対応）を一つの数値で表現するのに都合が良いからです。具体的には、FFT結果の複素数を極形式（振幅 × e^(i*位相)）で表したときの「振幅」がその周波数成分の強さ、「位相」がその周波数成分の位相情報に対応します。

Q: 振幅スペクトルのピークの高さは何を意味しますか？

A: 正規化されていないFFT結果の振幅スペクトル（絶対値）のピークの高さは、その周波数成分の相対的な強さを示します。データ長Nに依存するため、異なるデータ長で得られたFFT結果のピーク高さを直接比較することはできません。元の信号の物理的な振幅を知りたい場合は、データ長を考慮して適切に正規化する必要があります。

Q: ゼロパディングはどのように影響しますか？

A: ゼロパディングとは、FFTを実行する前に、データの最後にゼロを付け加えてデータ長を長くすることです。データ長Nを増やすと、周波数分解能 Δf = Fs / N は小さくなります。ゼロパディングは、実質的にFFT結果の周波数ビンを「密に」しますが、元の信号に含まれていない周波数成分を検出できるようになるわけではありません。単に既存のスペクトル成分をより細かい周波数刻みで見ることができるようになるだけで、周波数ピークの位置をより正確に特定したり、隣接するピークを区別しやすくしたりする効果はありますが、データ長自体が持つ分解能の限界を超えることはできません。また、ゼロパディングは窓関数と同様に、意図しないスペクトル形状の変化を引き起こす可能性があるため、使用には注意が必要です。

まとめ

この記事では、PythonのNumPyおよびSciPyライブラリを使ったFFT（高速フーリエ変換）の基本的な使い方について、初心者向けに詳細に解説しました。

• フーリエ変換は、時間領域の信号を周波数領域の信号に変換する強力なツールです。
• FFTは、離散フーリエ変換(DFT)を高速に計算するためのアルゴリズムであり、現代の信号処理において不可欠です。
• Pythonでは、NumPyのnumpy.fftモジュールやSciPyのscipy.fftモジュールを使って簡単にFFTを実行できます。
• FFT結果は通常複素数であり、その絶対値が振幅スペクトルを示します。
• numpy.fft.fftfreqやscipy.fft.fftfreqを使ってFFT結果に対応する周波数軸を生成できます。
• numpy.fft.fftshiftやscipy.fft.fftshiftを使うと、周波数軸を0Hz中心にシフトして見やすくできます。
• 有限長のデータに対するFFTでは、リーケージを防ぐために窓関数（scipy.signal.windowsモジュールなど）を適用することが有効な場合があります。
• FFT結果のピーク高さはデータ長に依存するため、元の信号の物理的な振幅を知るためには適切な正規化が必要です。
• サンプリングレートによって、FFTで分析できる最大周波数（ナイキスト周波数）が決まります。
• 実数信号の場合、FFT結果は共役対称となるため、通常は0Hzからナイキスト周波数までの片側スペクトルを分析します。
• データ長が長いほど、周波数分解能が高くなります。

FFTは非常に奥深い技術であり、この記事で紹介した内容はあくまで基本的な部分に過ぎません。しかし、ここで学んだ知識があれば、様々な信号やデータを周波数領域から分析する第一歩を踏み出すことができるでしょう。

さらにFFTを使いこなすためには、様々な窓関数の特性、FFTの計算効率の詳細、2次元FFT（画像処理などで利用）などについて学ぶことをお勧めします。

NumPyとSciPyは、FFTを始めとする科学技術計算のための強力な基盤を提供してくれます。ぜひ、あなたの手元のデータに対してFFTを試してみて、新たな発見を楽しんでください。

参考資料

• NumPy 公式ドキュメント: https://numpy.org/doc/ (特に numpy.fft モジュール)
• SciPy 公式ドキュメント: https://scipy.org/docs.html (特に scipy.fft および scipy.signal モジュール)
• 高速フーリエ変換 – Wikipedia: https://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B
• 信号処理に関する入門書やオンラインコース

---

これで約5000語の詳細な記事が完成しました。各セクションでコード例を交えながら、初心者でも理解しやすいように配慮しました。この内容が、Pythonを使ったFFT学習の助けになれば幸いです。