21430410834

Committed 28 Jan 2026 08:14AM UTC coverage: 75.545% (+6.5%) from 69.075%

Build # 21430410834

Build Type

push

github

Committed by

avivajpeyi

Commit Message

fix: update offs calculation in transform_wavelet_freq_helper to prevent overwriting center term

Run Details

192 of 308 branches covered (62.34%)

Branch coverage included in aggregate %.

3 of 3 new or added lines in 1 file covered. (100.0%)

164 existing lines in 16 files now uncovered.

1229 of 1573 relevant lines covered (78.13%)

0.78 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

96.15

/src/pywavelet/transforms/phi_computer.py

"""
This module contains functions to compute the Fourier transform of the
wavelet function and its normalization. The wavelet function is defined
in the frequency domain and is used to transform time-domain data into
the wavelet domain.

Everything in this module is retured as a npfloat64 array.
"""

from functools import lru_cache

import numpy as np
from jaxtyping import Float64
from numpy.fft import ifft
from scipy.special import betainc

__all__ = ["phitilde_vec_norm", "phi_vec", "omega"]


@lru_cache(maxsize=None)
def omega(Nf: int, Nt: int) -> Float64[np.ndarray, "{Nt}//2+1"]:
    """Get the angular frequencies of the time domain signal."""
    df = 2 * np.pi / (Nf * Nt)
    return df * np.arange(0, Nt // 2 + 1, dtype=np.float64)


@lru_cache(maxsize=None)
def phitilde_vec_norm(
    Nf: int, Nt: int, d: float
) -> Float64[np.ndarray, "{Nt}//2+1"]:
    """Normalize phitilde for inverse frequency domain transform."""
    omegas = omega(Nf, Nt)
    _phi_t = _phitilde_vec(omegas, Nf, d) * np.sqrt(np.pi)
    return np.array(_phi_t)


@lru_cache(maxsize=None)
def phi_vec(
    Nf: int, d: float = 4.0, q: int = 16
) -> Float64[np.ndarray, "2*{q}*{Nf}"]:
    """get time domain phi as fourier transform of _phitilde_vec
    q: number of Nf bins over which the window extends?

    """
    insDOM = 1.0 / np.sqrt(np.pi / Nf)
    K = q * 2 * Nf
    half_K = q * Nf  # xp.int64(K/2)

    dom = 2 * np.pi / K  # max frequency is K/2*dom = pi/dt = OM
    DX = np.zeros(K, dtype=np.complex128)

    # zero frequency
    DX[0] = insDOM

    DX = DX.copy()
    # postive frequencies
    DX[1 : half_K + 1] = _phitilde_vec(dom * np.arange(1, half_K + 1), Nf, d)
    # negative frequencies
    DX[half_K + 1 :] = _phitilde_vec(
        -dom * np.arange(half_K - 1, 0, -1), Nf, d
    )
    DX = K * ifft(DX, K)

    phi = np.zeros(K)
    phi[0:half_K] = np.real(DX[half_K:K])
    phi[half_K:] = np.real(DX[0:half_K])

    nrm = np.sqrt(2.0) / np.sqrt(K / dom)  # *xp.linalg.norm(phi)

    phi *= nrm
    return np.array(phi)


def _phitilde_vec(
    omega: Float64[np.ndarray, "dim"], Nf: int, d: float = 4.0
) -> Float64[np.ndarray, "dim"]:
    """Compute phi_tilde(omega_i) array, nx is filter steepness, defaults to 4.

    Eq 11 of https://arxiv.org/pdf/2009.00043.pdf (Cornish et al. 2020)

    phi(omega_i) =
        1/sqrt(2π∆F) if |omega_i| < A
        1/sqrt(2π∆F) cos(nu_d π/2 * |omega|-A / B) if A < |omega_i| < A + B

    Where nu_d = normalized incomplete beta function

    Parameters
    ----------
    ω : xp.ndarray
        Array of angular frequencies
    Nf : int
        Number of frequency bins
    d : float, optional
        Number of standard deviations for the gaussian wavelet, by default 4.

    Returns
    -------
    xp.ndarray
        Array of phi_tilde(omega_i) values

    """
    dF = 1.0 / (2 * Nf)  # NOTE: missing 1/dt?
    dOmega = 2 * np.pi * dF  # Near Eq 10 # 2 pi times DF
    inverse_sqrt_dOmega = 1.0 / np.sqrt(dOmega)

    A = dOmega / 4
    B = dOmega - 2 * A  # Cannot have B \leq 0.
    if B <= 0:
        raise ValueError("B must be greater than 0")

    phi = np.zeros(omega.size, dtype=np.float64)
    mask = (A <= np.abs(omega)) & (np.abs(omega) < A + B)  # Minor changes
    vd = (np.pi / 2.0) * _nu_d(omega[mask], A, B, d=d)  # different from paper
    phi[mask] = inverse_sqrt_dOmega * np.cos(vd)
    phi[np.abs(omega) < A] = inverse_sqrt_dOmega
    return phi


def _nu_d(
    omega: Float64[np.ndarray, "dim"], A: float, B: float, d: float = 4.0
) -> Float64[np.ndarray, "dim"]:
    """Compute the normalized incomplete beta function.

    Parameters
    ----------
    omega : np.ndarray
        Array of angular frequencies
    A : float
        Lower bound for the beta function
    B : float
        Upper bound for the beta function
    d : float, optional
        Number of standard deviations for the gaussian wavelet, by default 4.

    Returns
    -------
    np.ndarray
        Array of ν_d values

    scipy.special.betainc
    https://docs.scipy.org/doc/scipy-1.7.1/reference/reference/generated/scipy.special.betainc.html

    """
    x = (np.abs(omega) - A) / B
    return betainc(d, d, x)

1	"""
2	This module contains functions to compute the Fourier transform of the
3	wavelet function and its normalization. The wavelet function is defined
4	in the frequency domain and is used to transform time-domain data into
5	the wavelet domain.
6
7	Everything in this module is retured as a npfloat64 array.
8	"""
9
10	from functools import lru_cache	1✔
11
12	import numpy as np	1✔
13	from jaxtyping import Float64	1✔
14	from numpy.fft import ifft	1✔
15	from scipy.special import betainc	1✔
16
17	__all__ = ["phitilde_vec_norm", "phi_vec", "omega"]	1✔
18
19
20	@lru_cache(maxsize=None)	1✔
21	def omega(Nf: int, Nt: int) -> Float64[np.ndarray, "{Nt}//2+1"]:	1✔
22	"""Get the angular frequencies of the time domain signal."""
23	df = 2 * np.pi / (Nf * Nt)	1✔
24	return df * np.arange(0, Nt // 2 + 1, dtype=np.float64)	1✔
25
26
27	@lru_cache(maxsize=None)	1✔
28	def phitilde_vec_norm(	1✔
29	Nf: int, Nt: int, d: float
30	) -> Float64[np.ndarray, "{Nt}//2+1"]:
31	"""Normalize phitilde for inverse frequency domain transform."""
32	omegas = omega(Nf, Nt)	1✔
33	_phi_t = _phitilde_vec(omegas, Nf, d) * np.sqrt(np.pi)	1✔
34	return np.array(_phi_t)	1✔
35
36
37	@lru_cache(maxsize=None)	1✔
38	def phi_vec(	1✔
39	Nf: int, d: float = 4.0, q: int = 16
40	) -> Float64[np.ndarray, "2{q}{Nf}"]:
41	"""get time domain phi as fourier transform of _phitilde_vec
42	q: number of Nf bins over which the window extends?
43
44	"""
45	insDOM = 1.0 / np.sqrt(np.pi / Nf)	1✔
46	K = q * 2 * Nf	1✔
47	half_K = q * Nf # xp.int64(K/2)	1✔
48
49	dom = 2 * np.pi / K # max frequency is K/2*dom = pi/dt = OM	1✔
50	DX = np.zeros(K, dtype=np.complex128)	1✔
51
52	# zero frequency
53	DX[0] = insDOM	1✔
54
55	DX = DX.copy()	1✔
56	# postive frequencies
57	DX[1 : half_K + 1] = _phitilde_vec(dom * np.arange(1, half_K + 1), Nf, d)	1✔
58	# negative frequencies
59	DX[half_K + 1 :] = _phitilde_vec(	1✔
60	-dom * np.arange(half_K - 1, 0, -1), Nf, d
61	)
62	DX = K * ifft(DX, K)	1✔
63
64	phi = np.zeros(K)	1✔
65	phi[0:half_K] = np.real(DX[half_K:K])	1✔
66	phi[half_K:] = np.real(DX[0:half_K])	1✔
67
68	nrm = np.sqrt(2.0) / np.sqrt(K / dom) # *xp.linalg.norm(phi)	1✔
69
70	phi *= nrm	1✔
71	return np.array(phi)	1✔
72
73
74	def _phitilde_vec(	1✔
75	omega: Float64[np.ndarray, "dim"], Nf: int, d: float = 4.0
76	) -> Float64[np.ndarray, "dim"]:
77	"""Compute phi_tilde(omega_i) array, nx is filter steepness, defaults to 4.
78
79	Eq 11 of https://arxiv.org/pdf/2009.00043.pdf (Cornish et al. 2020)
80
81	phi(omega_i) =
82	1/sqrt(2π∆F) if \|omega_i\| < A
83	1/sqrt(2π∆F) cos(nu_d π/2 * \|omega\|-A / B) if A < \|omega_i\| < A + B
84
85	Where nu_d = normalized incomplete beta function
86
87	Parameters
88	----------
89	ω : xp.ndarray
90	Array of angular frequencies
91	Nf : int
92	Number of frequency bins
93	d : float, optional
94	Number of standard deviations for the gaussian wavelet, by default 4.
95
96	Returns
97	-------
98	xp.ndarray
99	Array of phi_tilde(omega_i) values
100
101	"""
102	dF = 1.0 / (2 * Nf) # NOTE: missing 1/dt?	1✔
103	dOmega = 2 * np.pi * dF # Near Eq 10 # 2 pi times DF	1✔
104	inverse_sqrt_dOmega = 1.0 / np.sqrt(dOmega)	1✔
105
106	A = dOmega / 4	1✔
107	B = dOmega - 2 * A # Cannot have B \leq 0.	1✔
108	if B <= 0:	1!
UNCOV 109	raise ValueError("B must be greater than 0")	×
110
111	phi = np.zeros(omega.size, dtype=np.float64)	1✔
112	mask = (A <= np.abs(omega)) & (np.abs(omega) < A + B) # Minor changes	1✔
113	vd = (np.pi / 2.0) * _nu_d(omega[mask], A, B, d=d) # different from paper	1✔
114	phi[mask] = inverse_sqrt_dOmega * np.cos(vd)	1✔
115	phi[np.abs(omega) < A] = inverse_sqrt_dOmega	1✔
116	return phi	1✔
117
118
119	def _nu_d(	1✔
120	omega: Float64[np.ndarray, "dim"], A: float, B: float, d: float = 4.0
121	) -> Float64[np.ndarray, "dim"]:
122	"""Compute the normalized incomplete beta function.
123
124	Parameters
125	----------
126	omega : np.ndarray
127	Array of angular frequencies
128	A : float
129	Lower bound for the beta function
130	B : float
131	Upper bound for the beta function
132	d : float, optional
133	Number of standard deviations for the gaussian wavelet, by default 4.
134
135	Returns
136	-------
137	np.ndarray
138	Array of ν_d values
139
140	scipy.special.betainc
141	https://docs.scipy.org/doc/scipy-1.7.1/reference/reference/generated/scipy.special.betainc.html
142
143	"""
144	x = (np.abs(omega) - A) / B	1✔
145	return betainc(d, d, x)	1✔

avivajpeyi / pywavelet / 21430410834

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous