ｚ変換とデジタルフィルタ

1. 定義

Z変換は、離散信号 $x[n]$ を複素数平面に写像する方法です。

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

ここでの $z$ は複素数変数です。
よく次のように置き換えられます：

$$
z = r e^{j\omega}
$$

• $r$：複素平面での半径（大きさ）
• $\omega$：角周波数（単位ラジアン/サンプル）

---

2. 遅延素子と $z^{-1}$

• 1サンプル遅延 を表す演算子が $z^{-1}$ です。

• つまり：

  $$
  X(z) \cdot z^{-1} \quad \Longleftrightarrow \quad x[n-1]
  $$

例：
入力 $x[n]$ を1サンプル遅らせた信号は、Z領域で $X(z)z^{-1}$ になります。

---

3. 単位円上の周波数応答

Z変換の変数を $z = e^{j\omega}$ と置くと、周波数応答になります。

• $|z| = 1$（単位円）上に制限すると → 実際のフィルタの振幅・位相特性が見られる。
• アナログの $s = j\omega$ に対応するのが、デジタルの $z = e^{j\omega}$。

---

4. デジタル微分器・積分器との関係

■ 微分器

• 前進差分:

  $$
  H(z) = 1 - z^{-1}
  $$

  周波数領域で

  $$
  H(e^{j\omega}) = 1 - e^{-j\omega}
  $$

  → $j\omega$ に比例する特性（高域で強調）

• 中心差分:

  $$
  H(z) = \dfrac{z - z^{-1}}{2}
  $$

  周波数領域で

  $$
  H(e^{j\omega}) = j \sin(\omega)
  $$

  → より対称的で精度が高い

---

■ 積分器

• 単純積分器:

  $$
  H(z) = \dfrac{1}{1 - z^{-1}}
  $$

  周波数領域で

  $$
  H(e^{j\omega}) = \dfrac{1}{1 - e^{-j\omega}}
  $$

  → $1/(j\omega)$ に似た特性（低域で強調）

• 遅延付き積分器:

  $$
  H(z) = \dfrac{z^{-1}}{1 - z^{-1}}
  $$

  → 実装が安定で、因果性が明確。

---

import numpy as np
import matplotlib.pyplot as plt

# ================= Parameters / パラメータ =================
fs = 100       # Sampling frequency [Hz] / サンプリング周波数
Ts = 1/fs      # Sampling period [s] / サンプリング周期
f_in = 5       # Input sine frequency [Hz] / 入力サイン波の周波数
N = 200        # Number of samples / サンプル数
t = np.arange(N) * Ts  # Time vector / 時間ベクトル

# 入力サイン波 Input sine wave
x = np.sin(2*np.pi*f_in*t)

# ================= Differentiators / 微分器 =================
# 前進差分 Forward difference
y_fwd = np.zeros(N)
y_fwd[1:] = x[1:] - x[:-1]

# 中心差分 Central difference
y_ctr = np.zeros(N)
y_ctr[1:-1] = (x[2:] - x[:-2]) / 2

# ================= Integrators / 積分器 =================
# 単純積分器 Simple accumulator
y_int = np.zeros(N)
for n in range(1, N):
    y_int[n] = y_int[n-1] + x[n]

# 遅延付き積分器 Integrator with unit delay
y_int_delay = np.zeros(N)
for n in range(1, N):
    y_int_delay[n] = y_int_delay[n-1] + x[n-1]

# ================= Plot time-domain responses / 時間応答 =================
plt.figure(figsize=(12,8))
plt.subplot(3,1,1)
plt.plot(t, x, 'k', label='Input sine')
plt.plot(t, y_fwd, 'r', label='Forward diff')
plt.plot(t, y_ctr, 'b', label='Central diff')
plt.title('Digital Differentiators (Time Domain)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

plt.subplot(3,1,2)
plt.plot(t, x, 'k', label='Input sine')
plt.plot(t, y_int, 'g', label='Simple integrator')
plt.plot(t, y_int_delay, 'm', label='Integrator with delay')
plt.title('Digital Integrators (Time Domain)')
plt.xlabel('Time [s]')
plt.ylabel('Amplitude')
plt.legend()
plt.grid(True)

# ================= Frequency response / 周波数応答 =================
w = np.linspace(0, np.pi, 512)  # normalized frequency / 正規化角周波数
H_fwd = 1 - np.exp(-1j*w)
H_ctr = (np.exp(1j*w) - np.exp(-1j*w)) / 2
H_int = 1 / (1 - np.exp(-1j*w))
H_int_delay = np.exp(-1j*w) / (1 - np.exp(-1j*w))

plt.subplot(3,1,3)
plt.plot(w/np.pi, 20*np.log10(np.abs(H_fwd)), 'r', label='Forward diff')
plt.plot(w/np.pi, 20*np.log10(np.abs(H_ctr)), 'b', label='Central diff')
plt.plot(w/np.pi, 20*np.log10(np.abs(H_int)), 'g', label='Simple int')
plt.plot(w/np.pi, 20*np.log10(np.abs(H_int_delay)), 'm', label='Int with delay')
plt.title('Frequency Responses |H(e^{jω})|')
plt.xlabel('Normalized frequency (×π rad/sample)')
plt.ylabel('Magnitude [dB]')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()