tcad-bodeplot/dea_api.py

import numpy as np
import pandas as pd
from scipy import signal
from fastapi import FastAPI, File, UploadFile, Form, HTTPException, Request
from fastapi.staticfiles import StaticFiles
from pydantic import BaseModel, Field
from typing import List
import io
from fastapi.responses import StreamingResponse, RedirectResponse, JSONResponse
import ipaddress

app = FastAPI(title="Difference Equation Analyzer API")

MAX_COEFFICIENTS = 64
MAX_CSV_BYTES = 32 * 1024 * 1024
MAX_PLOT_POINTS = 5000
BODE_POINTS = 500
BODE_MAX_MULTIPLIER = 3.162
NYQUIST_MARGIN = 0.999999
SECURITY_HEADERS = {
    "Content-Security-Policy": (
        "default-src 'self'; "
        "script-src 'self' https://cdn.plot.ly; "
        "style-src 'self' 'unsafe-inline'; "
        "img-src 'self' data:; "
        "connect-src 'self'; "
        "font-src 'self' data:; "
        "object-src 'none'; "
        "base-uri 'self'; "
        "frame-ancestors 'none'"
    ),
    "Referrer-Policy": "no-referrer",
    "X-Content-Type-Options": "nosniff",
    "X-Frame-Options": "DENY",
}

def add_security_headers(response):
    for key, value in SECURITY_HEADERS.items():
        response.headers[key] = value
    return response

@app.middleware("http")
async def restrict_to_lan(request: Request, call_next):
    client_ip = request.client.host
    if client_ip:
        try:
            ip_obj = ipaddress.ip_address(client_ip)
            if not (ip_obj.is_private or ip_obj.is_loopback):
                return add_security_headers(JSONResponse(status_code=403, content={"detail": "Access Denied: Only LAN connections are allowed."}))
        except ValueError:
            pass
    response = await call_next(request)
    return add_security_headers(response)

# 掛載靜態網頁檔案，將預設首頁導向到 /ui/
app.mount("/ui", StaticFiles(directory="static", html=True), name="static")

@app.get("/")
def redirect_to_ui():
    return RedirectResponse(url="/ui/")

class DesignParams(BaseModel):
    filter_type: str
    fs: float = Field(gt=0)
    lp_fc: float = Field(default=1000.0, gt=0)
    lp_order: int = Field(default=1, ge=1, le=8)
    hp_fc: float = Field(default=1000.0, gt=0)
    hp_order: int = Field(default=1, ge=1, le=8)
    bp_f_low: float = Field(default=500.0, gt=0)
    bp_f_high: float = Field(default=2000.0, gt=0)
    bp_order: int = Field(default=1, ge=1, le=8)
    notch_f0: float = Field(default=120.0, gt=0)
    notch_q: float = Field(default=1.0, gt=0)
    opoz_fz: float = Field(default=15000.0, gt=0)
    opoz_fp: float = Field(default=10.0, gt=0)
    tp1z_fz: float = Field(default=200.0, gt=0)
    tp1z_fp1: float = Field(default=10.0, gt=0)
    tp1z_fp2: float = Field(default=5000.0, gt=0)
    tptz_fz1: float = Field(default=200.0, gt=0)
    tptz_fz2: float = Field(default=25000.0, gt=0)
    tptz_fp1: float = Field(default=10.0, gt=0)
    tptz_fp2: float = Field(default=5000.0, gt=0)
    pid_kp: float = 0.003
    pid_ki: float = 10.0
    pid_kd: float = 0.000016
    sogi_f0: float = Field(default=60.0, gt=0)
    sogi_k: float = Field(default=1.414, gt=0)

class BodeParams(BaseModel):
    b: List[float] = Field(min_length=1, max_length=MAX_COEFFICIENTS)
    a: List[float] = Field(min_length=1, max_length=MAX_COEFFICIENTS)
    fs: float = Field(gt=0)

def require_finite(value, name):
    if not np.isfinite(value):
        raise HTTPException(status_code=400, detail=f"{name} 必須是有限數值")
    return float(value)

