🔷 矩形波とは？
矩形波（square wave）は、ある周期で「+A」と「–A」を交互に取る周期的な非正弦波です。
たとえば周期 $T$、振幅 $A$、デューティ比50%の対称な波を考えます。
🔷 フーリエ級数の基本形
任意の周期関数 $f(t)$ は次のように展開できます：
$$
f(t) = \frac{a_0}{2} + \sum_{n=1}^{\infty} \left[ a_n \cos\left( \frac{2\pi n t}{T} \right) + b_n \sin\left( \frac{2\pi n t}{T} \right) \right]
$$
矩形波のような奇関数（左右対称に反転）では、余弦項 $a_n$ はすべて 0 になります。
したがって、矩形波のフーリエ級数は次のような正弦波の無限和になります。
🔷 矩形波のフーリエ級数（デューティ50%）
振幅 $A$、周期 $T$、中心が 0 の対称な矩形波は以下のように表せます：
$$
f(t) = \frac{4A}{\pi} \sum_{n=1,3,5,\dots}^{\infty} \frac{1}{n} \sin\left( \frac{2\pi n t}{T} \right)
$$
または、角周波数 $\omega = \frac{2\pi}{T}$ を使うと：
$$
f(t) = \frac{4A}{\pi} \sum_{k=0}^{\infty} \frac{1}{2k+1} \sin\left( (2k+1)\omega t \right)
$$
🔷 各成分の意味
- 基本周波数 $f = \frac{1}{T}$：1つ目の成分（n=1）
- 高調波（odd harmonics）：n = 3, 5, 7, ...（奇数倍の周波数）
- 各高調波の振幅は $\frac{1}{n}$ に比例して小さくなる
- 偶数次成分（n = 2, 4, ...）は含まれない
🔷 グラフでの意味（ギブズ現象）
このフーリエ級数を何項か足し合わせていくと、矩形波に近づきますが、角のところでオーバーシュート（ギブズ現象）が残るのが特徴です。
# Program Name: square_wave_fft_N100.py
# Creation Date: 20250801
# Overview: Generate square wave using Fourier series (N=100 terms), then perform FFT and plot results
# Usage: Run to construct square wave (N=100), perform FFT, and visualize time/frequency domains
!pip install numpy matplotlib scipy
import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq
# -----------------------------
# パラメータ設定 / Parameters
# -----------------------------
A = 1.0 # 振幅 / Amplitude
T = 1.0 # 周期 / Period [sec]
f0 = 1.0 / T # 基本周波数 / Fundamental frequency [Hz]
fs = 1000 # サンプリング周波数 / Sampling rate [Hz]
N_samples = 1000 # サンプル数 / Number of time-domain samples
N_terms = 100 # フーリエ級数の項数 / Number of Fourier terms (odd harmonics only)
t = np.linspace(0, T, N_samples, endpoint=False) # 時間軸 / Time axis
omega = 2 * np.pi / T # 角周波数 / Angular frequency
# -----------------------------
# フーリエ級数で矩形波生成 / Generate square wave from Fourier series
# -----------------------------
square_wave = np.zeros_like(t) # 初期化 / initialize
for k in range(N_terms):
n = 2 * k + 1 # 奇数次のみ / odd harmonics only
square_wave += (1 / n) * np.sin(n * omega * t)
square_wave *= (4 * A / np.pi) # スケーリング / amplitude scale
# -----------------------------
# FFT解析 / Perform FFT
# -----------------------------
yf = fft(square_wave) # 複素スペクトル / Complex spectrum
xf = fftfreq(N_samples, 1/fs) # 周波数軸 / Frequency axis
amplitude = 2.0 / N_samples * np.abs(yf) # 振幅スペクトル / Amplitude spectrum
# -----------------------------
# 可視化 / Visualization
# -----------------------------
plt.figure(figsize=(12, 5))
# 時間波形 / Time domain
plt.subplot(1, 2, 1)
plt.plot(t, square_wave)
plt.title("Square Wave by Fourier Series (N=100)")
plt.xlabel("Time [s]")
plt.ylabel("Amplitude")
plt.grid(True)
# 周波数スペクトル / Frequency domain
plt.subplot(1, 2, 2)
plt.stem(xf[:N_samples//2], amplitude[:N_samples//2], basefmt=" ")
plt.title("FFT Amplitude Spectrum")
plt.xlabel("Frequency [Hz]")
plt.ylabel("Amplitude")
plt.grid(True)
plt.tight_layout()
plt.show()