記事概要

本記事では、Pythonを使用してさまざまな信号処理技術を適用し、フィルタリング技術、周波数解析、増幅回路の解析を行う方法について詳述します。特に、直列・並列接続の抵抗計算、電圧計算、デシベル換算、ボード線図、FFT（高速フーリエ変換）を利用した信号解析の基本概念を紹介します。

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

# SI接頭辞：記号, 10の指数, カタカナ読み
si_prefixes = [
    ("Y",  24, "ヨタ"),     # yotta
    ("Z",  21, "ゼタ"),     # zetta
    ("E",  18, "エクサ"),   # exa
    ("P",  15, "ペタ"),     # peta
    ("T",  12, "テラ"),     # tera
    ("G",   9, "ギガ"),     # giga
    ("M",   6, "メガ"),     # mega
    ("k",   3, "キロ"),     # kilo
    ("h",   2, "ヘクト"),   # hecto
    ("da",  1, "デカ"),     # deca
    ("",    0, "（なし）"), # no prefix
    ("d",  -1, "デシ"),     # deci
    ("c",  -2, "センチ"),   # centi
    ("m",  -3, "ミリ"),     # milli
    ("µ",  -6, "マイクロ"), # micro
    ("n",  -9, "ナノ"),     # nano
    ("p", -12, "ピコ"),     # pico
    ("f", -15, "フェムト"), # femto
    ("a", -18, "アト"),     # atto
    ("z", -21, "ゼプト"),   # zepto
    ("y", -24, "ヨクト")    # yocto
]

# ヘッダー表示
print("SI接頭辞 一覧")
print("{:<4} {:>6} {:<10}".format("記号", "10の指数", "カタカナ読み"))
print("-" * 26)

# 各行の出力
for symbol, exponent, katakana in si_prefixes:
    print("{:<4} {:>6} {:<10}".format(symbol, f"10^{exponent}", katakana))
# --- 基本設定 ---
fs = 1000  # サンプリング周波数 (Hz)
t = np.arange(0, 1, 1/fs)  # 時間軸 (1秒間)
f1 = 50  # サイン波1の周波数 (Hz)
f2 = 150  # サイン波2の周波数 (Hz)


# --- 直列接続の抵抗計算 ---
def series_resistance(R1, R2):
    return R1 + R2

# --- 並列接続の抵抗計算 ---
def parallel_resistance(R1, R2):
    return (R1 * R2) / (R1 + R2)

R1, R2 = 100, 200  # 抵抗値 (Ω)
print(f"直列接続の抵抗: {series_resistance(R1, R2)} Ω")
print(f"並列接続の抵抗: {parallel_resistance(R1, R2)} Ω")

# --- 電圧の計算 (オームの法則) ---
def calculate_voltage(I, R):
    return I * R

I = 0.01  # 電流 (A)
R = 100  # 抵抗 (Ω)
V = calculate_voltage(I, R)
print(f"電圧: {V} V")

# --- デシベルの計算 (増幅) ---
def calculate_db(amplification_factor):
    return 20 * np.log10(amplification_factor)

# 増幅倍数からデシベル計算
amplification_factor = 2  # 増幅倍数
dB = calculate_db(amplification_factor)
print(f"増幅 {amplification_factor}倍に対するデシベル: {dB} dB")

# デシベルから増幅倍数に変換
def db_to_amplification(db):
    return 10 ** (db / 20)

amplification_from_db = db_to_amplification(dB)
print(f"デシベル {dB} からの増幅倍数: {amplification_from_db}倍")


# --- サンプル信号の生成 ---
# サイン波の生成
sin_wave = np.sin(2 * np.pi * f1 * t) + np.sin(2 * np.pi * f2 * t)
# 方形波の生成
square_wave = np.sign(np.sin(2 * np.pi * f1 * t))
# 三角波の生成
triangular_wave = 2 * (t * f1 - np.floor(t * f1 + 0.5))

# --- 増幅回路の伝達関数 ---
def transfer_function(k, tau, s):
    """ 増幅回路と一次遅れ系の伝達関数 """
    return k / (tau * s + 1)

k = 1  # ゲイン
tau = 0.1  # 時定数
s = 1j  # 複素数
H = transfer_function(k, tau, s)
print(f"一次遅れ系伝達関数の値: {H}")