def validate_frequency(name, value, fs, *, below_nyquist=True):
    value = require_finite(value, name)
    if value <= 0:
        raise HTTPException(status_code=400, detail=f"{name} 必須大於 0")
    if below_nyquist and value >= (fs / 2) * NYQUIST_MARGIN:
        raise HTTPException(status_code=400, detail=f"{name} 必須小於 Nyquist 頻率 ({fs / 2:g} Hz)")
    return value

def parse_coefficients(raw, name):
    try:
        values = [float(x.strip()) for x in raw.replace(',', ' ').split() if x.strip()]
    except ValueError:
        raise HTTPException(status_code=400, detail=f"{name} 係數格式錯誤")
    return validate_coefficients(values, name)

def validate_coefficients(values, name):
    if len(values) == 0:
        raise HTTPException(status_code=400, detail=f"{name} 係數不可為空")
    if len(values) > MAX_COEFFICIENTS:
        raise HTTPException(status_code=400, detail=f"{name} 係數最多 {MAX_COEFFICIENTS} 個")
    coeffs = np.asarray(values, dtype=float)
    if not np.all(np.isfinite(coeffs)):
        raise HTTPException(status_code=400, detail=f"{name} 係數必須都是有限數值")
    if name == "a" and np.isclose(coeffs[0], 0.0):
        raise HTTPException(status_code=400, detail="a[0] 不可為 0")
    return coeffs

def normalize_coefficients(b_values, a_values):
    b_arr = validate_coefficients(b_values, "b")
    a_arr = validate_coefficients(a_values, "a")
    if not np.isclose(a_arr[0], 1.0):
        b_arr = b_arr / a_arr[0]
        a_arr = a_arr / a_arr[0]
    return b_arr, a_arr

def downsample_for_plot(index, original, filtered):
    total = len(index)
    if total <= MAX_PLOT_POINTS:
        return index, original, filtered, 1
    step = int(np.ceil(total / MAX_PLOT_POINTS))
    return index[::step], original[::step], filtered[::step], step

async def read_csv_upload(file):
    filename = (file.filename or "").lower()
    if filename and not filename.endswith(".csv"):
        raise HTTPException(status_code=400, detail="請上傳 CSV 檔案")
    contents = await file.read(MAX_CSV_BYTES + 1)
    if len(contents) > MAX_CSV_BYTES:
        raise HTTPException(status_code=413, detail=f"CSV 檔案不可超過 {MAX_CSV_BYTES // (1024 * 1024)}MB")
    if not contents.strip():
        raise HTTPException(status_code=400, detail="CSV 不可為空")
    try:
        df = pd.read_csv(io.BytesIO(contents))
    except Exception as e:
        raise HTTPException(status_code=400, detail=f"CSV 讀取失敗: {str(e)}")
    if df.empty:
        raise HTTPException(status_code=400, detail="CSV 不可為空")
    return df

def filtered_csv_data(df, b_vals, a_vals, col_idx):
    if col_idx < 0 or len(df.columns) <= col_idx:
        raise HTTPException(status_code=400, detail="欄位索引超出範圍")
    col_to_filter = df.columns[col_idx]
    x_signal = pd.to_numeric(df[col_to_filter], errors="coerce")
    if x_signal.isna().any():
        raise HTTPException(status_code=400, detail=f"欄位 {col_to_filter} 含有非數值資料")
    x_values = x_signal.to_numpy(dtype=float)
    if not np.isfinite(x_values).all():
        raise HTTPException(status_code=400, detail=f"欄位 {col_to_filter} 含有非有限數值")
    y_signal = signal.lfilter(b_vals, a_vals, x_values)
    return col_to_filter, x_values, y_signal

