前腕8ch表面筋電位（sEMG）と手話非手指動作（NMM）のリアルタイム信号処理：
双二次IIRノッチフィルタと256点STFT実装

#![no_std]
#![forbid(unsafe_code)]

pub const SEMG_CHANNELS: usize = 8;
pub const STFT_WINDOW_LEN: usize = 256;

/// Normalized Biquad IIR filter coefficients (a0 normalized to 1.0).
#[derive(Clone, Copy, Debug)]
pub struct BiquadCoeffs {
    pub b0: f32,
    pub b1: f32,
    pub b2: f32,
    pub a1: f32,
    pub a2: f32,
}

impl BiquadCoeffs {
    /// Precomputed 60.0 Hz Notch Filter at f_s = 1000.0 Hz, Q = 30.0 (West Japan / Takatsuki)
    pub const NOTCH_60HZ_FS1000_Q30: Self = Self {
        b0: 0.993_912_4,
        b1: -1.847_590_6,
        b2: 0.993_912_4,
        a1: -1.847_590_6,
        a2: 0.987_824_8,
    };

    /// Precomputed 50.0 Hz Notch Filter at f_s = 1000.0 Hz, Q = 30.0 (East Japan)
    pub const NOTCH_50HZ_FS1000_Q30: Self = Self {
        b0: 0.994_885_1,
        b1: -1.892_384_5,
        b2: 0.994_885_1,
        a1: -1.892_384_5,
        a2: 0.989_770_2,
    };
}

/// Direct Form II Transposed state registers for 8 synchronous sEMG channels.
#[derive(Clone, Copy, Debug)]
pub struct MultiChannelBiquad {
    coeffs: BiquadCoeffs,
    w1: [f32; SEMG_CHANNELS],
    w2: [f32; SEMG_CHANNELS],
}

/// Real-time extracted feature frame for Sign Language Handshape & NMM tension inference.
#[derive(Clone, Copy, Debug)]
pub struct SemgFeatureFrame {
    pub mean_abs_value: [f32; SEMG_CHANNELS],
    pub waveform_length: [f32; SEMG_CHANNELS],
    pub zero_crossings: [u16; SEMG_CHANNELS],
    pub isometric_co_contraction: f32,
}

impl MultiChannelBiquad {
    pub const fn new(coeffs: BiquadCoeffs) -> Self {
        Self {
            coeffs,
            w1: [0.0; SEMG_CHANNELS],
            w2: [0.0; SEMG_CHANNELS],
        }
    }

    /// Processes a single 8-channel sample vector `x[0..8]` in O(1) time using DF-II Transposed.
    #[inline]
    pub fn process_sample_8ch(&mut self, input: &[f32; SEMG_CHANNELS]) -> [f32; SEMG_CHANNELS] {
        let mut out = [0.0f32; SEMG_CHANNELS];
        let mut ch = 0;
        while ch < SEMG_CHANNELS {
            let x = input[ch];
            let y = self.coeffs.b0 * x + self.w1[ch];
            self.w1[ch] = self.coeffs.b1 * x - self.coeffs.a1 * y + self.w2[ch];
            self.w2[ch] = self.coeffs.b2 * x - self.coeffs.a2 * y;
            out[ch] = y;
            ch += 1;
        }
        out
    }
}

/// Applies a 256-point window and extracts MAV, Waveform Length (WL), Zero Crossings (ZC),
/// and antagonist co-contraction index for NMM grammatical intensity estimation.
pub fn extract_window_features_256(
    filtered_window: &[[f32; SEMG_CHANNELS]; STFT_WINDOW_LEN],
    hanning_weights: &[f32; STFT_WINDOW_LEN],
    zc_deadzone_mv: f32,
) -> SemgFeatureFrame {
    let mut mav = [0.0f32; SEMG_CHANNELS];
    let mut wl = [0.0f32; SEMG_CHANNELS];
    let mut zc = [0u16; SEMG_CHANNELS];

    for n in 0..STFT_WINDOW_LEN {
        let w = hanning_weights[n];
        for ch in 0..SEMG_CHANNELS {
            let s_curr = filtered_window[n][ch] * w;
            let abs_s = if s_curr < 0.0 { -s_curr } else { s_curr };
            mav[ch] += abs_s;

            if n > 0 {
                let s_prev = filtered_window[n - 1][ch] * hanning_weights[n - 1];
                let diff = s_curr - s_prev;
                wl[ch] += if diff < 0.0 { -diff } else { diff };

                let sign_change = (s_curr > 0.0 && s_prev < 0.0) || (s_curr < 0.0 && s_prev > 0.0);
                let amp_diff = if diff < 0.0 { -diff } else { diff };
                if sign_change && amp_diff >= zc_deadzone_mv {
                    zc[ch] += 1;
                }
            }
        }
    }

    let inv_n = 1.0 / (STFT_WINDOW_LEN as f32);
    let mut flexor_sum = 0.0f32;
    let mut extensor_sum = 0.0f32;
    for ch in 0..SEMG_CHANNELS {
        mav[ch] *= inv_n;
        if ch < 4 {
            flexor_sum += mav[ch];
        } else {
            extensor_sum += mav[ch];
        }
    }

    let denom = flexor_sum + extensor_sum + 1e-6;
    let min_group = if flexor_sum < extensor_sum { flexor_sum } else { extensor_sum };
    let isometric_co_contraction = (2.0 * min_group) / denom;

    SemgFeatureFrame {
        mean_abs_value: mav,
        waveform_length: wl,
        zero_crossings: zc,
        isometric_co_contraction,
    }
}