信号処理・制御システム完全ガイド: フーリエ変換からPID制御まで

• 1. 信号とシステムの基礎
  • 1.1 信号の分類
  • 1.2 基本信号
  • 1.3 線形時不変(LTI)システム
  • 1.4 畳み込み(Convolution)
• 2. フーリエ級数 (Fourier Series)
  • 2.1 周期信号のフーリエ級数
  • 2.2 矩形波(Square Wave)の例
  • 2.3 パーセバルの定理
• 3. フーリエ変換 (Fourier Transform)
  • 3.1 連続時間フーリエ変換 (CTFT)
  • 3.2 重要な変換ペア
  • 3.3 フーリエ変換の性質
  • 3.4 PythonによるFFT計算
• 4. 離散フーリエ変換 (DFT) と FFT
  • 4.1 離散時間フーリエ変換 (DTFT)
  • 4.2 DFTの定義
  • 4.3 FFTアルゴリズム
  • 4.4 スペクトル解析の例
  • 4.5 ナイキスト・シャノンのサンプリング定理
• 5. ラプラス変換 (Laplace Transform)
  • 5.1 定義
  • 5.2 重要な変換ペア
  • 5.3 伝達関数 (Transfer Function)
  • 5.4 2次標準系
• 6. Z変換
  • 6.1 定義
  • 6.2 Z変換とラプラス変換の関係
  • 6.3 安定性解析
  • 6.4 重要なZ変換ペア
• 7. デジタルフィルタ設計
  • 7.1 FIRフィルタ vs IIRフィルタ
  • 7.2 フィルタの種類
  • 7.3 古典的フィルタの種類
  • 7.4 scipy.signalによるフィルタ設計
  • 7.5 FIRフィルタ設計 (窓関数法)
• 8. 制御システムの基礎
  • 8.1 開ループ vs 閉ループ制御
  • 8.2 安定性: ラウス・フルビッツ基準
  • 8.3 ボード線図 (Bode Plot)
• 9. PID制御器
  • 9.1 PID制御動作
  • 9.2 ジーグラー・ニコルス調整法
  • 9.3 PythonによるPIDシミュレーション
• 10. 状態空間表現 (State Space)
  • 10.1 状態方程式
  • 10.2 可制御性と可観測性
  • 10.3 線形二次レギュレータ (LQR)
• 11. AI/MLと信号処理
  • 11.1 1D CNNによる信号分類
  • 11.2 STFTとメルスペクトログラム
  • 11.3 信号処理 + ディープラーニングの応用例
• 12. 確認問題
• 参考文献

1. 信号とシステムの基礎

1.1 信号の分類

信号(Signal)とは情報を伝達する物理量の変化です。主な分類を以下に示します。

連続時間信号 vs 離散時間信号

連続時間(CT)信号はすべての $t \in \mathbb{R}$ で定義されます。例: アナログ音声信号 $x(t) = A\cos(\omega_0 t + \phi)$

離散時間(DT)信号は整数インデックス $n \in \mathbb{Z}$ でのみ定義されます。例: デジタル音声サンプル列 $x[n]$

エネルギー信号と電力信号

信号のエネルギー $E$ と平均電力 $P$:

$E = \int_{-\infty}^{\infty} |x(t)|^2 \, dt$

$P = \lim_{T \to \infty} \frac{1}{2T} \int_{-T}^{T} |x(t)|^2 \, dt$

$E < \infty$ ならエネルギー信号、$0 < P < \infty$ なら電力信号です。

1.2 基本信号

単位インパルス(Unit Impulse)

$\delta(t) = \begin{cases} \infty & t = 0 \\ 0 & t \neq 0 \end{cases}, \quad \int_{-\infty}^{\infty} \delta(t) \, dt = 1$

サンプリング性質: $\int_{-\infty}^{\infty} x(t)\delta(t - t_0) \, dt = x(t_0)$

単位ステップ(Unit Step)

$u(t) = \begin{cases} 1 & t \geq 0 \\ 0 & t < 0 \end{cases}$

関係式 $u(t) = \int_{-\infty}^{t} \delta(\tau) \, d\tau$ が成立します。

1.3 線形時不変(LTI)システム

