sva_formula/
calculus.rs

1// Concern: differentiates, Hilbert-transforms and envelopes a SpectralSum, and bounds a slope | Non-concern: the transform (fourier_dual/) | IO: (&SpectralSum) -> SpectralSum or Left; (&Body) -> f64
2
3use crate::closed_form::{Body, Edge, Part};
4use crate::complex::C64;
5use crate::fourier_dual::{dual_read, reflect_read};
6use crate::refusal::{AtomSketch, Factor, Left, LeftReason};
7use crate::spectral_sum::atom::{Exp, Factors, Gauss, Indicator, Pole, Singular, SpectralAtom};
8use crate::spectral_sum::merge::simplify;
9use crate::spectral_sum::product::times;
10use crate::spectral_sum::{Lane, SpectralSum};
11use crate::through::Reads;
12
13/// The square of the analytic signal's modulus, which is in A; the root is the observation's
14/// job, so A stays closed.
15#[derive(Clone, Debug, PartialEq)]
16pub struct Envelope {
17    pub squared: SpectralSum,
18}
19
20/// A is closed under `d/dx` up to an order or degree past `u16`, which refuses. A series
21/// differentiates termwise, lazily where its term is not one of the shapes the Fourier dual rules read.
22pub fn d_dt(n: &SpectralSum) -> Result<SpectralSum, Left> {
23    let mut lanes = Vec::with_capacity(n.lanes.len());
24    for lane in &n.lanes {
25        let mut out = Lane {
26            atoms: Vec::new(),
27            series: lane.series.iter().map(derive_series).collect(),
28            modal: Vec::new(),
29        };
30        for atom in &lane.atoms {
31            out.atoms.extend(atom.derivative()?);
32        }
33        for bank in &lane.modal {
34            let expanded = crate::modal::atoms(bank, crate::origin::Origin::UNKNOWN);
35            for atom in &expanded {
36                out.atoms.extend(atom.derivative()?);
37            }
38        }
39        simplify(&mut out);
40        lanes.push(out);
41    }
42    Ok(SpectralSum::of(n.var, lanes))
43}
44
45/// At least `sup |f'(t)|` over every instant: each atom's own sup, summed. `None` where an
46/// atom's slope grows without bound either way in time, or `f` has no atom sum.
47pub fn steepest(f: &Body) -> Option<f64> {
48    steepest_read(f, &crate::through::Opaque)
49}
50
51/// The same, each ref read as `reads` answers it.
52pub fn steepest_read(f: &Body, reads: &dyn crate::through::Reads) -> Option<f64> {
53    let normalized =
54        crate::spectral_sum::build::normalize_read(f, crate::closed_form::Var::T, reads);
55    let slope = d_dt(&normalized.ok()?).ok()?;
56    let mut held = 0.0;
57    for lane in &slope.lanes {
58        if !lane.series.is_empty() || !lane.modal.is_empty() {
59            return None;
60        }
61        for atom in &lane.atoms {
62            held += crate::spectral_sum::sup::magnitude_upper_bound_from_instant(
63                atom,
64                f64::NEG_INFINITY,
65            )?;
66        }
67    }
68    held.is_finite().then_some(held)
69}
70
71fn derive_series(s: &crate::closed_form::Series) -> crate::closed_form::Series {
72    crate::closed_form::Series {
73        term: Part::new(
74            s.term.origin,
75            Body::Deriv {
76                order: 1,
77                of: s.term.clone(),
78            },
79        ),
80        ..s.clone()
81    }
82}
83
84/// Dual, multiply by `-i*sgn`, dual, reflect. It is never a convolution with `pv(t)/pi`.
85pub fn hilbert(n: &SpectralSum) -> Result<SpectralSum, Left> {
86    hilbert_read(n, &crate::through::Opaque)
87}
88
89pub fn hilbert_read(n: &SpectralSum, reads: &dyn Reads) -> Result<SpectralSum, Left> {
90    let spectrum = dual_read(n, reads)?;
91    let mut lanes = Vec::with_capacity(spectrum.lanes.len());
92    for lane in &spectrum.lanes {
93        if !lane.is_finite_sum() {
94            return Err(Left::new(
95                crate::origin::Origin::UNKNOWN,
96                AtomSketch::of(Factor::Value),
97                LeftReason::SeriesNonUniform,
98            ));
99        }
100        let mut atoms = Vec::new();
101        for a in &lane.atoms {
102            atoms.extend(times_signum(a)?);
103        }
104        let mut out = Lane::of(atoms);
105        simplify(&mut out);
106        lanes.push(out);
107    }
108    reflect_read(
109        &dual_read(&SpectralSum::of(spectrum.var, lanes), reads)?,
110        reads,
111    )
112}
113
114/// `sgn` is smooth away from the origin, so it reads a delta's own sign; a delta AT the
115/// origin is `H{1} = 0`, and its derivatives are the duals of the polynomials.
116fn times_signum(a: &SpectralAtom) -> Result<Vec<SpectralAtom>, Left> {
117    let minus_i = C64::new(0.0, -1.0);
118    if let Singular::Delta { at, order } = a.sing {
119        if at == 0.0 {
120            return match order {
121                0 => Ok(Vec::new()),
122                _ => Err(Left::new(
123                    a.origin,
124                    AtomSketch::of(Factor::Delta),
125                    LeftReason::HilbertOfPolynomial,
126                )),
127            };
128        }
129        let sign = if at > 0.0 { 1.0 } else { -1.0 };
130        return Ok(vec![SpectralAtom::new(
131            a.c * minus_i.scale(sign),
132            Factors::NONE,
133            a.sing,
134            a.origin,
135        )]);
136    }
137    if a.pole.is_some() {
138        return Err(Left::new(
139            a.origin,
140            AtomSketch::pair(
141                if a.pole.is_some_and(|p| p.pv) {
142                    Factor::PrincipalValue
143                } else {
144                    Factor::Pole
145                },
146                Factor::Indicator,
147            ),
148            LeftReason::PoleTimesIndicator,
149        ));
150    }
151    if a.gauss.is_some() {
152        return Err(Left::new(
153            a.origin,
154            AtomSketch::pair(Factor::Gaussian, Factor::Indicator),
155            LeftReason::GaussianTimesIndicator,
156        ));
157    }
158    let halves = [
159        (1.0, Edge::at(0.0), Edge::PosInf),
160        (-1.0, Edge::NegInf, Edge::at(0.0)),
161    ];
162    Ok(halves
163        .into_iter()
164        .filter_map(|(sign, l, r)| {
165            let window = a
166                .ind
167                .map_or(Indicator { l, r }, |w| w.meet(Indicator { l, r }));
168            (!window.is_empty()).then(|| {
169                SpectralAtom::new(
170                    a.c * minus_i.scale(sign),
171                    Factors {
172                        ind: Some(window),
173                        ..a.factors()
174                    },
175                    Singular::Regular,
176                    a.origin,
177                )
178            })
179        })
180        .collect())
181}
182
183pub fn analytic(n: &SpectralSum) -> Result<SpectralSum, Left> {
184    analytic_read(n, &crate::through::Opaque)
185}
186
187pub fn analytic_read(n: &SpectralSum, reads: &dyn Reads) -> Result<SpectralSum, Left> {
188    let quadrature = hilbert_read(n, reads)?;
189    let mut lanes = Vec::with_capacity(n.lanes.len());
190    for (real, imaginary) in n.lanes.iter().zip(quadrature.lanes.iter()) {
191        let mut out = Lane::of(
192            real.atoms
193                .iter()
194                .copied()
195                .chain(
196                    imaginary
197                        .atoms
198                        .iter()
199                        .map(|a| a.with(a.c * C64::I, a.factors())),
200                )
201                .collect(),
202        );
203        simplify(&mut out);
204        lanes.push(out);
205    }
206    Ok(SpectralSum::of(n.var, lanes))
207}
208
209pub fn envelope(n: &SpectralSum) -> Result<Envelope, Left> {
210    envelope_read(n, &crate::through::Opaque)
211}
212
213pub fn envelope_read(n: &SpectralSum, reads: &dyn Reads) -> Result<Envelope, Left> {
214    let signal = analytic_read(n, reads)?;
215    let mut lanes = Vec::with_capacity(signal.lanes.len());
216    for lane in &signal.lanes {
217        let mut atoms = Vec::new();
218        for a in &lane.atoms {
219            for b in &lane.atoms {
220                atoms.extend(times(a, &conjugate(b))?);
221            }
222        }
223        let mut out = Lane::of(atoms);
224        simplify(&mut out);
225        lanes.push(out);
226    }
227    Ok(Envelope {
228        squared: SpectralSum::of(signal.var, lanes),
229    })
230}
231
232fn conjugate(a: &SpectralAtom) -> SpectralAtom {
233    if a.is_delta() {
234        return SpectralAtom::new(a.c.conj(), Factors::NONE, a.sing, a.origin);
235    }
236    let f = a.factors();
237    SpectralAtom::new(
238        a.c.conj(),
239        Factors {
240            exp: f.exp.map(|e| Exp {
241                omega: -e.omega,
242                ..e
243            }),
244            gauss: f.gauss.map(|g| Gauss { ..g }),
245            pole: f.pole.map(|p| Pole {
246                at: p.at.conj(),
247                ..p
248            }),
249            ..f
250        },
251        Singular::Regular,
252        a.origin,
253    )
254}
sva_formula/calculus.rs

sva_formula/
calculus.rs
