記事概要
本記事では、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()