LTIシステムは2つの核心的な性質を満たします。

• 線形性(Linearity): $\mathcal{T}\{ax_1(t) + bx_2(t)\} = a\mathcal{T}\{x_1(t)\} + b\mathcal{T}\{x_2(t)\}$
• 時不変性(Time Invariance): $y(t) = \mathcal{T}\{x(t)\} \Rightarrow \mathcal{T}\{x(t-t_0)\} = y(t-t_0)$

1.4 畳み込み(Convolution)

LTIシステムの出力は入力とインパルス応答の畳み込みで表されます。

$y(t) = x(t) * h(t) = \int_{-\infty}^{\infty} x(\tau) h(t - \tau) \, d\tau$

離散時間の場合:

$y[n] = x[n] * h[n] = \sum_{k=-\infty}^{\infty} x[k] h[n - k]$

畳み込みは交換律・結合律・分配律を満たします。

2. フーリエ級数 (Fourier Series)

2.1 周期信号のフーリエ級数

周期 $T_0$ を持つ信号 $x(t)$ は基本周波数 $\omega_0 = 2\pi/T_0$ の複素指数関数の和で表せます。

$x(t) = \sum_{k=-\infty}^{\infty} c_k e^{jk\omega_0 t}$

フーリエ係数:

$c_k = \frac{1}{T_0} \int_{T_0} x(t) e^{-jk\omega_0 t} \, dt$

2.2 矩形波(Square Wave)の例

振幅 A、周期 T、デューティサイクル 50% の矩形波のフーリエ係数:

$c_k = \begin{cases} A/2 & k = 0 \\ \frac{A}{k\pi} \sin\left(\frac{k\pi}{2}\right) & k \neq 0 \end{cases}$

奇数次高調波のみが存在し、高次になるほど振幅は減少します。不連続点付近ではギブス現象(Gibbs Phenomenon)が現れます。

2.3 パーセバルの定理

$\frac{1}{T_0} \int_{T_0} |x(t)|^2 \, dt = \sum_{k=-\infty}^{\infty} |c_k|^2$

時間領域の平均電力と周波数領域の係数の二乗和が等しくなります。

3. フーリエ変換 (Fourier Transform)

3.1 連続時間フーリエ変換 (CTFT)

非周期信号のフーリエ変換ペア:

$X(j\omega) = \int_{-\infty}^{\infty} x(t) e^{-j\omega t} \, dt \quad \text{(正変換)}$

$x(t) = \frac{1}{2\pi} \int_{-\infty}^{\infty} X(j\omega) e^{j\omega t} \, d\omega \quad \text{(逆変換)}$

3.2 重要な変換ペア

3.3 フーリエ変換の性質

• 線形性: $ax_1(t) + bx_2(t) \leftrightarrow aX_1(j\omega) + bX_2(j\omega)$
• 時間シフト: $x(t-t_0) \leftrightarrow e^{-j\omega t_0} X(j\omega)$
• 周波数シフト(変調): $e^{j\omega_0 t} x(t) \leftrightarrow X(j(\omega - \omega_0))$
• 畳み込み定理: $x(t)*h(t) \leftrightarrow X(j\omega)H(j\omega)$
• 乗算定理: $x(t)h(t) \leftrightarrow \frac{1}{2\pi}X(j\omega)*H(j\omega)$
• 微分: $\frac{dx}{dt} \leftrightarrow j\omega X(j\omega)$
• パーセバル: $\int|x(t)|^2 dt = \frac{1}{2\pi}\int|X(j\omega)|^2 d\omega$

3.4 PythonによるFFT計算

import numpy as np
from scipy.fft import fft, ifft, fftfreq, fftshift
import matplotlib.pyplot as plt

# 矩形波信号の生成
fs = 1000       # サンプリング周波数 (Hz)
T = 1.0         # 信号長 (秒)
t = np.linspace(0, T, int(fs * T), endpoint=False)

from scipy.signal import square
x = square(2 * np.pi * 10 * t)   # 10Hz 矩形波

# FFT計算
N = len(x)
X = fft(x)
freqs = fftfreq(N, 1/fs)

# 両側スペクトル (正規化)
magnitude = np.abs(fftshift(X)) / N
freq_shifted = fftshift(freqs)

plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(t[:100], x[:100])
plt.title('Square Wave (time domain)')
plt.xlabel('Time (s)')

plt.subplot(1, 2, 2)
plt.plot(freq_shifted, magnitude)
plt.xlim(-100, 100)
plt.title('Fourier Transform (frequency domain)')
plt.xlabel('Frequency (Hz)')
plt.tight_layout()
plt.show()

4. 離散フーリエ変換 (DFT) と FFT

4.1 離散時間フーリエ変換 (DTFT)

離散時間信号 $x[n]$ の DTFT:

$X(e^{j\omega}) = \sum_{n=-\infty}^{\infty} x[n] e^{-j\omega n}$

逆変換:

$x[n] = \frac{1}{2\pi} \int_{-\pi}^{\pi} X(e^{j\omega}) e^{j\omega n} \, d\omega$

DTFTは連続な周波数 $\omega$ を持ち、$2\pi$ 周期です。

4.2 DFTの定義

N点DFTはDTFTをN個の周波数で等間隔サンプリングしたものです。

$X[k] = \sum_{n=0}^{N-1} x[n] e^{-j2\pi kn/N}, \quad k = 0, 1, \ldots, N-1$

逆DFT:

$x[n] = \frac{1}{N} \sum_{k=0}^{N-1} X[k] e^{j2\pi kn/N}, \quad n = 0, 1, \ldots, N-1$

4.3 FFTアルゴリズム

直接DFT計算の計算量は $O(N^2)$ ですが、FFT(クーリー・テューキーアルゴリズム)は $N = 2^m$ のとき分割統治法で $O(N \log N)$ を実現します。

旋回因子(Twiddle Factor): $W_N^k = e^{-j2\pi k/N}$

4.4 スペクトル解析の例

import numpy as np
from scipy.fft import fft, fftfreq
import matplotlib.pyplot as plt

# 複合信号の生成: 50Hz + 120Hz
fs = 1000           # サンプリング周波数 (Hz)
N = fs              # サンプル数 (1秒)
t = np.linspace(0, 1, N, endpoint=False)

# ノイズ付き複合信号
signal = (np.sin(2 * np.pi * 50 * t)
          + 0.5 * np.sin(2 * np.pi * 120 * t)
          + 0.2 * np.random.randn(N))

# FFT計算
yf = fft(signal)
xf = fftfreq(N, 1/fs)

# 片側スペクトル (正の周波数のみ)
half = N // 2
amplitude = 2.0 / N * np.abs(yf[:half])
freq_pos = xf[:half]

plt.figure(figsize=(10, 4))
plt.plot(freq_pos, amplitude)
plt.xlabel('Frequency (Hz)')
plt.ylabel('Amplitude')
plt.title('Single-Sided Amplitude Spectrum')
plt.axvline(x=50,  color='r', linestyle='--', label='50 Hz')
plt.axvline(x=120, color='g', linestyle='--', label='120 Hz')
plt.legend()
plt.grid(True)
plt.show()

print(f"周波数分解能: {fs/N} Hz")
print(f"ナイキスト周波数: {fs/2} Hz")

4.5 ナイキスト・シャノンのサンプリング定理

エイリアシング(折り返し歪み)なしに信号を復元するには、サンプリング周波数 $f_s$ が信号の最大周波数 $f_{max}$ の2倍以上必要です。

$f_s \geq 2 f_{max} \quad \text{(ナイキスト条件)}$

5. ラプラス変換 (Laplace Transform)

5.1 定義

両側ラプラス変換:

$X(s) = \int_{-\infty}^{\infty} x(t) e^{-st} \, dt, \quad s = \sigma + j\omega$

片側ラプラス変換:

$X(s) = \int_{0^-}^{\infty} x(t) e^{-st} \, dt$

5.2 重要な変換ペア

5.3 伝達関数 (Transfer Function)

LTIシステムの入出力関係をs領域で表現:

$H(s) = \frac{Y(s)}{X(s)} = \frac{b_m s^m + b_{m-1} s^{m-1} + \cdots + b_0}{a_n s^n + a_{n-1} s^{n-1} + \cdots + a_0}$

