Skip to content

Spectral utils

Utility functions for spectral analysis. Defaults are chosen to replicate MATLAB's default behavior.

csd 

csd(x, y, fs, num_windows=8, window_type='hamming', window_len=None, nfft=None, detrend=False, onesided=True)

Cross spectral density via Welch's method.

Parameters:

Name Type Description Default
x ndarray

Two input signals of identical shape; the longest axis is time.

 required
y ndarray

Two input signals of identical shape; the longest axis is time.

 required
fs float

Sampling frequency (Hz).

 required
num_windows int

Number of (50%-overlapping) Welch windows when window_len is not given.

 8
window_type str

Window passed to scipy.signal.csd.

 'hamming'
window_len int

Window length in samples. If None, derived from num_windows and N.

 None
nfft int

FFT length. Defaults to window_len.

 None
detrend bool

If True, detrend each segment before transforming.

 False
onesided bool

If True, return the one-sided spectrum.

 True

Returns:

Name Type Description
f ndarray

Frequency vector (Hz).
Pxy ndarray

Complex cross spectral density (units of x*y / Hz).
Source code in src/pytoast/utils/spectral_utils.py 
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
def csd(
    x: np.ndarray,
    y: np.ndarray,
    fs: float,
    num_windows: int = 8,
    window_type: str = "hamming",
    window_len: int | None = None,
    nfft: int | None = None,
    detrend: bool = False,
    onesided: bool = True,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Cross spectral density via Welch's method.

    Parameters
    ----------
    x, y : np.ndarray
        Two input signals of identical shape; the longest axis is time.
    fs : float
        Sampling frequency (Hz).
    num_windows : int, optional
        Number of (50%-overlapping) Welch windows when ``window_len`` is not given.
    window_type : str, optional
        Window passed to ``scipy.signal.csd``.
    window_len : int, optional
        Window length in samples. If None, derived from ``num_windows`` and ``N``.
    nfft : int, optional
        FFT length. Defaults to ``window_len``.
    detrend : bool, optional
        If True, detrend each segment before transforming.
    onesided : bool, optional
        If True, return the one-sided spectrum.

    Returns
    -------
    f : np.ndarray
        Frequency vector (Hz).
    Pxy : np.ndarray
        Complex cross spectral density (units of x*y / Hz).
    """
    N = max(x.shape)
    if window_len is None:
        window_len = get_window_len(N, num_windows)
    if nfft is None:
        nfft = window_len

    f, Pxy = sig.csd(
        x=x,
        y=y,
        fs=fs,
        window=window_type,
        nperseg=window_len,
        nfft=nfft,
        detrend=detrend,
        return_onesided=onesided,
    )

    return f, Pxy

get_frequency_range 

get_frequency_range(f, f_low=None, f_high=None)

Index range into f covering [f_low, f_high].

Parameters:

Name Type Description Default
f ndarray

Monotonically increasing frequency vector (Hz).

 required
f_low float

Lower frequency bound (Hz). If None, start at index 0.

 None
f_high float

Upper frequency bound (Hz). If None, end at len(f).

 None

Returns:

Type Description
tuple of int

(start_index, end_index) into f.
Source code in src/pytoast/utils/spectral_utils.py 
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
def get_frequency_range(
    f: np.ndarray, f_low: float | None = None, f_high: float | None = None
) -> tuple[int, int]:
    """
    Index range into ``f`` covering [f_low, f_high].

    Parameters
    ----------
    f : np.ndarray
        Monotonically increasing frequency vector (Hz).
    f_low : float, optional
        Lower frequency bound (Hz). If None, start at index 0.
    f_high : float, optional
        Upper frequency bound (Hz). If None, end at ``len(f)``.

    Returns
    -------
    tuple of int
        ``(start_index, end_index)`` into ``f``.
    """
    if f_low is not None:
        start_index = int(np.argmin(np.abs(f - f_low)))
    else:
        start_index = 0

    if f_high is not None:
        end_index = int(np.argmin(np.abs(f - f_high)))
    else:
        end_index = len(f)

    return start_index, end_index

get_window_len 

get_window_len(N, num_windows)

Welch-method window length.

Parameters:

Name Type Description Default
N int

Number of samples in the time series.

 required
num_windows int

Number of (50%-overlapping) windows desired.

 required

Returns:

Type Description
int

Window length in samples.
Source code in src/pytoast/utils/spectral_utils.py 
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
def get_window_len(N: int, num_windows: int) -> int:
    """
    Welch-method window length.

    Parameters
    ----------
    N : int
        Number of samples in the time series.
    num_windows : int
        Number of (50%-overlapping) windows desired.

    Returns
    -------
    int
        Window length in samples.
    """
    return int(2 * N / (num_windows + 1))

psd 

psd(x, fs, num_windows=8, window_type='hamming', window_len=None, nfft=None, detrend=False, onesided=True)

Power spectral density via Welch's method.

Parameters:

Name Type Description Default
x ndarray

Input signal. The longest axis is treated as time.

 required
fs float

Sampling frequency (Hz).

 required
num_windows int

Number of (50%-overlapping) Welch windows when window_len is not given.

 8
window_type str

Window passed to scipy.signal.welch (default 'hamming').

 'hamming'
window_len int

Window length in samples. If None, derived from num_windows and N.

 None
nfft int

FFT length. Defaults to window_len.

 None
detrend bool

If True, detrend each segment before transforming.

 False
onesided bool

If True, return the one-sided spectrum.

 True

Returns:

Name Type Description
f ndarray

Frequency vector (Hz).
Pxx ndarray

Power spectral density (units of x^2 / Hz).
Source code in src/pytoast/utils/spectral_utils.py 
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
def psd(
    x: np.ndarray,
    fs: float,
    num_windows: int = 8,
    window_type: str = "hamming",
    window_len: int | None = None,
    nfft: int | None = None,
    detrend: bool = False,
    onesided: bool = True,
) -> tuple[np.ndarray, np.ndarray]:
    """
    Power spectral density via Welch's method.

    Parameters
    ----------
    x : np.ndarray
        Input signal. The longest axis is treated as time.
    fs : float
        Sampling frequency (Hz).
    num_windows : int, optional
        Number of (50%-overlapping) Welch windows when ``window_len`` is not given.
    window_type : str, optional
        Window passed to ``scipy.signal.welch`` (default ``'hamming'``).
    window_len : int, optional
        Window length in samples. If None, derived from ``num_windows`` and ``N``.
    nfft : int, optional
        FFT length. Defaults to ``window_len``.
    detrend : bool, optional
        If True, detrend each segment before transforming.
    onesided : bool, optional
        If True, return the one-sided spectrum.

    Returns
    -------
    f : np.ndarray
        Frequency vector (Hz).
    Pxx : np.ndarray
        Power spectral density (units of ``x``^2 / Hz).
    """
    N = max(x.shape)
    if window_len is None:
        window_len = get_window_len(N, num_windows)
    if nfft is None:
        nfft = window_len

    f, Pxx = sig.welch(
        x=x,
        fs=fs,
        window=window_type,
        nperseg=window_len,
        nfft=nfft,
        detrend=detrend,
        return_onesided=onesided,
    )

    return f, Pxx