Coverage for src/chebpy/trigtech.py: 100%
320 statements
« prev ^ index » next coverage.py v7.15.2, created at 2026-07-22 21:50 +0000
1"""Trigonometric (Fourier) technology for periodic function approximation.
3This module provides the Trigtech class, which represents smooth periodic functions
4on [-1, 1] using truncated Fourier series. It is the trigonometric analogue of
5Chebtech and sits in the same class hierarchy:
7 Onefun → Smoothfun → Trigtech
9Coefficient storage convention (NumPy-native / FFT order)
10----------------------------------------------------------
11Given n equispaced sample points x_j = -1 + 2j/n (j = 0, …, n-1), the stored
12coefficients are
14 coeffs[k] = (1/n) * sum_j f(x_j) * exp(-2*pi*i*j*k/n)
15 = (numpy.fft.fft(values) / n)[k]
17This is exactly the output of ``numpy.fft.fft(values) / n``, i.e. NumPy-native
18(FFT) ordering: DC at index 0, positive frequencies 1 … n//2, then negative
19frequencies -(n//2)+1 … -1.
21Use ``_coeffs_to_plotorder()`` to obtain the human-readable DC-centred ordering
22(equivalent to ``numpy.fft.fftshift``).
24Evaluation
25----------
26Any point x ∈ [-1, 1] is evaluated via the DFT summation formula:
28 f(x) = Σ_k coeffs[k] * exp(i*π*ω_k*(x+1))
30where ω_k = numpy.fft.fftfreq(n)*n gives the integer frequencies in FFT order.
32References:
33----------
34* Trefethen, "Spectral Methods in MATLAB" (SIAM 2000)
35* Chebfun @trigtech (github.com/chebfun/chebfun)
36"""
38import warnings
39from abc import ABC
40from typing import Any, cast
42import matplotlib.pyplot as plt
43import numpy as np
45from .decorators import self_empty
46from .plotting import plotfun, plotfuncoeffs
47from .settings import _preferences as prefs
48from .smoothfun import Smoothfun
49from .utilities import Interval, coerce_list
52def _trig_adaptive(
53 cls: Any,
54 fun: Any,
55 hscale: float = 1,
56 maxpow2: int | None = None,
57) -> np.ndarray:
58 """Adaptively determine the Fourier coefficients needed to represent *fun*.
60 Uses successively finer equispaced grids (sizes 2**k) until the
61 high-frequency Fourier modes decay below tolerance. Convergence is
62 assessed via the one-sided symmetric maximum of the DC-centred coefficient
63 magnitudes: ``abs_sym[k] = max(|c_k|, |c_{-k}|) / vscale``. The series
64 is considered converged when the Nyquist/highest-frequency mode
65 ``abs_sym[-1]`` falls below *tol*.
67 Args:
68 cls: Trigtech class (provides ``_trigpts`` and ``_vals2coeffs``).
69 fun: Callable to approximate.
70 hscale: Horizontal scale for tolerance adjustment.
71 maxpow2: Maximum power of 2 to try (defaults to ``prefs.maxpow2``).
73 Returns:
74 numpy.ndarray: Fourier coefficients in NumPy FFT order.
75 """
76 minpow2 = 3 # start at n = 8
77 maxpow2 = maxpow2 if maxpow2 is not None else prefs.maxpow2
78 tol = prefs.eps * max(hscale, 1)
79 coeffs: np.ndarray = np.array([])
80 for k in range(minpow2, max(minpow2, maxpow2) + 1):
81 n = 2**k
82 points = cls._trigpts(n)
83 values = fun(points)
84 coeffs = cls._vals2coeffs(values)
85 vscale = float(np.max(np.abs(values)))
86 if vscale <= tol:
87 return np.array([0.0])
89 # Build one-sided symmetric maximum:
90 # abs_sym[ki] = max(|c_{ki}|, |c_{-ki}|) / vscale for ki = 0…n//2
91 centered = np.fft.fftshift(coeffs)
92 dc_idx = n // 2
93 abs_sym = np.zeros(dc_idx + 1)
94 for ki in range(dc_idx + 1):
95 p = centered[dc_idx + ki] if dc_idx + ki < n else 0.0
96 q = centered[dc_idx - ki]
97 abs_sym[ki] = max(abs(p), abs(q)) / vscale
99 # Convergence: the Nyquist/highest-frequency mode is negligible.
100 if abs_sym[-1] <= tol:
101 above = np.where(abs_sym > tol)[0]
102 if len(above) == 0: # pragma: no cover
103 # Defensive: the normalised peak coefficient is >= 1/n >> tol
104 # whenever vscale > tol, so 'above' is never empty here.
105 return np.array([0.0])
106 max_k = int(above[-1]) # highest significant frequency index
107 start = dc_idx - max_k
108 end = dc_idx + max_k + 1
109 return np.fft.ifftshift(centered[start:end])
111 if k == maxpow2:
112 warnings.warn(
113 f"The {cls.__name__} constructor did not converge: using {n} points",
114 stacklevel=3,
115 )
116 break
117 return coeffs
120class Trigtech(Smoothfun, ABC):
121 """Trigonometric (Fourier) function approximation on [-1, 1].
123 Represents a smooth periodic function f: [-1, 1] -> R (or C) as a
124 truncated Fourier series. Coefficients are stored in NumPy FFT order;
125 see module docstring for the precise convention.
127 This class is ``ABC`` so that it cannot be instantiated directly—exactly
128 mirroring Chebtech, which is also abstract (concrete only through the
129 ``Chebtech`` name used everywhere). In practice ``Trigtech`` is both the
130 abstract base and the concrete class: it is not further subclassed, but
131 the ABC marker prevents accidental bare construction without going through
132 a named constructor.
133 """
135 # ------------------------------------------------------------------
136 # alternative constructors
137 # ------------------------------------------------------------------
139 @classmethod
140 def initconst(cls, c: Any = None, *, interval: Any = None) -> "Trigtech":
141 """Initialise a Trigtech from a constant *c*."""
142 if not np.isscalar(c):
143 raise ValueError(c)
144 if isinstance(c, int):
145 c = float(c)
146 return cls(np.array([c]), interval=interval)
148 @classmethod
149 def initempty(cls, *, interval: Any = None) -> "Trigtech":
150 """Initialise an empty Trigtech."""
151 return cls(np.array([]), interval=interval)
153 @classmethod
154 def initidentity(cls, *, interval: Any = None) -> "Trigtech":
155 """Trigtech approximation of the identity f(x) = x on [-1, 1].
157 Note: f(x) = x is *not* periodic on [-1, 1], so this will not converge
158 to machine precision. It is provided for interface compatibility with
159 Chebtech; in practice ``Classicfun.initidentity`` is used instead.
160 """
161 interval = interval if interval is not None else prefs.domain
162 return cls.initfun_adaptive(lambda x: x, interval=interval)
164 @classmethod
165 def initfun(cls, fun: Any = None, n: Any = None, *, interval: Any = None) -> "Trigtech":
166 """Convenience constructor: adaptive if *n* is None, fixed-length otherwise."""
167 if n is None:
168 return cls.initfun_adaptive(fun, interval=interval)
169 return cls.initfun_fixedlen(fun, n, interval=interval)
171 @classmethod
172 def initfun_fixedlen(cls, fun: Any = None, n: Any = None, *, interval: Any = None) -> "Trigtech":
173 """Initialise a Trigtech from callable *fun* using *n* equispaced points."""
174 if n is None:
175 raise ValueError("initfun_fixedlen requires the n parameter to be specified") # noqa: TRY003
176 points = cls._trigpts(int(n))
177 values = fun(points)
178 coeffs = cls._vals2coeffs(values)
179 return cls(coeffs, interval=interval)
181 @classmethod
182 def initfun_adaptive(cls, fun: Any = None, *, interval: Any = None) -> "Trigtech":
183 """Initialise a Trigtech from callable *fun* using the adaptive constructor."""
184 interval = interval if interval is not None else prefs.domain
185 interval = Interval(*interval)
186 coeffs = _trig_adaptive(cls, fun, hscale=interval.hscale)
187 return cls(coeffs, interval=interval)
189 @classmethod
190 def initvalues(cls, values: Any = None, *, interval: Any = None) -> "Trigtech":
191 """Initialise a Trigtech from function values at equispaced points."""
192 return cls(cls._vals2coeffs(np.asarray(values)), interval=interval)
194 # ------------------------------------------------------------------
195 # core dunder methods
196 # ------------------------------------------------------------------
198 def __init__(self, coeffs: Any, interval: Any = None) -> None:
199 """Initialise a Trigtech with FFT-order *coeffs* on *interval*.
201 Coefficients are always stored as complex128. The :attr:`iscomplex`
202 property returns True only when the function *values* are complex
203 (i.e., the coefficients do **not** satisfy the conjugate-symmetry
204 condition C_{n-k} ≈ conj(C_k)).
206 Args:
207 coeffs: 1-D array of Fourier coefficients in NumPy FFT order.
208 interval: Two-element interval [a, b]. Defaults to ``prefs.domain``.
209 """
210 interval = interval if interval is not None else prefs.domain
211 self._coeffs = np.array(coeffs, dtype=complex)
212 self._interval = Interval(*interval)
214 def __call__(self, x: Any, how: str = "fft") -> Any: # noqa: ARG002 (how kept for Chebtech interface parity)
215 """Evaluate the Trigtech at points *x* via the DFT summation formula.
217 f(x) = Σ_k coeffs[k] * exp(i*π*ω_k*(x+1))
219 where ω_k = fftfreq(n)*n gives integer frequencies in FFT order.
220 For real-valued functions the imaginary part of the result is discarded.
222 Args:
223 x: Evaluation points in [-1, 1].
224 how: Ignored; present for interface compatibility with Chebtech.
225 """
226 if self.isempty:
227 return np.array([])
228 scalar = np.isscalar(x)
229 x = np.atleast_1d(np.asarray(x, dtype=float)).ravel()
231 if self.isconst:
232 c0 = self._coeffs[0].real if not self.iscomplex else self._coeffs[0]
233 out = c0 * np.ones(x.size)
234 return float(out[0]) if scalar else out
236 n = self.size
237 freqs = np.fft.fftfreq(n) * n # [0, 1, …, n//2, -(n//2)+1, …, -1]
238 # shape: (len(x), n) @ (n,) → (len(x),)
239 phases = np.exp(1j * np.pi * np.outer(x + 1.0, freqs))
240 result = phases @ self._coeffs
241 if not self.iscomplex:
242 result = result.real
243 return float(result[0]) if scalar else result
245 def __repr__(self) -> str: # pragma: no cover
246 """Return a concise string representation."""
247 return f"<{self.__class__.__name__}{{{self.size}}}>"
249 # ------------------------------------------------------------------
250 # properties
251 # ------------------------------------------------------------------
253 @property
254 def coeffs(self) -> np.ndarray:
255 """Fourier coefficients in NumPy FFT order (always complex128)."""
256 return self._coeffs
258 @property
259 def interval(self) -> Interval:
260 """Interval that the Trigtech is mapped to."""
261 return self._interval
263 @property
264 def size(self) -> int:
265 """Number of stored Fourier coefficients."""
266 return self._coeffs.size
268 @property
269 def isempty(self) -> bool:
270 """True if the Trigtech has no coefficients."""
271 return self.size == 0
273 @property
274 def iscomplex(self) -> bool:
275 """True if the function is complex-valued (values have a non-negligible imaginary part).
277 This is determined by checking whether the Fourier coefficients violate
278 the conjugate-symmetry condition C_{n-k} ≈ conj(C_k) that holds for
279 every real-valued periodic function.
280 """
281 n = self.size
282 if n <= 1:
283 return bool(np.any(np.abs(np.imag(self._coeffs)) > 0))
284 abs_max = float(np.max(np.abs(self._coeffs)))
285 if abs_max == 0.0:
286 return False
287 tol = 1e-8 * abs_max
288 # mirror[k-1] = conj(C_{n-k}) for k = 1,...,n-1
289 mirror = np.conj(self._coeffs[-1:0:-1])
290 return bool(np.any(np.abs(self._coeffs[1:] - mirror) > tol))
292 @property
293 def isconst(self) -> bool:
294 """True if the Trigtech represents a constant (single coefficient)."""
295 return self.size == 1
297 @property
298 def isperiodic(self) -> bool:
299 """Always True: Trigtech always represents a periodic function."""
300 return True
302 @property
303 @self_empty(0.0)
304 def vscale(self) -> float:
305 """Estimate the vertical scale (max |f|)."""
306 return float(np.abs(np.asarray(coerce_list(self.values()))).max())
308 # ------------------------------------------------------------------
309 # utilities
310 # ------------------------------------------------------------------
312 def copy(self) -> "Trigtech":
313 """Return a deep copy."""
314 return self.__class__(self._coeffs.copy(), interval=self._interval.copy())
316 def imag(self) -> "Trigtech":
317 """Return the imaginary part of the function as a real-valued Trigtech.
319 For a complex function f(x) = g(x) + i·h(x), the Fourier coefficients
320 of h(x) are H[k] = (D[k] - conj(D[n-k])) / (2i) for k ≥ 1,
321 and H[0] = Im(D[0]).
322 """
323 if not self.iscomplex:
324 return self.initconst(0.0, interval=self._interval)
325 n = self.size
326 c = self._coeffs
327 imag_c = np.zeros(n, dtype=complex)
328 imag_c[0] = np.imag(c[0])
329 if n > 1:
330 mirror = np.conj(c[-1:0:-1]) # conj(c[n-1]), ..., conj(c[1])
331 imag_c[1:] = (c[1:] - mirror) / (2j)
332 return self.__class__(imag_c, self._interval)
334 def prolong(self, n: int) -> "Trigtech":
335 """Return a Trigtech of length *n* (truncate or zero-pad in frequency space).
337 The operation aligns DC components of the source and target DC-centred
338 representations, then either pads with zeros (n > m) or slices (n < m).
339 This correctly handles the asymmetry between even- and odd-length arrays.
340 """
341 m = self.size
342 if n == m:
343 return self.copy()
345 centered = np.fft.fftshift(self._coeffs)
346 dc_src = m // 2
347 dc_tgt = n // 2
349 if n > m:
350 padded = np.zeros(n, dtype=centered.dtype)
351 start = dc_tgt - dc_src
352 padded[start : start + m] = centered
353 return self.__class__(np.fft.ifftshift(padded), interval=self._interval)
354 else:
355 start = dc_src - dc_tgt
356 truncated = centered[start : start + n]
357 return self.__class__(np.fft.ifftshift(truncated), interval=self._interval)
359 def real(self) -> "Trigtech":
360 """Return the real part of the function as a real-valued Trigtech.
362 For a complex function f(x) = g(x) + i·h(x), the Fourier coefficients
363 of g(x) are G[k] = (D[k] + conj(D[n-k])) / 2 for k ≥ 1,
364 and G[0] = Re(D[0]).
365 """
366 if not self.iscomplex:
367 return self
368 n = self.size
369 c = self._coeffs
370 real_c = np.zeros(n, dtype=complex)
371 real_c[0] = np.real(c[0])
372 if n > 1:
373 mirror = np.conj(c[-1:0:-1]) # conj(c[n-1]), ..., conj(c[1])
374 real_c[1:] = (c[1:] + mirror) / 2
375 return self.__class__(real_c, self._interval)
377 def simplify(self) -> "Trigtech":
378 """Truncate high-frequency Fourier coefficients that are below tolerance.
380 Uses the same one-sided symmetric-maximum criterion as the adaptive
381 constructor: the highest-frequency mode retained is the one where
382 ``max(|c_k|, |c_{-k}|) / vscale > tol``.
383 """
384 oldlen = len(self._coeffs)
385 longself = self.prolong(max(17, oldlen))
386 n = longself.size
387 tol = prefs.eps * max(self._interval.hscale, 1)
389 centered = np.fft.fftshift(longself._coeffs)
390 dc_idx = n // 2
391 abs_max = float(np.max(np.abs(centered)))
392 if abs_max == 0.0:
393 return self.initconst(0.0, interval=self._interval)
395 abs_sym = np.zeros(dc_idx + 1)
396 for ki in range(dc_idx + 1):
397 p = centered[dc_idx + ki] if dc_idx + ki < n else 0.0
398 q = centered[dc_idx - ki]
399 abs_sym[ki] = max(abs(p), abs(q)) / abs_max
401 above = np.where(abs_sym > tol)[0]
402 if len(above) == 0: # pragma: no cover
403 # Defensive: with abs_max > 0 the normalised peak equals 1 > tol,
404 # so 'above' always contains at least the peak index.
405 return self.initconst(0.0, interval=self._interval)
406 max_k = int(above[-1])
407 max_k = min(max_k, oldlen // 2) # don't exceed original size
409 start = dc_idx - max_k
410 end = dc_idx + max_k + 1
411 return self.__class__(np.fft.ifftshift(centered[start:end]), interval=self._interval)
413 def values(self) -> np.ndarray:
414 """Function values at the n equispaced points x_j = -1 + 2j/n."""
415 return self._coeffs2vals(self._coeffs)
417 def _coeffs_to_plotorder(self) -> np.ndarray:
418 """Return coefficients in DC-centred (human-readable) order.
420 Equivalent to ``numpy.fft.fftshift(self.coeffs)``:
421 ordering is [c_{-n//2}, …, c_{-1}, c_0, c_1, …, c_{n//2-1}].
422 """
423 return np.fft.fftshift(self._coeffs)
425 # ------------------------------------------------------------------
426 # algebra
427 # ------------------------------------------------------------------
429 @self_empty()
430 def __add__(self, f: Any) -> "Trigtech":
431 """Add a scalar or another Trigtech."""
432 cls = self.__class__
433 if np.isscalar(f):
434 dtype: Any = complex if np.iscomplexobj(f) else self._coeffs.dtype
435 cfs = np.array(self._coeffs, dtype=dtype)
436 cfs[0] += f # add to DC component
437 return cls(cfs, interval=self._interval)
438 if f.isempty:
439 return cast("Trigtech", f.copy())
440 g = self
441 n, m = g.size, f.size
442 if n < m:
443 g = g.prolong(m)
444 elif m < n:
445 f = f.prolong(n)
446 cfs = f.coeffs + g.coeffs
447 eps = prefs.eps
448 tol = 0.5 * eps * max(f.vscale, g.vscale)
449 if np.all(np.abs(cfs) < tol):
450 return cls.initconst(0.0, interval=self._interval)
451 return cls(cfs, interval=self._interval)
453 @self_empty()
454 def __div__(self, f: Any) -> "Trigtech":
455 """Divide this Trigtech by a scalar or another Trigtech."""
456 cls = self.__class__
457 if np.isscalar(f):
458 return cls(self._coeffs / np.asarray(f), interval=self._interval)
459 if f.isempty:
460 return cast("Trigtech", f.copy())
461 return cls.initfun_adaptive(lambda x: self(x) / f(x), interval=self._interval)
463 __truediv__ = __div__
465 @self_empty()
466 def __mul__(self, g: Any) -> "Trigtech":
467 """Multiply this Trigtech by a scalar or another Trigtech.
469 Trig-polynomial multiplication is circular convolution in frequency
470 space. We implement this cleanly by evaluating both on a grid of
471 size n1 + n2 (sufficient to avoid aliasing), multiplying pointwise,
472 and taking the FFT.
473 """
474 cls = self.__class__
475 if np.isscalar(g):
476 return cls(g * self._coeffs, interval=self._interval)
477 if g.isempty:
478 return cast("Trigtech", g.copy())
479 n = self.size + g.size
480 f_vals = self.prolong(n).values()
481 g_vals = g.prolong(n).values()
482 return cls(cls._vals2coeffs(f_vals * g_vals), interval=self._interval)
484 def __neg__(self) -> "Trigtech":
485 """Return the negation."""
486 return self.__class__(-self._coeffs, interval=self._interval)
488 def __pos__(self) -> "Trigtech":
489 """Return self (unary plus)."""
490 return self
492 @self_empty()
493 def __pow__(self, f: Any) -> "Trigtech":
494 """Raise this Trigtech to a power *f* (scalar or Trigtech)."""
496 def powfun(fn: Any, x: Any) -> Any:
497 return fn if np.isscalar(fn) else fn(x)
499 return self.__class__.initfun_adaptive(
500 lambda x: np.power(self(x), powfun(f, x)),
501 interval=self._interval,
502 )
504 def __rdiv__(self, f: Any) -> "Trigtech":
505 """Compute f / self where *f* is a scalar."""
506 return self.__class__.initfun_adaptive(
507 lambda x: (0.0 * x + f) / self(x),
508 interval=self._interval,
509 )
511 __radd__ = __add__
512 __rmul__ = __mul__
513 __rtruediv__ = __rdiv__
515 def __rsub__(self, f: Any) -> "Trigtech":
516 """Compute f - self."""
517 return cast("Trigtech", -(self - f))
519 @self_empty()
520 def __rpow__(self, f: Any) -> "Trigtech":
521 """Compute f ** self."""
522 return self.__class__.initfun_adaptive(
523 lambda x: np.power(f, self(x)),
524 interval=self._interval,
525 )
527 def __sub__(self, f: Any) -> "Trigtech":
528 """Subtract *f* (scalar or Trigtech) from this Trigtech."""
529 return cast("Trigtech", self + (-f))
531 # ------------------------------------------------------------------
532 # rootfinding
533 # ------------------------------------------------------------------
535 def roots(self, sort: bool | None = None) -> np.ndarray:
536 """Find the roots of this Trigtech on [-1, 1].
538 Converts to a Chebyshev representation via re-sampling on Chebyshev
539 points and delegates to the Chebtech colleague-matrix root-finder.
541 Args:
542 sort: If True, sort the roots in ascending order. Defaults to
543 ``prefs.sortroots``.
544 """
545 from .algorithms import newtonroots, rootsunit
546 from .chebtech import Chebtech
548 sort = sort if sort is not None else prefs.sortroots
550 if self.isempty:
551 return np.array([])
553 # Sample on a Chebyshev grid and fit a Chebtech of the same resolution
554 n = max(2 * self.size + 1, 33)
555 cheb_pts = Chebtech._chebpts(n)
556 vals = self(cheb_pts)
557 ct = Chebtech(Chebtech._vals2coeffs(vals))
558 rts = rootsunit(ct.coeffs)
559 rts = newtonroots(ct, rts)
560 rts = np.clip(rts, -1.0, 1.0)
561 return np.sort(rts) if sort else rts
563 # ------------------------------------------------------------------
564 # calculus
565 # ------------------------------------------------------------------
567 @self_empty(resultif=0.0)
568 def sum(self) -> Any:
569 """Definite integral of the Trigtech over [-1, 1].
571 Only the DC coefficient contributes:
572 ∫_{-1}^{1} exp(i*π*k*(x+1)) dx = 0 for k ≠ 0
573 ∫_{-1}^{1} 1 dx = 2 for k = 0
574 """
575 return 2.0 * float(np.real(self._coeffs[0]))
577 @self_empty()
578 def cumsum(self) -> "Trigtech":
579 """Indefinite integral, zero at x = -1, in Fourier coefficient space.
581 For mode k ≠ 0: antiderivative coefficient = c_k / (i*π*ω_k)
582 For mode k = 0: set to the constant needed so that F(-1) = 0.
584 Note: if the DC component (self.coeffs[0]) is non-zero the true
585 antiderivative contains a linear trend and is not periodic. We still
586 return a Trigtech representing the *periodic* part, adjusted so that
587 the result evaluates to 0 at x = -1.
588 """
589 n = self.size
590 c = self._coeffs.copy()
591 freqs = np.fft.fftfreq(n) * n # FFT-order integer frequencies
593 int_c = np.zeros(n, dtype=complex)
594 mask = freqs != 0
595 int_c[mask] = c[mask] / (1j * np.pi * freqs[mask])
597 # Enforce F(-1) = 0.
598 # F(x) = Σ_k int_c[k] * exp(i*π*ω_k*(x+1))
599 # At x = -1: exp(i*π*ω_k*0) = 1 for all k, so F(-1) = Σ int_c
600 # Set int_c[0] so that sum(int_c) = 0.
601 int_c[0] = -np.sum(int_c[1:])
602 return self.__class__(int_c, interval=self._interval)
604 @self_empty()
605 def diff(self) -> "Trigtech":
606 """Derivative via the Fourier multiplier i*π*ω_k.
608 d/dx [c_k * exp(i*π*ω_k*(x+1))] = i*π*ω_k * c_k * exp(i*π*ω_k*(x+1))
609 """
610 if self.isconst:
611 return self.__class__(np.array([0.0 + 0.0j]), interval=self._interval)
612 n = self.size
613 freqs = np.fft.fftfreq(n) * n
614 d_coeffs = (1j * np.pi * freqs) * self._coeffs
615 return self.__class__(d_coeffs, interval=self._interval)
617 # ------------------------------------------------------------------
618 # static helpers (FFT ↔ values)
619 # ------------------------------------------------------------------
621 @staticmethod
622 def _trigpts(n: int) -> np.ndarray:
623 """Return *n* equispaced points on [-1, 1)."""
624 if n == 0:
625 return np.array([])
626 return -1.0 + 2.0 * np.arange(n) / n
628 @staticmethod
629 def _vals2coeffs(vals: Any) -> np.ndarray:
630 """Convert values at equispaced points to FFT coefficients (divided by n).
632 Always returns complex128, even for real-valued inputs, because Fourier
633 coefficients for functions such as sin are purely imaginary and would be
634 discarded if forced to real.
636 Inverse of ``_coeffs2vals``.
637 """
638 vals = np.asarray(vals)
639 n = vals.size
640 if n == 0:
641 return np.array([], dtype=complex)
642 return cast(np.ndarray, np.fft.fft(vals) / n)
644 @staticmethod
645 def _coeffs2vals(coeffs: Any) -> np.ndarray:
646 """Convert FFT coefficients (divided by n) to values at equispaced points.
648 Inverse of ``_vals2coeffs``.
649 """
650 coeffs = np.asarray(coeffs, dtype=complex)
651 n = coeffs.size
652 if n == 0:
653 return np.array([], dtype=float)
654 vals = n * np.fft.ifft(coeffs)
655 # Discard negligible imaginary parts for conjugate-symmetric coefficients
656 max_real = float(np.max(np.abs(np.real(vals))))
657 if float(np.max(np.abs(np.imag(vals)))) < 1e-10 * max(max_real, 1.0):
658 return np.real(vals)
659 return cast(np.ndarray, vals)
661 # ------------------------------------------------------------------
662 # plotting
663 # ------------------------------------------------------------------
665 def plot(self, ax: Any = None, **kwargs: Any) -> Any:
666 """Plot the Trigtech over [-1, 1].
668 Args:
669 ax: Matplotlib axes. If None, uses the current axes.
670 **kwargs: Forwarded to matplotlib.
672 Returns:
673 The axes on which the plot was drawn.
674 """
675 return plotfun(self, (-1, 1), ax=ax, **kwargs)
677 def plotcoeffs(self, ax: Any = None, **kwargs: Any) -> Any:
678 """Plot the absolute Fourier coefficient magnitudes in DC-centred order.
680 Uses ``_coeffs_to_plotorder()`` so the horizontal axis runs from
681 the most-negative frequency on the left to the most-positive on
682 the right, with DC in the centre.
684 Args:
685 ax: Matplotlib axes. If None, uses the current axes.
686 **kwargs: Forwarded to matplotlib.
688 Returns:
689 The axes on which the plot was drawn.
690 """
691 ax = ax or plt.gca()
692 return plotfuncoeffs(np.abs(self._coeffs_to_plotorder()), ax=ax, **kwargs)