極(Poles): $H(s) \to \infty$ となる $s$ の値(分母 = 0)

零点(Zeros): $H(s) = 0$ となる $s$ の値(分子 = 0)

すべての極がs平面の左半面(LHP)にあればBIBO安定です。

5.4 2次標準系

$H(s) = \frac{\omega_n^2}{s^2 + 2\zeta\omega_n s + \omega_n^2}$

• $\omega_n$: 固有角周波数(Natural Frequency)
• $\zeta$: 減衰比(Damping Ratio)
  • $\zeta < 1$: 不足減衰(Underdamped) — 振動応答
  • $\zeta = 1$: 臨界減衰(Critically Damped) — 最速無振動
  • $\zeta > 1$: 過減衰(Overdamped) — 緩慢無振動

6. Z変換

6.1 定義

離散時間信号 $x[n]$ のZ変換:

$X(z) = \sum_{n=-\infty}^{\infty} x[n] z^{-n}, \quad z \in \mathbb{C}$

$z = e^{j\omega}$ (単位円上)と置換するとDTFTに一致します。

6.2 Z変換とラプラス変換の関係

連続時間システムをサンプリング周期 $T_s$ で離散化すると:

$z = e^{sT_s}$

この写像によりs平面とz平面が対応します。

• s平面の虚軸 $\leftrightarrow$ z平面の単位円
• s平面の左半面 $\leftrightarrow$ z平面の単位円内部

6.3 安定性解析

離散時間システムが安定であるためには、すべての極が単位円の内部になければなりません。

$|z_i| < 1 \quad \text{すべての極に対して}$

6.4 重要なZ変換ペア

7. デジタルフィルタ設計

7.1 FIRフィルタ vs IIRフィルタ

FIR(有限インパルス応答)フィルタ

$y[n] = \sum_{k=0}^{M} b_k x[n-k]$

• インパルス応答が有限長
• 常に安定(フィードバックなし)
• 線形位相を実現可能
• 急峻な遮断特性には高次数が必要

IIR(無限インパルス応答)フィルタ

$y[n] = \sum_{k=0}^{M} b_k x[n-k] - \sum_{k=1}^{N} a_k y[n-k]$

• フィードバックを含み、インパルス応答は無限
• 低次数で急峻な遮断特性
• 位相が非線形(位相歪みの可能性)
• 設計不良の場合は不安定になりうる

7.2 フィルタの種類

• 低域通過フィルタ(LPF): 低周波を通過、高周波を遮断
• 高域通過フィルタ(HPF): 高周波を通過、低周波を遮断
• 帯域通過フィルタ(BPF): 特定の周波数帯域のみ通過
• 帯域阻止フィルタ(BEF/Notch): 特定の周波数帯域を遮断

7.3 古典的フィルタの種類

バタワース(Butterworth): 通過帯域で最大限平坦な応答

チェビシェフ I 型: 通過帯域にリップルを許容し、より急峻な遮断特性

楕円(Elliptic): 通過帯域と阻止帯域の両方にリップルを持ち、最小次数で最大減衰

7.4 scipy.signalによるフィルタ設計

import numpy as np
from scipy import signal
import matplotlib.pyplot as plt

fs = 1000  # サンプリング周波数 (Hz)

# バタワース低域通過フィルタ (4次、遮断周波数100Hz)
nyq = fs / 2
cutoff = 100 / nyq  # 正規化周波数 (0〜1)
b_butter, a_butter = signal.butter(N=4, Wn=cutoff, btype='low')

# チェビシェフ I 型 (4次、1dBリップル、遮断周波数100Hz)
b_cheby, a_cheby = signal.cheby1(N=4, rp=1, Wn=cutoff, btype='low')

# 周波数応答の計算
w, h_butter = signal.freqz(b_butter, a_butter, fs=fs)
w, h_cheby  = signal.freqz(b_cheby,  a_cheby,  fs=fs)

plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(w, 20 * np.log10(np.abs(h_butter)), label='Butterworth')
plt.plot(w, 20 * np.log10(np.abs(h_cheby)),  label='Chebyshev I')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Magnitude (dB)')
plt.title('Filter Frequency Response')
plt.legend()
plt.grid(True)
plt.ylim(-80, 5)

