ノッチフィルタの設計とPython実装：電源ノイズ除去の実践

import numpy as np
from scipy import signal

fs = 1000
f0 = 50

for Q_set in [5, 15, 30, 60, 100]:
    b, a = signal.iirnotch(f0, Q_set, fs)
    _, poles, _ = signal.tf2zpk(b, a)
    r = np.abs(poles[0])

    # 式(3)による近似Q
    omega0 = 2 * np.pi * f0 / fs
    Q_formula = omega0 / (2 * (1 - r))

    # freqzで-3dB帯域幅を実測してQを逆算
    w, h = signal.freqz(b, a, worN=400_000, fs=fs)
    mag_db = 20 * np.log10(np.maximum(np.abs(h), 1e-16))
    li = np.where((mag_db[:-1] > -3) & (mag_db[1:] <= -3) & (w[:-1] < f0))[0]
    ri = np.where((mag_db[:-1] <= -3) & (mag_db[1:] > -3) & (w[:-1] > f0))[0]
    bw_measured = w[ri[0]] - w[li[-1]]
    Q_measured = f0 / bw_measured

    err_pct = (Q_formula - Q_measured) / Q_measured * 100
    print(f"Q_set={Q_set:4d}  r={r:.6f}  Q_formula={Q_formula:.4f}  "
          f"Q_measured={Q_measured:.4f}  誤差={err_pct:+.3f}%")

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

# --- パラメータ ---
fs = 1000       # サンプリング周波数 [Hz]
f0 = 50         # ノッチ周波数 [Hz]
Q = 30          # 品質係数

# --- ノッチフィルタの設計 ---
b, a = signal.iirnotch(f0, Q, fs)

# --- 零点・極の計算 ---
zeros, poles, _ = signal.tf2zpk(b, a)

# --- 零点・極配置図のプロット ---
fig, ax = plt.subplots(figsize=(6, 6))

# 単位円を描画
theta = np.linspace(0, 2 * np.pi, 200)
ax.plot(np.cos(theta), np.sin(theta), 'k-', linewidth=0.5)

# 零点と極をプロット
ax.scatter(zeros.real, zeros.imag, marker='o', s=100,
           facecolors='none', edgecolors='blue', linewidths=2, label='Zeros')
ax.scatter(poles.real, poles.imag, marker='x', s=100,
           color='red', linewidths=2, label='Poles')

ax.set_xlabel('Real')
ax.set_ylabel('Imaginary')
ax.set_title(f'Pole-Zero Plot (f0={f0} Hz, Q={Q})')
ax.set_aspect('equal')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

# 極の半径を表示
print(f"極の半径 r = {np.abs(poles[0]):.6f}")
print(f"零点の半径 = {np.abs(zeros[0]):.6f}")

極の半径 r = 0.994778
零点の半径 = 1.000000

from scipy.signal import iirnotch

# パラメータ
f0 = 50   # ノッチ周波数 [Hz]
Q = 30    # 品質係数
fs = 1000 # サンプリング周波数 [Hz]

# フィルタ係数の計算
b, a = iirnotch(f0, Q, fs)

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

fs = 1000
f0 = 50

fig, ax = plt.subplots(figsize=(10, 6))

for Q in [5, 15, 30, 60]:
    b, a = signal.iirnotch(f0, Q, fs)
    w, h = signal.freqz(b, a, worN=4096, fs=fs)
    ax.plot(w, 20 * np.log10(np.maximum(np.abs(h), 1e-12)),
            label=f'Q = {Q}')

ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Notch Filter Frequency Response (f0 = 50 Hz)')
ax.set_xlim(0, 200)
ax.set_ylim(-60, 5)
ax.axvline(f0, color='gray', linestyle='--', alpha=0.5, label=f'f0 = {f0} Hz')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

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

# --- テスト信号の生成 ---
np.random.seed(42)
fs = 1000           # サンプリング周波数 [Hz]
T = 2.0             # 信号長 [s]
t = np.arange(0, T, 1/fs)
N = len(t)

# 有用信号: 10Hzと30Hzの正弦波
useful_signal = np.sin(2 * np.pi * 10 * t) + 0.5 * np.sin(2 * np.pi * 30 * t)

# 電源ノイズ: 50Hz基本波 + 100Hz, 150Hzの高調波
power_noise = (0.8 * np.sin(2 * np.pi * 50 * t)
               + 0.4 * np.sin(2 * np.pi * 100 * t)
               + 0.2 * np.sin(2 * np.pi * 150 * t))

# ランダムノイズ
random_noise = 0.3 * np.random.randn(N)

# 観測信号
observed = useful_signal + power_noise + random_noise

# --- カスケード接続のノッチフィルタ設計 ---
Q = 30
notch_freqs = [50, 100, 150]  # 基本波 + 高調波

# 各ノッチ周波数のフィルタ係数をSOS形式で結合
sos_list = []
for f_notch in notch_freqs:
    b, a = signal.iirnotch(f_notch, Q, fs)
    sos = signal.tf2sos(b, a)
    sos_list.append(sos[0])

sos_cascade = np.array(sos_list)

# --- フィルタ適用（ゼロ位相フィルタリング） ---
filtered = signal.sosfiltfilt(sos_cascade, observed)

# --- 周波数スペクトルの計算 ---
freqs = np.fft.rfftfreq(N, 1/fs)
spectrum_before = np.abs(np.fft.rfft(observed)) / N
spectrum_after = np.abs(np.fft.rfft(filtered)) / N

# --- プロット ---
fig, axes = plt.subplots(3, 1, figsize=(10, 10))