# --- ボード線図（Bode Plot）の作成 ---
# 増幅回路の伝達関数を使用してボード線図を描く
# W = 2πf として周波数を定義
frequencies = np.logspace(0, 3, 100)  # 10^0 to 10^3 Hzの対数間隔
s = 1j * 2 * np.pi * frequencies  # jω
H_w = transfer_function(k, tau, s)

# ゲインのプロット
gain = 20 * np.log10(np.abs(H_w))  # ゲインをdBで計算

# 位相のプロット
phase = np.angle(H_w, deg=True)  # 位相を度で計算

# ボード線図を描画
fig, axs = plt.subplots(2, 1, figsize=(10, 8))

# ゲインプロット
axs[0].semilogx(frequencies, gain)
axs[0].set_title("Bode Plot: Gain")
axs[0].set_xlabel("Frequency (Hz)")
axs[0].set_ylabel("Gain (dB)")
axs[0].grid(True)

# 位相プロット
axs[1].semilogx(frequencies, phase)
axs[1].set_title("Bode Plot: Phase")
axs[1].set_xlabel("Frequency (Hz)")
axs[1].set_ylabel("Phase (Degrees)")
axs[1].grid(True)

plt.tight_layout()
plt.show()

# --- FFT（高速フーリエ変換）を使用した周波数解析 ---
def fft_analysis(signal):
    """ FFT解析 """
    n = len(signal)
    frequencies = np.fft.fftfreq(n, 1/fs)
    fft_values = np.fft.fft(signal)
    return frequencies, fft_values

frequencies, fft_values = fft_analysis(sin_wave)

# --- プロット設定 ---
# 1. 時間領域の信号波形
fig, axs = plt.subplots(2, 1, figsize=(10, 8))

# サイン波、方形波、三角波をプロット
axs[0].plot(t, sin_wave, label="Sin Wave")
axs[0].plot(t, square_wave, label="Square Wave", alpha=0.6)
axs[0].plot(t, triangular_wave, label="Triangle Wave", alpha=0.6)
axs[0].set_title("Time Domain Signals")
axs[0].set_xlabel("Time (s)")
axs[0].set_ylabel("Amplitude")
axs[0].legend()
axs[0].grid(True)

# フィルタ後の信号をプロット
axs[1].plot(t, sin_wave, label="Original Signal")
axs[1].set_title("Original Signal (No Filter)")
axs[1].set_xlabel("Time (s)")
axs[1].set_ylabel("Amplitude")
axs[1].legend()
axs[1].grid(True)

plt.tight_layout()
plt.show()

# --- 周波数領域の解析（FFT） ---
plt.figure(figsize=(10, 6))
plt.plot(frequencies[:len(sin_wave)//2], np.abs(fft_values)[:len(sin_wave)//2])  # 正の周波数成分のみ表示
plt.title("FFT of the Signal")
plt.xlabel("Frequency (Hz)")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()

# --- 逆ラプラス変換の実施 ---
# 逆ラプラス変換（ここでは一次遅れ系のインパルス応答を逆ラプラス変換）
# scipy.signalには直接的な逆ラプラス変換の関数はないので、インパルス応答を利用

time_response = signal.lti([k], [tau, 1])  # 伝達関数 H(s) = k / (tau * s + 1)
t, y = signal.step(time_response)  # インパルス応答の取得

plt.figure(figsize=(10, 6))
plt.plot(t, y)
plt.title("Impulse Response of First-Order System")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.grid(True)
plt.show()

# --- 畳み込み積分によるフィルタの適用 ---
def convolution_integral(signal, kernel):
    """ 畳み込み積分による信号フィルタリング """
    return np.convolve(signal, kernel, mode='same')

kernel = np.ones(50) / 50  # 移動平均フィルタ
filtered_signal = convolution_integral(sin_wave, kernel)

# フィルタリングされた信号をプロット
plt.figure(figsize=(10, 6))
plt.plot(t, filtered_signal, label="Filtered Signal (Convolution)")
plt.title("Signal after Convolution Integral")
plt.xlabel("Time (s)")
plt.ylabel("Amplitude")
plt.grid(True)
plt.legend()
plt.show()