# ノイズ信号のフィルタリング
t = np.linspace(0, 1, fs, endpoint=False)
clean    = np.sin(2 * np.pi * 50 * t)
noisy    = clean + 0.5 * np.sin(2 * np.pi * 300 * t)
filtered = signal.filtfilt(b_butter, a_butter, noisy)

plt.subplot(1, 2, 2)
plt.plot(t[:200], noisy[:200],    alpha=0.5, label='Noisy')
plt.plot(t[:200], filtered[:200],            label='Filtered')
plt.xlabel('Time (s)')
plt.legend()
plt.title('Filtering Result')
plt.tight_layout()
plt.show()

7.5 FIRフィルタ設計 (窓関数法)

from scipy.signal import firwin, freqz

# Hamming窓を用いたFIR LPF
# タップ数 = 51、遮断周波数 = 100Hz
numtaps = 51
cutoff_hz = 100
b_fir = firwin(numtaps, cutoff_hz / (fs / 2), window='hamming')

w_fir, h_fir = freqz(b_fir, [1.0], fs=fs)
print(f"FIRフィルタ次数: {numtaps - 1}")
print(f"線形位相: True (対称係数)")

8. 制御システムの基礎

8.1 開ループ vs 閉ループ制御

開ループ(Open-Loop)制御: 出力が制御動作に影響を与えません。シンプルですが外乱やモデル誤差に弱いです。

閉ループ(Closed-Loop)制御: 出力を計測して基準入力と比較し、誤差を最小化します。フィードバックにより堅牢性が高まります。

基本的なフィードバックループの伝達関数:

$\frac{Y(s)}{R(s)} = \frac{G(s)}{1 + G(s)H(s)}$

ここで $G(s)$ はプラント、$H(s)$ はセンサの伝達関数です。

8.2 安定性: ラウス・フルビッツ基準

特性方程式 $1 + G(s)H(s) = 0$ のすべての根がs平面の左半面にあれば安定です。

ラウス配列(Routh Array)を用いると、係数だけから極の位置を判定できます。第1列の符号変化の回数が右半面の極の数と等しくなります。

8.3 ボード線図 (Bode Plot)

周波数応答 $H(j\omega)$ を対数スケールで表現します。

• ゲイン線図: $20\log_{10}|H(j\omega)|$ [dB] vs $\log_{10}\omega$
• 位相線図: $\angle H(j\omega)$ [度] vs $\log_{10}\omega$

安定余裕:

• 位相余裕(Phase Margin, PM): ゲインが0 dBのときの位相が -180度からどれだけ余裕があるか
• ゲイン余裕(Gain Margin, GM): 位相が -180度のときのゲインが 0 dB からどれだけ余裕があるか

PM > 45度、GM > 6 dB が一般的に良好な設計の目安です。

from scipy import signal
import matplotlib.pyplot as plt
import numpy as np

# 2次伝達関数: H(s) = 10 / (s^2 + 2s + 10)
num = [10]
den = [1, 2, 10]
sys = signal.TransferFunction(num, den)

w, mag, phase = signal.bode(sys)

fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6))
ax1.semilogx(w, mag)
ax1.set_ylabel('Magnitude (dB)')
ax1.set_title('Bode Plot')
ax1.grid(True, which='both')

ax2.semilogx(w, phase)
ax2.set_ylabel('Phase (degrees)')
ax2.set_xlabel('Frequency (rad/s)')
ax2.grid(True, which='both')

plt.tight_layout()
plt.show()

9. PID制御器

9.1 PID制御動作

PID制御器は3つの制御動作の和です。

$u(t) = K_p e(t) + K_i \int_0^t e(\tau) d\tau + K_d \frac{de(t)}{dt}$

ここで $e(t) = r(t) - y(t)$ は誤差信号です。

• 比例(P)制御: 現在の誤差に比例した制御出力。速い応答が得られますが定常偏差が残ります。
• 積分(I)制御: 累積誤差を解消して定常偏差をゼロに。応答が遅くなります。
• 微分(D)制御: 誤差の変化率に反応してオーバーシュートを抑制。ノイズに敏感。

伝達関数:

