NumPyの離散フーリエ変換：fft.rfft() 関数とは？

NumPy の離散フーリエ変換: fft.rfft() の詳細解説

NumPy の fft.rfft() 関数は、実数配列の離散フーリエ変換 (DFT) を計算します。これは、フーリエ変換の高速アルゴリズムである Cooley-Tukey アルゴリズム を用いて効率的に計算されます。

用途

fft.rfft() は、音声処理、画像処理、信号処理など、さまざまな分野で広く使用されています。具体的には、以下のような用途があります。

• 周波数分析: 音声や信号の周波数成分を分析
• フィルタリング: 特定の周波数成分を除去したり、強調したりする
• スペクトル推定: 信号のスペクトルを推定

入力と出力

fft.rfft() 関数は、以下の入力と出力を持ちます。

• 入力: 1次元または多次元の numpy.ndarray 型の実数配列
• 出力: 入力配列と同じ形状の numpy.ndarray 型の複素数配列

出力データの解釈

fft.rfft() 関数の出力データは、以下の形式で解釈できます。

• 出力データの最初の要素は、入力配列の直流成分 (DC成分) を表します。
• 出力データの残りの要素は、負の周波数成分と正の周波数成分が交互に並んでいます。
• 負の周波数成分は、正の周波数成分と対称的な関係にあります。

以下のサンプルコードは、fft.rfft() 関数を使用して、1次元配列の DFT を計算し、結果をグラフで表示します。

import numpy as np
import matplotlib.pyplot as plt

# 入力データ
x = np.linspace(0, 2*np.pi, 100)
y = np.sin(x)

# DFT を計算
yf = np.fft.rfft(y)

# 周波数軸
freq = np.fft.rfftfreq(x.size)

# 結果を表示
plt.plot(freq, np.abs(yf))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()

fft.rfft() 関数に関する詳細は、以下の NumPy ドキュメントを参照してください。

補足

• fft.rfft() 関数は、fft.fft() 関数と比べて、以下の利点があります。
  • 入力データが実数の場合、計算量が半分になる。
  • 出力データが半分になるため、メモリ使用量が少なくて済む。
• fft.rfft() 関数は、ifft.irfft() 関数と組み合わせて、逆フーリエ変換を計算することができます。

fft.rfft() 関数のサンプルコード

このサンプルコードは、fft.rfft() 関数を使用して、音声データの周波数分析を行います。

import numpy as np
import matplotlib.pyplot as plt
import soundfile as sf

# 音声データを読み込み
data, samplerate = sf.read("audio.wav")

# DFT を計算
yf = np.fft.rfft(data)

# 周波数軸
freq = np.fft.rfftfreq(data.size, samplerate)

# 結果を表示
plt.plot(freq, np.abs(yf))
plt.xlabel("Frequency (Hz)")
plt.ylabel("Magnitude")
plt.show()

画像データのフィルタリング

このサンプルコードは、fft.rfft() 関数を使用して、画像データの特定の周波数成分を除去するフィルタリングを行います。

import numpy as np
import matplotlib.pyplot as plt
from PIL import Image

# 画像データを読み込み
img = Image.open("image.png").convert("L")
data = np.array(img)

# DFT を計算
yf = np.fft.rfft2(data)

# 低周波成分のみを残す
yf_filtered = np.copy(yf)
yf_filtered[freq > cutoff_freq] = 0

# 逆フーリエ変換を計算
data_filtered = np.fft.irfft2(yf_filtered)

# 結果を表示
plt.subplot(121)
plt.imshow(data, cmap="gray")
plt.title("Original")
plt.subplot(122)
plt.imshow(data_filtered, cmap="gray")
plt.title("Filtered")
plt.show()

信号処理におけるスペクトル推定

このサンプルコードは、fft.rfft() 関数を使用して、信号処理におけるスペクトル推定を行います。

import numpy as np
import matplotlib.pyplot as plt

# 信号データ
x = np.sin(2*np.pi*10*t) + np.random.randn(t.size)

# DFT を計算
yf = np.fft.rfft(x)

# 周波数軸
freq = np.fft.rfftfreq(t.size)

# 平均スペクトル
psd = np.abs(yf)**2 / t.size

# 結果を表示
plt.plot(freq, psd)
plt.xlabel("Frequency (Hz)")
plt.ylabel("Power Spectral Density")
plt.show()

その他

上記以外にも、fft.rfft() 関数はさまざまな用途で使用できます。詳細は、NumPy ドキュメントやチュートリアルを参照してください。

fft.rfft() 関数の代替方法

fft.fft() 関数は、複素数配列の DFT を計算します。実数配列の DFT を計算するためには、まず入力配列を複素数配列に変換する必要があります。

import numpy as np

# 入力データ
x = np.array([1, 2, 3, 4])

# 複素数配列に変換
x_complex = np.zeros(x.size, dtype=complex)
x_complex.real = x

# DFT を計算
yf = np.fft.fft(x_complex)

直接計算

DFT は直接計算することもできます。ただし、この方法は計算量が非常に多いため、大きな配列には適していません。

import numpy as np

def DFT(x):
  N = x.size
  y = np.zeros(N, dtype=complex)
  for k in range(N):
    for n in range(N):
      y[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
  return y

# 入力データ
x = np.array([1, 2, 3, 4])

# DFT を計算
yf = DFT(x)

その他のライブラリ

NumPy 以外にも、SciPy や PyTorch などのライブラリで DFT を計算する関数を提供しています。

どの方法を選択するかは、計算量、精度、使いやすさなどの要件によって異なります。