@app.post("/api/design")
def design_filter(params: DesignParams):
    fs_val = require_finite(params.fs, "fs")
    f_type = params.filter_type
    
    try:
        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)
            dc_gain = np.sum(b_new) / np.sum(a_new)
            b_new = b_new / dc_gain
        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)
            dc_gain = np.sum(b_new) / np.sum(a_new)
            b_new = b_new / dc_gain
        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)
            dc_gain = np.sum(b_new) / np.sum(a_new)
            b_new = b_new / dc_gain
        elif f_type == "PID 控制器":
            ki_term = params.pid_ki / (2.0 * fs_val)
            kd_term = params.pid_kd * fs_val
            b0 = params.pid_kp + ki_term + kd_term
            b1 = -params.pid_kp + ki_term - 2.0 * kd_term
            b2 = kd_term
            b_new = [b0, b1, b2]
            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="未知的濾波器類型")

        b_arr, a_arr = normalize_coefficients(b_new, a_new)
            
        return {"b": b_arr.tolist(), "a": a_arr.tolist()}
    except HTTPException:
        raise
    except Exception as e:
        raise HTTPException(status_code=400, detail=str(e))

@app.post("/api/bode")
def calculate_bode(params: BodeParams):
    try:
        fs_val = require_finite(params.fs, "fs")
        b_vals = validate_coefficients(params.b, "b")
        a_vals = validate_coefficients(params.a, "a")
        f_min = params.fs / 50000.0
        f_max = params.fs * BODE_MAX_MULTIPLIER  # 修復為 Legacy 版本的顯示範圍上限 (原為 params.fs)
        f_eval = np.logspace(np.log10(f_min), np.log10(f_max), BODE_POINTS)
        
        # WorN determines frequencies to evaluate at
        w, h = signal.freqz(b_vals, a_vals, worN=f_eval, fs=fs_val)
        
        mag = 20 * np.log10(np.abs(h) + 1e-12)
        phase = np.angle(h, deg=True)
        
        return {
            "freq": f_eval.tolist(),
            "mag": mag.tolist(),
            "phase": phase.tolist()
        }
    except HTTPException:
        raise
    except Exception as e:
        raise HTTPException(status_code=400, detail=str(e))

@app.post("/api/filter")
async def filter_csv(
    file: UploadFile = File(...),
    b: str = Form(...),
    a: str = Form(...),
    col_idx: int = Form(0)
):
    try:
        b_vals = parse_coefficients(b, "b")
        a_vals = parse_coefficients(a, "a")
        df = await read_csv_upload(file)
        col_to_filter, x_signal, y_signal = filtered_csv_data(df, b_vals, a_vals, col_idx)
        index, original, filtered, step = downsample_for_plot(df.index.to_numpy(), x_signal, y_signal)

        return {
            "index": index.tolist(),
            "original": original.tolist(),
            "filtered": filtered.tolist(),
            "col_name": col_to_filter,
            "total_points": int(len(df.index)),
            "plot_points": int(len(index)),
            "downsample_step": int(step)
        }
    except HTTPException:
        raise
    except Exception as e:
        raise HTTPException(status_code=400, detail=f"CSV處理失敗: {str(e)}")

@app.post("/api/filter/download")
async def filter_csv_download(
    file: UploadFile = File(...),
    b: str = Form(...),
    a: str = Form(...),
    col_idx: int = Form(0)
):
    # 此端點專為產生包含輸出結果的 CSV 檔供下載
    try:
        b_vals = parse_coefficients(b, "b")
        a_vals = parse_coefficients(a, "a")
        df = await read_csv_upload(file)
        col_to_filter, x_signal, y_signal = filtered_csv_data(df, b_vals, a_vals, col_idx)
        
        output_col_name = f"{col_to_filter}_filtered"
        df[output_col_name] = y_signal
        
        csv_buffer = io.StringIO()
        df.to_csv(csv_buffer, index=False)
        
        return StreamingResponse(
            iter([csv_buffer.getvalue()]),
            media_type="text/csv",
            headers={"Content-Disposition": "attachment; filename=filtered_output.csv"}
        )
    except HTTPException:
        raise
    except Exception as e:
        raise HTTPException(status_code=400, detail=f"CSV處理失敗: {str(e)}")