$C(s) = K_p + \frac{K_i}{s} + K_d s = K_p\left(1 + \frac{1}{T_i s} + T_d s\right)$

9.2 ジーグラー・ニコルス調整法

方法1(開ループステップ応答):

プロセス反応曲線から無駄時間 $L$ と時定数 $T$ を求めます。

方法2(閉ループ限界感度法):

$K_i = K_d = 0$ として $K_p$ を増加させ、持続振動が起きる限界ゲイン $K_u$ と振動周期 $T_u$ を求めます。

9.3 PythonによるPIDシミュレーション

import numpy as np
import matplotlib.pyplot as plt

class PIDController:
    def __init__(self, kp, ki, kd, dt):
        self.kp = kp
        self.ki = ki
        self.kd = kd
        self.dt = dt
        self.integral = 0.0
        self.prev_error = 0.0

    def update(self, setpoint, measured):
        error = setpoint - measured
        self.integral += error * self.dt
        derivative = (error - self.prev_error) / self.dt
        output = self.kp * error + self.ki * self.integral + self.kd * derivative
        self.prev_error = error
        return output

# 1次プラント: dy/dt = -y + u (時定数1秒)
def plant_step(y, u, dt):
    return y + dt * (-y + u)

dt = 0.01
t = np.arange(0, 10, dt)
setpoint = 1.0

pid = PIDController(kp=2.0, ki=1.0, kd=0.5, dt=dt)
y = 0.0
y_list = []

for _ in t:
    u = pid.update(setpoint, y)
    u = np.clip(u, -10, 10)   # アクチュエータ飽和
    y = plant_step(y, u, dt)
    y_list.append(y)

plt.figure(figsize=(10, 4))
plt.plot(t, y_list, label='Plant Output')
plt.axhline(setpoint, color='r', linestyle='--', label='Setpoint')
plt.xlabel('Time (s)')
plt.ylabel('Output')
plt.title('PID Control Simulation')
plt.legend()
plt.grid(True)
plt.show()

10. 状態空間表現 (State Space)

10.1 状態方程式

連続時間LTIシステムの状態空間表現:

$\dot{\mathbf{x}}(t) = \mathbf{A}\mathbf{x}(t) + \mathbf{B}\mathbf{u}(t)$ $\mathbf{y}(t) = \mathbf{C}\mathbf{x}(t) + \mathbf{D}\mathbf{u}(t)$

• $\mathbf{x}$: 状態ベクトル ($n \times 1$)
• $\mathbf{u}$: 入力ベクトル ($m \times 1$)
• $\mathbf{y}$: 出力ベクトル ($p \times 1$)
• $\mathbf{A}$: システム行列 ($n \times n$)
• $\mathbf{B}$: 入力行列 ($n \times m$)
• $\mathbf{C}$: 出力行列 ($p \times n$)
• $\mathbf{D}$: 直達行列 ($p \times m$)

伝達関数との関係: $H(s) = \mathbf{C}(s\mathbf{I} - \mathbf{A})^{-1}\mathbf{B} + \mathbf{D}$

10.2 可制御性と可観測性

可制御性(Controllability): 任意の初期状態から有限時間内に任意の目標状態へ移行できるか？

可制御性行列 $\mathcal{C} = [\mathbf{B} \; \mathbf{AB} \; \mathbf{A}^2\mathbf{B} \; \cdots \; \mathbf{A}^{n-1}\mathbf{B}]$ のランクが $n$ なら完全可制御。

可観測性(Observability): 出力だけから初期状態を決定できるか？

可観測性行列 $\mathcal{O} = [\mathbf{C}^T \; (\mathbf{CA})^T \; \cdots \; (\mathbf{CA}^{n-1})^T]^T$ のランクが $n$ なら完全可観測。

10.3 線形二次レギュレータ (LQR)

二次コスト関数を最小化する最適状態フィードバック制御:

$J = \int_0^{\infty} (\mathbf{x}^T \mathbf{Q} \mathbf{x} + \mathbf{u}^T \mathbf{R} \mathbf{u}) \, dt$

最適制御入力: $\mathbf{u}^* = -\mathbf{K}\mathbf{x}$

