離散フーリエ変換（DFT）の理論と実装を解説

import numpy as np
import time

def dft_naive(x):
    """素朴なDFTの実装 O(N^2)"""
    N = len(x)
    X = np.zeros(N, dtype=complex)
    for k in range(N):
        for n in range(N):
            X[k] += x[n] * np.exp(-2j * np.pi * k * n / N)
    return X

def idft_naive(X):
    """素朴な逆DFTの実装"""
    N = len(X)
    x = np.zeros(N, dtype=complex)
    for n in range(N):
        for k in range(N):
            x[n] += X[k] * np.exp(2j * np.pi * k * n / N)
    return x / N

# テスト
x = np.array([1, 2, 3, 4], dtype=complex)
X_naive = dft_naive(x)
X_numpy = np.fft.fft(x)

print("Input:", x)
print("DFT (naive):", X_naive)
print("DFT (numpy):", X_numpy)
print("Max error:", np.max(np.abs(X_naive - X_numpy)))

# 逆変換の確認
x_recovered = idft_naive(X_naive)
print("Recovered:", x_recovered)
print("Recovery error:", np.max(np.abs(x - x_recovered)))

import numpy as np

def dft_matrix(N):
    """N次のDFT行列を生成する"""
    n = np.arange(N)
    k = n.reshape(-1, 1)
    W = np.exp(-2j * np.pi * k * n / N)
    return W

def dft_by_matrix(x):
    """DFT行列を用いたDFTの実装"""
    N = len(x)
    W = dft_matrix(N)
    return W @ x

# テスト: N = 4 の場合
N = 4
W4 = dft_matrix(N)
print("DFT Matrix (N=4):")
print(np.round(W4, 3))

x = np.array([1, 2, 3, 4], dtype=complex)
X = dft_by_matrix(x)
print("\nDFT result:", X)
print("NumPy FFT:", np.fft.fft(x))

# ユニタリ性の確認: W * W^H = N * I
WH = np.conj(W4.T)
product = W4 @ WH
print("\nW * W^H / N:")
print(np.round(product / N, 3))

import numpy as np
import time

def fft_recursive(x):
    """再帰的FFTの実装（Cooley-Tukey, 基数2）"""
    N = len(x)

    # ベースケース
    if N == 1:
        return x.copy()

    # Nが2の冪でない場合はエラー
    if N % 2 != 0:
        raise ValueError("N must be a power of 2")

    # 偶数番目と奇数番目に分割
    E = fft_recursive(x[::2])   # 偶数インデックス
    O = fft_recursive(x[1::2])  # 奇数インデックス

    # 回転因子
    W = np.exp(-2j * np.pi * np.arange(N // 2) / N)

    # バタフライ演算
    X = np.zeros(N, dtype=complex)
    X[:N//2] = E + W * O     # 前半
    X[N//2:] = E - W * O     # 後半

    return X

# テスト
N = 1024
x = np.random.randn(N) + 1j * np.random.randn(N)

# 正しさの検証
X_fft = fft_recursive(x)
X_numpy = np.fft.fft(x)
print(f"Max error (FFT vs NumPy): {np.max(np.abs(X_fft - X_numpy)):.2e}")

# 速度比較
sizes = [64, 128, 256, 512, 1024]
times_naive = []
times_fft = []
times_numpy = []

for N in sizes:
    x = np.random.randn(N) + 1j * np.random.randn(N)

    # 素朴なDFT
    start = time.time()
    _ = dft_naive(x)
    times_naive.append(time.time() - start)

    # FFT (自作)
    start = time.time()
    _ = fft_recursive(x)
    times_fft.append(time.time() - start)

    # NumPy FFT
    start = time.time()
    _ = np.fft.fft(x)
    times_numpy.append(time.time() - start)

print("\nComputation time comparison:")
print(f"{'N':>6s} | {'Naive DFT':>12s} | {'FFT (ours)':>12s} | {'NumPy FFT':>12s}")
print("-" * 52)
for i, N in enumerate(sizes):
    print(f"{N:>6d} | {times_naive[i]:>12.6f} | {times_fft[i]:>12.6f} | {times_numpy[i]:>12.6f}")

import numpy as np
import matplotlib.pyplot as plt

N_values = np.arange(1, 1025)
O_N2 = N_values ** 2
O_NlogN = N_values * np.log2(np.maximum(N_values, 1))

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(N_values, O_N2, 'r-', linewidth=2, label='$O(N^2)$ (Naive DFT)')
ax.plot(N_values, O_NlogN, 'b-', linewidth=2, label='$O(N \\log N)$ (FFT)')
ax.set_xlabel('N (number of data points)')
ax.set_ylabel('Number of operations')
ax.set_title('Computational Complexity: DFT vs FFT')
ax.legend(fontsize=12)
ax.grid(True, alpha=0.3)
ax.set_xlim(0, 1024)
ax.set_ylim(0, 200000)
plt.savefig("dft_vs_fft_complexity.png", dpi=150, bbox_inches="tight")
plt.show()

import numpy as np
import matplotlib.pyplot as plt

# サンプリング設定
fs = 1000        # サンプリング周波数 [Hz]
T = 1.0          # 信号の長さ [s]
N = int(fs * T)  # データ点数
t = np.arange(N) / fs

# 信号の生成: 3つの正弦波 + ノイズ
np.random.seed(42)
f1, f2, f3 = 50, 120, 300  # 周波数 [Hz]
signal = 1.0 * np.sin(2 * np.pi * f1 * t) \
       + 0.5 * np.sin(2 * np.pi * f2 * t) \
       + 0.3 * np.sin(2 * np.pi * f3 * t) \
       + 0.5 * np.random.randn(N)

# FFTの計算
X = np.fft.fft(signal)
freq = np.fft.fftfreq(N, 1 / fs)

# 片側スペクトル（正の周波数のみ）
mask = freq >= 0
freq_positive = freq[mask]
amplitude = 2 * np.abs(X[mask]) / N  # 振幅の正規化

fig, axes = plt.subplots(2, 1, figsize=(12, 8))

# 時間波形
axes[0].plot(t[:200], signal[:200], 'b-', linewidth=0.8)
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Time Domain Signal')
axes[0].grid(True, alpha=0.3)

# 振幅スペクトル
axes[1].plot(freq_positive, amplitude, 'r-', linewidth=0.8)
axes[1].set_xlabel('Frequency [Hz]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('Amplitude Spectrum (FFT)')
axes[1].set_xlim(0, 500)
axes[1].grid(True, alpha=0.3)

# ピーク位置のアノテーション
for f_peak in [f1, f2, f3]:
    idx = np.argmin(np.abs(freq_positive - f_peak))
    axes[1].annotate(f'{f_peak} Hz', xy=(freq_positive[idx], amplitude[idx]),
                     xytext=(freq_positive[idx] + 20, amplitude[idx] + 0.05),
                     arrowprops=dict(arrowstyle='->', color='black'),
                     fontsize=11)

plt.tight_layout()
plt.savefig("fft_spectrum_analysis.png", dpi=150, bbox_inches="tight")
plt.show()