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

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}")

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()

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()