$\mathbf{Q}$ は状態誤差の重み、$\mathbf{R}$ は制御入力エネルギーの重みです。

import numpy as np
from scipy import linalg

# 二重積分器: x_ddot = u
A = np.array([[0, 1], [0, 0]])
B = np.array([[0], [1]])

# LQR重み: 位置誤差を重視
Q = np.diag([10, 1])
R = np.array([[1]])

# 連続時間代数リカッチ方程式を解く
P = linalg.solve_continuous_are(A, B, Q, R)
K = np.linalg.inv(R) @ B.T @ P

print("LQRゲインK:", K)
print("閉ループ極:", np.linalg.eigvals(A - B @ K))

11. AI/MLと信号処理

11.1 1D CNNによる信号分類

1次元畳み込みニューラルネットワークはECG、地震波、音声などの時系列信号のパターン認識に優れています。

import torch
import torch.nn as nn

class Signal1DCNN(nn.Module):
    def __init__(self, num_classes=5):
        super().__init__()
        self.conv_layers = nn.Sequential(
            nn.Conv1d(1, 32, kernel_size=7, padding=3),
            nn.BatchNorm1d(32),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Conv1d(32, 64, kernel_size=5, padding=2),
            nn.BatchNorm1d(64),
            nn.ReLU(),
            nn.MaxPool1d(2),
            nn.Conv1d(64, 128, kernel_size=3, padding=1),
            nn.ReLU(),
            nn.AdaptiveAvgPool1d(8)
        )
        self.classifier = nn.Sequential(
            nn.Flatten(),
            nn.Linear(128 * 8, 256),
            nn.ReLU(),
            nn.Dropout(0.5),
            nn.Linear(256, num_classes)
        )

    def forward(self, x):
        x = self.conv_layers(x)
        return self.classifier(x)

# 使用例
model = Signal1DCNN(num_classes=5)
# 入力形状: (batch_size, channels=1, sequence_length=1000)
x = torch.randn(8, 1, 1000)
output = model(x)
print("出力形状:", output.shape)  # (8, 5)

11.2 STFTとメルスペクトログラム

音声AIシステムでは信号を時間-周波数表現に変換してからCNNで処理する手法が広く使われています。

import librosa
import numpy as np

# STFT (Short-Time Fourier Transform)
# D = librosa.stft(y)
# S_db = librosa.amplitude_to_db(np.abs(D), ref=np.max)

# メルスペクトログラム
# mel_spec = librosa.feature.melspectrogram(y=y, sr=sr, n_mels=128)
# mel_db   = librosa.power_to_db(mel_spec, ref=np.max)

# MFCC (メル周波数ケプストラム係数)
# mfcc = librosa.feature.mfcc(y=y, sr=sr, n_mfcc=13)

print("STFT: ハン窓, hop_length=512, n_fft=2048")
print("メルフィルタバンク: 128バンド, 0Hz〜ナイキスト")
print("MFCC: 音声認識では13次元がよく使われる")

11.3 信号処理 + ディープラーニングの応用例

• ECG不整脈分類: 1D CNN + LSTMによる心拍リズム異常の検出
• 異常振動検知: FFT特徴抽出後にオートエンコーダで異常検出
• 音声認識: MFCC + Transformerアーキテクチャ(Whisper, wav2vec 2.0)
• 地震検知: 地震計信号のP波・S波の自動分類
• 脳波(EEG)解析: 周波数帯域別パワースペクトルで集中度・睡眠状態を分類

12. 確認問題

参考文献

1. Oppenheim, A. V. & Schafer, R. W. — Discrete-Time Signal Processing (3rd ed.), Prentice Hall, 2010
2. Ogata, K. — Modern Control Engineering (5th ed.), Prentice Hall, 2010
3. Proakis, J. G. & Manolakis, D. G. — Digital Signal Processing, Prentice Hall, 2007
4. SciPy Signal Processing Documentation — https://docs.scipy.org/doc/scipy/reference/signal.html
5. MIT OpenCourseWare 6.003 — Signals and Systems — https://ocw.mit.edu/courses/6-003-signals-and-systems-fall-2011/
6. NumPy FFT ドキュメント — https://numpy.org/doc/stable/reference/routines.fft.html