ハイパスフィルタの設計とPython実装

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

# --- パラメータ設定 ---
fs = 1000.0      # サンプリング周波数 [Hz]
fc = 100.0       # 遮断周波数 [Hz]
order = 4        # フィルタ次数

# --- バターワース ハイパスフィルタの設計 ---
nyq = fs / 2.0   # ナイキスト周波数
wn = fc / nyq    # 正規化遮断周波数

b, a = signal.butter(order, wn, btype='high')

# --- 周波数応答の計算 ---
w, h = signal.freqz(b, a, worN=8000, fs=fs)

# --- プロット ---
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))

# ゲイン特性 (dB)
ax1.semilogx(w, 20 * np.log10(np.abs(h) + 1e-10))
ax1.axvline(fc, color='r', linestyle='--', label=f'$f_c$ = {fc} Hz')
ax1.axhline(-3, color='gray', linestyle=':', label='-3 dB')
ax1.set_xlabel('Frequency [Hz]')
ax1.set_ylabel('Gain [dB]')
ax1.set_title(f'Butterworth Highpass Filter (order={order})')
ax1.legend()
ax1.grid(True, which='both', alpha=0.3)
ax1.set_xlim([1, nyq])
ax1.set_ylim([-80, 5])

# 位相特性
ax2.semilogx(w, np.angle(h, deg=True))
ax2.axvline(fc, color='r', linestyle='--')
ax2.set_xlabel('Frequency [Hz]')
ax2.set_ylabel('Phase [degrees]')
ax2.set_title('Phase Response')
ax2.grid(True, which='both', alpha=0.3)
ax2.set_xlim([1, nyq])

plt.tight_layout()
plt.savefig('highpass_response.png', dpi=150, bbox_inches='tight')
plt.show()

from scipy.signal import firwin, freqz

# --- FIR ハイパスフィルタ ---
numtaps = 101    # フィルタ係数の数（奇数推奨）
fir_hpf = firwin(
    numtaps,
    fc,
    pass_zero=False,  # ハイパス指定
    fs=fs,
    window='hamming'
)

w_fir, h_fir = freqz(fir_hpf, worN=8000, fs=fs)

# IIRとFIRの比較プロット
plt.figure(figsize=(10, 5))
plt.semilogx(w, 20 * np.log10(np.abs(h) + 1e-10),
             label='IIR Butterworth (order=4)', linewidth=2)
plt.semilogx(w_fir, 20 * np.log10(np.abs(h_fir) + 1e-10),
             label=f'FIR Hamming (taps={numtaps})', linewidth=2, linestyle='--')
plt.axvline(fc, color='r', linestyle=':', alpha=0.7, label=f'$f_c$ = {fc} Hz')
plt.axhline(-3, color='gray', linestyle=':', alpha=0.7)
plt.xlabel('Frequency [Hz]')
plt.ylabel('Gain [dB]')
plt.title('Highpass Filter Comparison: IIR vs FIR')
plt.legend()
plt.grid(True, which='both', alpha=0.3)
plt.xlim([1, nyq])
plt.ylim([-80, 5])
plt.tight_layout()
plt.show()

# --- テスト信号の生成 ---
t = np.linspace(0, 1.0, int(fs), endpoint=False)  # 1秒分

# 目的信号: 200 Hz のサイン波（高周波成分）
signal_pure = np.sin(2 * np.pi * 200 * t)

# ノイズ: 5 Hz の低周波ドリフト + ホワイトノイズ
noise = (
    2.0 * np.sin(2 * np.pi * 5 * t) +   # 低周波ドリフト（大振幅）
    0.3 * np.random.randn(len(t))         # ホワイトノイズ
)

x_noisy = signal_pure + noise

# --- フィルタ適用 ---
# lfilter: 因果的フィルタ（リアルタイム処理向け、位相遅れあり）
y_lfilter = signal.lfilter(b, a, x_noisy)

# filtfilt: ゼロ位相フィルタ（オフライン処理向け、位相遅れなし）
y_filtfilt = signal.filtfilt(b, a, x_noisy)

# --- 比較プロット ---
fig, axes = plt.subplots(3, 1, figsize=(12, 8), sharex=True)

axes[0].plot(t[:300], x_noisy[:300], alpha=0.8, label='ノイズ混じり信号（低周波ドリフトあり）')
axes[0].set_ylabel('Amplitude')
axes[0].set_title('Input Signal (with low-frequency drift)')
axes[0].legend()
axes[0].grid(True, alpha=0.3)

axes[1].plot(t[:300], y_lfilter[:300], color='orange', label='lfilter (因果的)')
axes[1].plot(t[:300], signal_pure[:300], color='gray', linestyle='--', alpha=0.6, label='真の信号')
axes[1].set_ylabel('Amplitude')
axes[1].set_title('After lfilter (位相遅れあり)')
axes[1].legend()
axes[1].grid(True, alpha=0.3)

axes[2].plot(t[:300], y_filtfilt[:300], color='green', label='filtfilt (ゼロ位相)')
axes[2].plot(t[:300], signal_pure[:300], color='gray', linestyle='--', alpha=0.6, label='真の信号')
axes[2].set_ylabel('Amplitude')
axes[2].set_xlabel('Time [s]')
axes[2].set_title('After filtfilt (ゼロ位相)')
axes[2].legend()
axes[2].grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('highpass_filtering_result.png', dpi=150, bbox_inches='tight')
plt.show()

# 次数の違いによる周波数応答の比較
fig, ax = plt.subplots(figsize=(10, 6))

for ord_n in [1, 2, 4, 8]:
    b_n, a_n = signal.butter(ord_n, wn, btype='high')
    w_n, h_n = signal.freqz(b_n, a_n, worN=8000, fs=fs)
    ax.semilogx(w_n, 20 * np.log10(np.abs(h_n) + 1e-10), label=f'order={ord_n}')

ax.axhline(-3, color='gray', linestyle=':', label='-3 dB')
ax.axvline(fc, color='r', linestyle='--', alpha=0.5, label=f'$f_c$ = {fc} Hz')
ax.set_xlabel('Frequency [Hz]')
ax.set_ylabel('Gain [dB]')
ax.set_title('Effect of Filter Order on Highpass Response')
ax.legend()
ax.grid(True, which='both', alpha=0.3)
ax.set_xlim([1, nyq])
ax.set_ylim([-80, 5])
plt.tight_layout()
plt.show()