# (a) 時間領域: フィルタ適用前
axes[0].plot(t, observed, alpha=0.7, label='Observed')
axes[0].plot(t, useful_signal, 'k--', alpha=0.5, label='True signal')
axes[0].set_xlabel('Time [s]')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Before Notch Filtering')
axes[0].set_xlim(0, 0.5)
axes[0].legend()
axes[0].grid(True, alpha=0.3)

# (b) 時間領域: フィルタ適用後
axes[1].plot(t, filtered, alpha=0.7, label='Filtered')
axes[1].plot(t, useful_signal, 'k--', alpha=0.5, label='True signal')
axes[1].set_xlabel('Time [s]')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('After Notch Filtering (50, 100, 150 Hz)')
axes[1].set_xlim(0, 0.5)
axes[1].legend()
axes[1].grid(True, alpha=0.3)

# (c) 周波数領域: 前後比較
axes[2].plot(freqs, 20 * np.log10(np.maximum(spectrum_before, 1e-12)),
             alpha=0.7, label='Before')
axes[2].plot(freqs, 20 * np.log10(np.maximum(spectrum_after, 1e-12)),
             alpha=0.7, label='After')
for f_notch in notch_freqs:
    axes[2].axvline(f_notch, color='red', linestyle='--', alpha=0.3)
axes[2].set_xlabel('Frequency [Hz]')
axes[2].set_ylabel('Magnitude [dB]')
axes[2].set_title('Frequency Spectrum: Before vs After')
axes[2].set_xlim(0, 200)
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# --- 定量評価 ---
rmse_before = np.sqrt(np.mean((observed - useful_signal) ** 2))
rmse_after = np.sqrt(np.mean((filtered - useful_signal) ** 2))
snr_before = 10 * np.log10(np.mean(useful_signal**2) / np.mean((observed - useful_signal) ** 2))
snr_after = 10 * np.log10(np.mean(useful_signal**2) / np.mean((filtered - useful_signal) ** 2))
print(f"RMSE: before={rmse_before:.5f}, after={rmse_after:.5f}")
print(f"SNR : before={snr_before:.3f} dB, after={snr_after:.3f} dB")

for f_notch in notch_freqs:
    idx = np.argmin(np.abs(freqs - f_notch))
    atten_db = 20 * np.log10(spectrum_after[idx] / spectrum_before[idx])
    print(f"  {f_notch} Hz 減衰量 = {atten_db:.2f} dB")

for f_useful in [10, 30]:
    idx = np.argmin(np.abs(freqs - f_useful))
    change_db = 20 * np.log10(spectrum_after[idx] / spectrum_before[idx])
    print(f"  有用信号 {f_useful} Hz の変化 = {change_db:.3f} dB")

RMSE: before=0.70573, after=0.29683
SNR : before=0.986 dB, after=8.509 dB
  50 Hz 減衰量 = -26.67 dB
  100 Hz 減衰量 = -30.49 dB
  150 Hz 減衰量 = -33.32 dB
  有用信号 10 Hz の変化 = -0.001 dB
  有用信号 30 Hz の変化 = -0.011 dB

import numpy as np
from scipy import signal

fs = 1000
f0_design = 50.0
drift_range = np.linspace(-1.0, 1.0, 201)  # 設計周波数からのずれ [Hz]

for Q in [10, 30, 60, 100]:
    b, a = signal.iirnotch(f0_design, Q, fs)
    for d in [0.1, 0.2, 0.5, 1.0]:
        idx = np.argmin(np.abs(drift_range - d))
        w, h = signal.freqz(b, a, worN=[2 * np.pi * (f0_design + d) / fs])
        atten_db = 20 * np.log10(np.abs(h[0]))
        print(f"Q={Q:3d}, drift={d:+.1f}Hz -> 減衰量 = {atten_db:.2f} dB")

import numpy as np
from scipy import signal

fs = 1000
f0 = 50

def quantize(coeffs, bits):
    """[-2, 2) の範囲を bits ビットの固定小数点に丸める簡易モデル"""
    scale = 2 ** (bits - 2)
    return np.round(coeffs * scale) / scale

for Q in [30, 100, 300]:
    b, a = signal.iirnotch(f0, Q, fs)
    r_orig = np.max(np.abs(np.roots(a)))
    print(f"Q={Q}: 元の極半径 r = {r_orig:.6f}")
    for bits in [16, 10, 8]:
        a_q = quantize(a, bits)
        r_q = np.max(np.abs(np.roots(a_q)))
        print(f"  {bits}bit量子化後: r = {r_q:.6f}, 安定 = {r_q < 1.0}")

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

fs = 1000
f0 = 50

# ノッチフィルタ（Q=30）
b_notch, a_notch = signal.iirnotch(f0, Q=30, fs=fs)
w_notch, h_notch = signal.freqz(b_notch, a_notch, worN=4096, fs=fs)

# バンドストップフィルタ（バターワース4次、40-60Hz除去）
sos_bs = signal.butter(4, [40, 60], btype='bandstop', fs=fs, output='sos')
w_bs, h_bs = signal.sosfreqz(sos_bs, worN=4096, fs=fs)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(w_notch, 20 * np.log10(np.maximum(np.abs(h_notch), 1e-12)),
        label='Notch (Q=30)')
ax.plot(w_bs, 20 * np.log10(np.maximum(np.abs(h_bs), 1e-12)),
        label='Bandstop Butterworth (40-60 Hz)')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Magnitude [dB]')
ax.set_title('Notch Filter vs Bandstop Filter')
ax.set_xlim(0, 200)
ax.set_ylim(-60, 5)
ax.axvline(f0, color='gray', linestyle='--', alpha=0.5)
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()