Skip to main content

qualia_core_db/solvers/transforms/
ztransform.rs

1//! Z-transform `X(z) = Σ_{n≥0} x[n] z^{−n}` of a finite causal sequence, evaluated at a
2//! complex `z`, plus the standard closed forms for the unit step and the geometric
3//! sequence.
4
5use super::fourier::Cplx;
6
7#[inline]
8fn cmul(a: Cplx, b: Cplx) -> Cplx {
9    (a.0 * b.0 - a.1 * b.1, a.0 * b.1 + a.1 * b.0)
10}
11/// Complex reciprocal `1/z`. `None` at `z = 0`.
12fn cinv(z: Cplx) -> Option<Cplx> {
13    let d = z.0 * z.0 + z.1 * z.1;
14    if d == 0.0 {
15        return None;
16    }
17    Some((z.0 / d, -z.1 / d))
18}
19fn csub(a: Cplx, b: Cplx) -> Cplx {
20    (a.0 - b.0, a.1 - b.1)
21}
22
23/// `X(z) = Σ_{n=0}^{N−1} x[n] z^{−n}` for a finite real sequence. `None` at `z = 0`.
24pub fn z_transform_finite(x: &[f64], z: Cplx) -> Option<Cplx> {
25    let zinv = cinv(z)?;
26    let mut acc = (0.0, 0.0);
27    let mut zpow = (1.0, 0.0); // z^{-n}, starts at n=0
28    for &xn in x {
29        acc = (acc.0 + xn * zpow.0, acc.1 + xn * zpow.1);
30        zpow = cmul(zpow, zinv);
31    }
32    Some(acc)
33}
34
35/// Closed form for the unit step `u[n]`: `X(z) = 1/(1 − z^{−1}) = z/(z−1)`, valid for
36/// `|z| > 1`. `None` at `z = 0` or `z = 1`.
37pub fn unit_step_z(z: Cplx) -> Option<Cplx> {
38    let zinv = cinv(z)?;
39    let denom = csub((1.0, 0.0), zinv); // 1 − z^{-1}
40    let id = cinv(denom)?;
41    Some(id)
42}
43
44/// Closed form for `a^n u[n]`: `X(z) = 1/(1 − a·z^{−1})`, valid for `|z| > |a|`.
45pub fn geometric_z(a: f64, z: Cplx) -> Option<Cplx> {
46    let zinv = cinv(z)?;
47    let azinv = (a * zinv.0, a * zinv.1);
48    let denom = csub((1.0, 0.0), azinv);
49    cinv(denom)
50}
51
52#[cfg(test)]
53mod tests {
54    use super::*;
55    const EPS: f64 = 1e-9;
56
57    #[test]
58    fn finite_sequence_evaluates() {
59        // x = [1,2,3] at z = 2 (real): 1 + 2/2 + 3/4 = 2.75
60        let v = z_transform_finite(&[1.0, 2.0, 3.0], (2.0, 0.0)).unwrap();
61        assert!((v.0 - 2.75).abs() < EPS && v.1.abs() < EPS);
62        // delta[n] = [1] → X(z) = 1 everywhere
63        let d = z_transform_finite(&[1.0], (3.0, -1.0)).unwrap();
64        assert!((d.0 - 1.0).abs() < EPS && d.1.abs() < EPS);
65        assert!(z_transform_finite(&[1.0], (0.0, 0.0)).is_none());
66    }
67
68    #[test]
69    fn closed_forms_match_truncated_sums() {
70        // Geometric a=0.5: closed form ≈ truncated finite sum for |z|>|a|.
71        let z = (2.0, 0.0);
72        let closed = geometric_z(0.5, z).unwrap();
73        let seq: Vec<f64> = (0..60).map(|n| 0.5_f64.powi(n)).collect();
74        let approx = z_transform_finite(&seq, z).unwrap();
75        assert!((closed.0 - approx.0).abs() < 1e-6 && (closed.1 - approx.1).abs() < 1e-6);
76        // Unit step closed form vs truncated.
77        let us = unit_step_z(z).unwrap();
78        let ones = vec![1.0; 60];
79        let ua = z_transform_finite(&ones, z).unwrap();
80        assert!((us.0 - ua.0).abs() < 1e-6);
81    }
82}