scientific-cal 0.2.4

use std::f64;
use crate::filters::{base::FilterBand, iir::lib::{besselap, butterap, cheb1ap, cheb2ap, zpk2sos, bilinear_zpk, lp2hp_zpk, lp2lp_zpk, lp2bp_zpk, lp2bs_zpk}, BesselCriterion};
/// iir 型滤波器
///
/// # 参数
///
/// - `N`：order
/// - `Wn`: cutoff / (fs / 2)。
/// - `btype`：'lowpass', 'highpass'。
/// - `analog`：默认false   true 启动模拟滤波器
/// - `ftype`：'phase', 'delay', 'mag'
/// - `fs`：
///
/// # 返回值
/// 
/// 滤波器设计参数 Vec<Sos>。
/// 
/// # 示例
///
/// ```
/// 无
/// ```
pub fn iir_filter(
    n: usize, 
    wn: &mut Vec<f64>, 
    btype: FilterBand, 
    analog: bool, 
    ftype: String, 
    mut fs: Option<f64>, 
    rp: Option<f64>,
    rs: Option<f64>
) -> Result<Vec<[f64; 6]>, String>{
    // Normalize Wn if fs is provided
    if let Some(fs_val) = fs {
        if analog {
            return Err("fs cannot be specified for an analog filter".to_string()); // fs 不能指定为模拟滤波器
        }
        for w in wn.iter_mut() {
            *w /= fs_val / 2.0;
        }
    }

    // Validate Wn
    if wn.iter().any(|&w| w < 0.0) {
        return Err("filter critical frequencies must be greater than 0".to_string()); // 滤波器临界频率必须大于 0
    }

    if wn.len() > 1 && wn[0] >= wn[1] {
        return Err("Wn[0] must be less than Wn[1]".to_string());   // Wn[0] 必须小于 Wn[1]
    }
    // Placeholder: analog lowpass prototype (besselap)
    let (mut z, mut p, mut k) = match ftype.split('_').collect::<Vec<_>>().as_slice() {
        ["bessel", "delay"] => besselap(n, BesselCriterion::Delay)?,
        ["bessel", "phase"] => besselap(n, BesselCriterion::Phase)?,
        ["bessel", "mag"] => besselap(n, BesselCriterion::Mag)?,
        ["butter"] => butterap(n)?,
        ["cheby1"] => {
            let p = if let Some(d) = rp {
                d
            }else{
               return Err("passband ripple (rp) must be provided to design a Chebyshev I filter.".to_string())
            };
            if p < 0.0 {
                return Err("passband ripple (rp) must be positive".to_string());
            }
            cheb1ap(n, p)?
        },
        ["cheby2"] => {
            let s = if let Some(d) = rs {
                d
            }else{
               return Err("stopband attenuation (rs) must be provided to design an Chebyshev II filter.".to_string())
            };
            if s < 0.0 {
                return Err("stopband attenuation (rs) must be positive".to_string());
            }
            cheb2ap(n, s)?
        },
        _ => return Err("filter error".to_string()),
    };
    let warped = if !analog {
        if wn.iter().any(|&w| w < 0.0 || w >= 1.0) {
            if let Some(fs_val) = fs {
                return Err(format!(
                    "Digital filter critical frequencies must be 0 < Wn < fs/2 (fs={} -> fs/2={})",  // 数字滤波器临界频率必须为 0 < Wn < fs/2
                    fs_val,
                    fs_val / 2.0
                ));
            }
            return Err("Digital filter critical frequencies must be 0 < Wn < 1".to_string()); // 数字滤波器临界频率必须为 0 < Wn < 1
        }
        fs = Some(2.0);
        wn.iter().map(|&w| 2.0 * 2.0 * (std::f64::consts::PI * w / 2.0).tan()).collect::<Vec<_>>()
    } else {
        wn.clone()
    };

    // Transform filter type
    match btype {
        FilterBand::LowPass => {
            if wn.len() != 1 {
                return Err("Must specify a single critical frequency for lowpass or highpass filter".to_string());  // 必须为低通或高通滤波器指定单个临界频率
            }
            (z, p, k) = lp2lp_zpk(&z, &p, k, warped[0]);
        }
        FilterBand::HighPass => {
            if wn.len() != 1 {
                return Err("Must specify a single critical frequency for lowpass or highpass filter".to_string());  // 必须为低通或高通滤波器指定单个临界频率
            }
            (z, p, k) = lp2hp_zpk(&z, &p, k, warped[0]);
        }
        FilterBand::BandPass => {
            if wn.len() != 2 {
                return Err("Must specify two single critical frequency for BandPass or BandStop filter".to_string());  
            }
            let bw = warped[1] - warped[0];
            let wo = (warped[0] * warped[1]).sqrt();
            (z, p, k) = lp2bp_zpk(&z, &p, k, wo, bw);
        },
        FilterBand::BandStop => {
            if wn.len() != 2 {
                return Err("Must specify two single critical frequency for BandPass or BandStop filter".to_string());  
            }
            let bw = warped[1] - warped[0];
            let wo = (warped[0] * warped[1]).sqrt();
            (z, p, k) = lp2bs_zpk(&z, &p, k, wo, bw);
        },
    }
    // Placeholder: Bilinear transform
    if !analog {
        (z, p, k) = bilinear_zpk(&z, &p, k, fs);
    }
    // Output format
    Ok(zpk2sos(&z, &p, k, analog, None)?)
}