サンプリング・z⁻¹遅延・FFT・z平面・デジタル微分と積分

# =========================================
# Demo: Sampling, z^-1 Delay, FFT, z-plane,
#       Digital Differentiation and Integration
# =========================================
import numpy as np
import matplotlib.pyplot as plt

# -------------------------------
# Parameter Settings
# -------------------------------
fs = 20          # Sampling frequency [Hz]
f = 1            # Signal frequency [Hz]
N = 40           # Number of samples
T = 1 / fs       # Sampling period

n = np.arange(N)                      # Discrete-time index
x = np.sin(2 * np.pi * f * n / fs)    # Input sine wave

# z^-1: One-sample delay
x_delay = np.roll(x, 1)
x_delay[0] = 0   # replace first element with zero

# -------------------------------
# Difference Operator (1 - z^-1)
# -------------------------------
x_diff = x - x_delay   # digital differentiation

# -------------------------------
# Integration Operator (1 / (1 - z^-1))
# -------------------------------
x_int = np.zeros_like(x)
for k in range(1, N):
    x_int[k] = x_int[k-1] + x[k]   # cumulative sum

# -------------------------------
# FFT Spectrum
# -------------------------------
X_fft = np.fft.fft(x, 512)
freqs = np.fft.fftfreq(len(X_fft), d=1/fs)

# -------------------------------
# z-plane points (unit circle)
# -------------------------------
z_points = [np.exp(-1j * 2 * np.pi * k / N) for k in range(N)]

# -------------------------------
# Plotting
# -------------------------------
fig, axs = plt.subplots(4, 1, figsize=(8, 12))

# (1) Time-domain: sine wave and delay
axs[0].stem(n, x, basefmt=" ", linefmt="C0-", markerfmt="C0o", label="x[n]")
axs[0].stem(n, x_delay, basefmt=" ", linefmt="C1--", markerfmt="C1s", label="z^-1 x[n]")
axs[0].set_title("Sampled Sine Wave and z^-1 Delay")
axs[0].set_xlabel("n (samples)")
axs[0].set_ylabel("Amplitude")
axs[0].legend()
axs[0].grid(True)

# (2) Difference and Integration
axs[1].plot(n, x, "k-", label="x[n] (input)")
axs[1].plot(n, x_diff, "r-", label="(1 - z^-1)x[n] (difference)")
axs[1].plot(n, x_int/np.max(np.abs(x_int)), "b-", label="Integrated x[n] (scaled)")
axs[1].set_title("Difference and Integration Operators")
axs[1].set_xlabel("n (samples)")
axs[1].set_ylabel("Amplitude")
axs[1].legend()
axs[1].grid(True)

# (3) FFT spectrum
axs[2].plot(freqs[:len(freqs)//2], np.abs(X_fft[:len(X_fft)//2]))
axs[2].set_title("FFT Spectrum of Input x[n]")
axs[2].set_xlabel("Frequency [Hz]")
axs[2].set_ylabel("|X(f)|")
axs[2].grid(True)

# (4) z-plane
theta = np.linspace(0, 2*np.pi, 200)
axs[3].plot(np.cos(theta), np.sin(theta), "k--", label="Unit Circle")
axs[3].scatter([z.real for z in z_points], [z.imag for z in z_points], color="red")
for k, z in enumerate(z_points):
    axs[3].text(z.real*1.1, z.imag*1.1, f"z^(-{k})", fontsize=8)
axs[3].set_title("z-plane: z^-k Locations")
axs[3].set_xlabel("Re")
axs[3].set_ylabel("Im")
axs[3].grid(True)
axs[3].set_aspect("equal", adjustable="box")
axs[3].legend()

plt.tight_layout()
plt.show()