tcad-bodeplot/dea/csv_processing.py

import io
import os
import uuid
import re

import numpy as np
import pandas as pd
from fastapi import HTTPException
from scipy import signal
import math

from .config import MAX_CSV_BYTES, MAX_PLOT_POINTS, MAX_ROWS

TEMP_CSV_DIR = "temp_csv"
os.makedirs(TEMP_CSV_DIR, exist_ok=True)

def get_cached_file_path(file_id):
    if not re.match(r'^[0-9a-f\-]{36}$', file_id):
        raise HTTPException(status_code=400, detail="無效的檔案ID")
    file_path = os.path.join(TEMP_CSV_DIR, f"{file_id}.csv")
    if not os.path.exists(file_path):
        raise HTTPException(status_code=404, detail="找不到快取的 CSV 檔案，請重新上傳")
    return file_path


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


async def save_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 不可為空")
    
    file_id = str(uuid.uuid4())
    file_path = os.path.join(TEMP_CSV_DIR, f"{file_id}.csv")
    with open(file_path, "wb") as f:
        f.write(contents)
    return file_id


def integer_lfilter(b_int, a_int, x_float, shift_in, shift_out, shift_b, shift_a, use_round=False):
    """
    全整數差分方程式模擬 (Mimic DSP hardware)
    高精度狀態變數架構：前饋不位移，保留小數精度於狀態變數中。
    支援硬體 Rounding (四捨五入) vs Floor (無條件捨去向下取整)。
    """
    b_int = np.asarray(b_int, dtype=np.int64)
    a_int = np.asarray(a_int, dtype=np.int64)
    x_int = np.round(x_float * (2**shift_in)).astype(np.int64)
    
    # y_hist 將保留在 Q_{in + b} 格式以降低 Truncation Error
    y_hist = np.zeros(len(x_int), dtype=np.int64)
    y_out = np.zeros(len(x_int), dtype=np.int64)
    
    nb = len(b_int)
    na = len(a_int)
    
    A0 = int(a_int[0])
    if A0 == 0: A0 = 1
    
    # 輸出所需的總位移量
    out_shift = shift_in + shift_b - shift_out
    
    # 預先計算四捨五入的補償值 (+0.5)
    # 韌體開發提示 (C Implementation / RISC-V):
    # 1. 標準 RISC-V (RV32I/IMAC) 的 SRA 指令是純 Floor (無條件捨去)，沒有硬體 Rounding shift。
    #    (除非具備 'P' DSP Extension 才可能有 1-cycle 的硬體 rounding shift)。
    # 2. 演算法秘技：在 C 語言中要實現 1-clock 的四捨五入，不要呼叫 float 的 round()。
    #    請使用 `y = (acc + (1 << (shift - 1))) >> shift`。
    #    編譯器會將 (1 << (shift - 1)) 編譯為常數，整體只消耗 1 個 ADD 指令，極大消除了 DC Bias。
    round_offset_a = (A0 >> 1) if use_round else 0
    round_offset_out = (1 << (out_shift - 1)) if (use_round and out_shift > 0) else 0

    for n in range(len(x_int)):
        acc = 0
        # Feedforward: 前饋完全不位移，結果為 Q_{in + b}
        for i in range(nb):
            if n - i >= 0:
                acc += b_int[i] * x_int[n - i]
                
        # Feedback: 歷史紀錄為 Q_{in + b}，係數為 Q_a，乘積為 Q_{in + b + a}
        fb_sum = 0
        for j in range(1, na):
            if n - j >= 0:
                fb_sum += a_int[j] * y_hist[n - j]
                
        # Feedback 縮放：除以 A0 (即 >> shift_a)，結果退回 Q_{in + b} 完美對齊
        # Python 的 // 等同於硬體的 SRA (Arithmetic Right Shift)，會向負無窮大 Floor
        fb_shifted = (fb_sum + round_offset_a) // A0
        
        acc -= fb_shifted
        
        # 將超高精度的 acc 直接存入歷史變數
        y_hist[n] = acc
        
        # 最終輸出再針對 Q_out 進行位移縮放
        if out_shift > 0:
            y_out[n] = (acc + round_offset_out) >> out_shift
        elif out_shift < 0:
            y_out[n] = acc << (-out_shift)
        else:
            y_out[n] = acc

    return y_out.astype(float) / (2**shift_out)


def filter_preview_response(file_id, b_vals, a_vals, col_idx, b_int=None, a_int=None, shift_in=14, shift_out=14, shift_b=14, shift_a=14, use_round=False):
    path = get_cached_file_path(file_id)
    
    # 預先讀取欄位名稱，避免用 usecols 讀取後找不到原始索引
    cols = pd.read_csv(path, nrows=0).columns.tolist()
    if col_idx < 0 or col_idx >= len(cols):
        raise HTTPException(status_code=400, detail="欄位索引超出範圍")
    col_to_filter = cols[col_idx]
    
    # 精準讀取單一欄位，並加上筆數限制 (極大降低記憶體用量與時間)
    df = pd.read_csv(path, usecols=[col_idx], nrows=MAX_ROWS)
    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} 含有非有限數值")

    # 路徑 1: 理想浮點數路徑
    y_float = signal.lfilter(b_vals, a_vals, x_values)

    # 路徑 2: 整數模擬路徑
    if b_int is not None and a_int is not None:
        y_fixed = integer_lfilter(b_int, a_int, x_values, shift_in, shift_out, shift_b, shift_a, use_round)
    else:
        y_fixed = y_float

    index, original, filtered_float, filtered_fixed, step = downsample_for_plot(df.index.to_numpy(), x_signal, y_float, y_fixed)

    return {
        "index": index.tolist(),
        "original": original.tolist(),
        "filtered": filtered_float.tolist(),
        "filtered_fixed": filtered_fixed.tolist(),
        "col_name": col_to_filter,
        "total_points": int(len(df.index)),
        "plot_points": int(len(index)),
        "downsample_step": int(step),
    }


def filtered_csv_text(file_id, b_vals, a_vals, col_idx, b_int=None, a_int=None, shift_in=14, shift_out=14, shift_b=14, shift_a=14, use_round=False):
    path = get_cached_file_path(file_id)
    # 匯出時需要原始所有欄位，但仍受限於 MAX_ROWS
    df = pd.read_csv(path, nrows=MAX_ROWS)
    
    if col_idx < 0 or col_idx >= len(df.columns):
        raise HTTPException(status_code=400, detail="欄位索引超出範圍")
    col_to_filter = df.columns[col_idx]
    x_signal = pd.to_numeric(df[col_to_filter], errors="coerce")
    x_values = x_signal.to_numpy(dtype=float)
    
    y_float = signal.lfilter(b_vals, a_vals, x_values)
    if b_int is not None and a_int is not None:
        y_fixed = integer_lfilter(b_int, a_int, x_values, shift_in, shift_out, shift_b, shift_a, use_round)
    else:
        y_fixed = y_float
        
    df[f"{col_to_filter}_filtered_ideal"] = y_float
    df[f"{col_to_filter}_filtered_fixed"] = y_fixed
    csv_buffer = io.StringIO()
    df.to_csv(csv_buffer, index=False)
    return csv_buffer.getvalue()