Files
tcad-bodeplot/dea/filter_design.py
T
2026-05-14 09:41:00 +08:00

93 lines
4.5 KiB
Python
Executable File

import numpy as np
from fastapi import HTTPException
from scipy import signal
from .schemas import DesignParams
from .validation import normalize_coefficients, require_finite, validate_frequency
def _normalize_dc_gain(b_values, a_values):
dc_gain = np.sum(b_values) / np.sum(a_values)
return b_values / dc_gain
def design_coefficients(params: DesignParams):
fs_val = require_finite(params.fs, "fs")
f_type = params.filter_type
if f_type == "Lowpass (低通)":
validate_frequency("低通截止頻率", params.lp_fc, fs_val)
b_new, a_new = signal.butter(params.lp_order, params.lp_fc, btype="low", fs=fs_val)
elif f_type == "Highpass (高通)":
validate_frequency("高通截止頻率", params.hp_fc, fs_val)
b_new, a_new = signal.butter(params.hp_order, params.hp_fc, btype="high", fs=fs_val)
elif f_type == "Bandpass (帶通)":
f_low, f_high = params.bp_f_low, params.bp_f_high
validate_frequency("帶通下截止頻率", f_low, fs_val)
validate_frequency("帶通上截止頻率", f_high, fs_val)
if f_low >= f_high:
raise HTTPException(status_code=400, detail="帶通下截止頻率必須小於上截止頻率")
b_new, a_new = signal.butter(params.bp_order, [f_low, f_high], btype="bandpass", fs=fs_val)
elif f_type == "Notch (陷波器)":
validate_frequency("陷波中心頻率", params.notch_f0, fs_val)
b_new, a_new = signal.iirnotch(params.notch_f0, params.notch_q, fs=fs_val)
elif f_type == "1P1Z (一極一零)":
validate_frequency("1P1Z 零點頻率", params.opoz_fz, fs_val, below_nyquist=False)
validate_frequency("1P1Z 極點頻率", params.opoz_fp, fs_val, below_nyquist=False)
b_s = [1.0, 2 * np.pi * params.opoz_fz]
a_s = [1.0, 2 * np.pi * params.opoz_fp]
b_new, a_new = signal.bilinear(b_s, a_s, fs=fs_val)
b_new = _normalize_dc_gain(b_new, a_new)
elif f_type == "2P1Z (二極一零)":
validate_frequency("2P1Z 零點頻率", params.tp1z_fz, fs_val, below_nyquist=False)
validate_frequency("2P1Z 極點頻率 1", params.tp1z_fp1, fs_val, below_nyquist=False)
validate_frequency("2P1Z 極點頻率 2", params.tp1z_fp2, fs_val, below_nyquist=False)
w_z = 2 * np.pi * params.tp1z_fz
w_p1, w_p2 = 2 * np.pi * params.tp1z_fp1, 2 * np.pi * params.tp1z_fp2
b_s = [1.0, w_z]
a_s = [1.0, w_p1 + w_p2, w_p1 * w_p2]
b_new, a_new = signal.bilinear(b_s, a_s, fs=fs_val)
b_new = _normalize_dc_gain(b_new, a_new)
elif f_type == "2P2Z (二極二零)":
validate_frequency("2P2Z 極點頻率 1", params.tptz_fp1, fs_val, below_nyquist=False)
validate_frequency("2P2Z 極點頻率 2", params.tptz_fp2, fs_val, below_nyquist=False)
validate_frequency("2P2Z 零點頻率 1", params.tptz_fz1, fs_val, below_nyquist=False)
validate_frequency("2P2Z 零點頻率 2", params.tptz_fz2, fs_val, below_nyquist=False)
w_p1, w_p2 = 2 * np.pi * params.tptz_fp1, 2 * np.pi * params.tptz_fp2
w_z1, w_z2 = 2 * np.pi * params.tptz_fz1, 2 * np.pi * params.tptz_fz2
b_s = [1.0, w_z1 + w_z2, w_z1 * w_z2]
a_s = [1.0, w_p1 + w_p2, w_p1 * w_p2]
b_new, a_new = signal.bilinear(b_s, a_s, fs=fs_val)
b_new = _normalize_dc_gain(b_new, a_new)
elif f_type == "PID 控制器":
ki_term = params.pid_ki / (2.0 * fs_val)
kd_term = params.pid_kd * fs_val
b_new = [
params.pid_kp + ki_term + kd_term,
-params.pid_kp + ki_term - 2.0 * kd_term,
kd_term,
]
a_new = [1.0, -1.0, 0.0]
elif f_type == "SOGI-Alpha (帶通)":
validate_frequency("SOGI 中心頻率", params.sogi_f0, fs_val)
w0 = 2 * np.pi * params.sogi_f0
b_s = [params.sogi_k * w0, 0.0]
a_s = [1.0, params.sogi_k * w0, w0**2]
b_new, a_new = signal.bilinear(b_s, a_s, fs=fs_val)
elif f_type == "SOGI-Beta (低通)":
validate_frequency("SOGI 中心頻率", params.sogi_f0, fs_val)
w0 = 2 * np.pi * params.sogi_f0
b_s = [params.sogi_k * (w0**2)]
a_s = [1.0, params.sogi_k * w0, w0**2]
b_new, a_new = signal.bilinear(b_s, a_s, fs=fs_val)
else:
raise HTTPException(status_code=400, detail="未知的濾波器類型")
return normalize_coefficients(b_new, a_new)
def design_response(params: DesignParams):
b_arr, a_arr = design_coefficients(params)
return {"b": b_arr.tolist(), "a": a_arr.tolist()}