25194239086

Committed 26 Apr 2026 08:07PM UTC coverage: 61.697% (-28.8%) from 90.54%

Build # 25194239086

Build Type

push

github

Committed by

jameswilburlewis

Commit Message

Added test for loading Cluster CODIF differential energy flux

Coverage Stats

0 of 7 new or added lines in 1 file covered. (0.0%)

19460 existing lines in 418 files now uncovered.

30204 of 48955 relevant lines covered (61.7%)

1.44 hits per line

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

8.0

/pyspedas/tplot_tools/tplot_math/pwrspc.py

import logging
import numpy as np
from scipy.stats import linregress


def pwrspc(time, quantity, noline=False, nohanning=False, bin=3, notperhz=False):
    """
    Compute the power spectrum of a given time series.

    Parameters
    ----------
        time (array):
            The time array.
        quantity (array):
            The data array for which the power spectrum is to be computed.
            Should be one dimensional and the same length as time.
        noline (bool):
            If True, straight line is not subtracted from the data.
        nohanning (bool):
            If True, no Hanning window is applied to the data.
        bin (int):
            Bin size for binning the data. Default is 3.
        notperhz (bool):
            If True, the output units are the square of the input units.

    Returns
    -------
    tuple
        Tuple containing:
            - freq (array):
                The frequency array.
            - power (array):
                The computed power spectrum.

    Notes:
        This is similar to IDL pwrspc.pro routine.

        A Hanning window is applied to the input data, and its power is divided out of the returned spectrum.
        A straight line is subtracted from the data to reduce spurious power due to sawtooth behavior of a background.
        Units are (units)^2 where units are the units of the input quantity. Frequency is in 1/time units.
        Thus, the output represents the mean squared amplitude of the signal at each specific frequency.
        The total (sum) power under the curve is equal to the mean (over time) power of the oscillation in the time domain.
        If the keyword notperhz is True, then power is in units^2. If notperhz is False (default), power is in units^2/Hz.


    Examples
    --------

        >>> # Compute the power spectrum of a given time series
        >>> from pyspedas import pwrspc
        >>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
        >>> quantity = [1,2,3,1,2,3,1,2,3,1]
        >>> freq, power = pwrspc(time, quantity)
        >>> print(freq, power)
    """

    t = np.array(time, dtype=np.float64)
    x = np.array(quantity, dtype=np.float64)

    # If the dimensions of the input arrays are not the same, and not one dimension, return
    if t.ndim != 1 or x.ndim != 1 or len(t) != len(x) or len(t) < 1:
        logging.error(
            "Both input arrays should be one dimensional and of the same length."
        )
        return np.array(None), np.array(None)

    # Subtract first point from time array
    t -= t[0]

    # Subtract straight line from data
    if not noline:
        slope, intercept, _, _, _ = linregress(t, x)
        x -= slope * t + intercept

    binsize = bin
    window = 0.0
    if not nohanning:
        window = np.hanning(len(x))
        x *= window

    nt = len(t)
    if nt % 2 != 0:
        logging.info("Needs an even number of data points, dropping last point...")
        t = t[:-1]
        x = x[:-1]
        nt -= 1

    xs2 = np.abs(np.fft.fft(x)) ** 2
    dbign = float(nt)
    logging.info("bign=" + str(dbign))

    k = np.arange(0, dbign // 2 + 1)
    tres = float(np.median(np.diff(t)))
    fk = k / (dbign * tres)

    pwr = np.zeros(nt // 2 + 1)
    pwr[0] = xs2[0] / dbign**2
    pwr[1 : nt // 2] = (xs2[1 : nt // 2] + xs2[nt : nt // 2 : -1]) / dbign**2
    pwr[-1] = xs2[-1] / dbign**2

    if not nohanning:
        wss = dbign * np.sum(window**2)
        pwr = pwr * dbign**2 / wss

    dfreq = binsize * (fk[1] - fk[0])
    npwr = len(pwr) - 1
    nfinal = int(npwr / binsize)
    iarray = np.arange(nfinal)
    power = np.zeros(nfinal)

    idx = (iarray + 0.5) * binsize + 1
    freq = [fk[int(i)] for i in idx]

    for i in range(binsize):
        power += pwr[iarray * binsize + i + 1]

    if not notperhz:
        power /= dfreq

    logging.info("dfreq=" + str(dfreq))

    return np.array(freq), np.array(power)

1	import logging	3✔
2	import numpy as np	3✔
3	from scipy.stats import linregress	3✔
4
5
6	def pwrspc(time, quantity, noline=False, nohanning=False, bin=3, notperhz=False):	3✔
7	"""
8	Compute the power spectrum of a given time series.
9
10	Parameters
11	----------
12	time (array):
13	The time array.
14	quantity (array):
15	The data array for which the power spectrum is to be computed.
16	Should be one dimensional and the same length as time.
17	noline (bool):
18	If True, straight line is not subtracted from the data.
19	nohanning (bool):
20	If True, no Hanning window is applied to the data.
21	bin (int):
22	Bin size for binning the data. Default is 3.
23	notperhz (bool):
24	If True, the output units are the square of the input units.
25
26	Returns
27	-------
28	tuple
29	Tuple containing:
30	- freq (array):
31	The frequency array.
32	- power (array):
33	The computed power spectrum.
34
35	Notes:
36	This is similar to IDL pwrspc.pro routine.
37
38	A Hanning window is applied to the input data, and its power is divided out of the returned spectrum.
39	A straight line is subtracted from the data to reduce spurious power due to sawtooth behavior of a background.
40	Units are (units)^2 where units are the units of the input quantity. Frequency is in 1/time units.
41	Thus, the output represents the mean squared amplitude of the signal at each specific frequency.
42	The total (sum) power under the curve is equal to the mean (over time) power of the oscillation in the time domain.
43	If the keyword notperhz is True, then power is in units^2. If notperhz is False (default), power is in units^2/Hz.
44
45
46	Examples
47	--------
48
49	>>> # Compute the power spectrum of a given time series
50	>>> from pyspedas import pwrspc
51	>>> time = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
52	>>> quantity = [1,2,3,1,2,3,1,2,3,1]
53	>>> freq, power = pwrspc(time, quantity)
54	>>> print(freq, power)
55	"""
56
UNCOV 57	t = np.array(time, dtype=np.float64)	×
UNCOV 58	x = np.array(quantity, dtype=np.float64)	×
59
60	# If the dimensions of the input arrays are not the same, and not one dimension, return
UNCOV 61	if t.ndim != 1 or x.ndim != 1 or len(t) != len(x) or len(t) < 1:	×
62	logging.error(	×
63	"Both input arrays should be one dimensional and of the same length."
64	)
65	return np.array(None), np.array(None)	×
66
67	# Subtract first point from time array
UNCOV 68	t -= t[0]	×
69
70	# Subtract straight line from data
UNCOV 71	if not noline:	×
UNCOV 72	slope, intercept, _, _, _ = linregress(t, x)	×
UNCOV 73	x -= slope * t + intercept	×
74
UNCOV 75	binsize = bin	×
UNCOV 76	window = 0.0	×
UNCOV 77	if not nohanning:	×
UNCOV 78	window = np.hanning(len(x))	×
UNCOV 79	x *= window	×
80
UNCOV 81	nt = len(t)	×
UNCOV 82	if nt % 2 != 0:	×
83	logging.info("Needs an even number of data points, dropping last point...")	×
84	t = t[:-1]	×
85	x = x[:-1]	×
86	nt -= 1	×
87
UNCOV 88	xs2 = np.abs(np.fft.fft(x)) ** 2	×
UNCOV 89	dbign = float(nt)	×
UNCOV 90	logging.info("bign=" + str(dbign))	×
91
UNCOV 92	k = np.arange(0, dbign // 2 + 1)	×
UNCOV 93	tres = float(np.median(np.diff(t)))	×
UNCOV 94	fk = k / (dbign * tres)	×
95
UNCOV 96	pwr = np.zeros(nt // 2 + 1)	×
UNCOV 97	pwr[0] = xs2[0] / dbign**2	×
UNCOV 98	pwr[1 : nt // 2] = (xs2[1 : nt // 2] + xs2[nt : nt // 2 : -1]) / dbign**2	×
UNCOV 99	pwr[-1] = xs2[-1] / dbign**2	×
100
UNCOV 101	if not nohanning:	×
UNCOV 102	wss = dbign * np.sum(window**2)	×
UNCOV 103	pwr = pwr * dbign**2 / wss	×
104
UNCOV 105	dfreq = binsize * (fk[1] - fk[0])	×
UNCOV 106	npwr = len(pwr) - 1	×
UNCOV 107	nfinal = int(npwr / binsize)	×
UNCOV 108	iarray = np.arange(nfinal)	×
UNCOV 109	power = np.zeros(nfinal)	×
110
UNCOV 111	idx = (iarray + 0.5) * binsize + 1	×
UNCOV 112	freq = [fk[int(i)] for i in idx]	×
113
UNCOV 114	for i in range(binsize):	×
UNCOV 115	power += pwr[iarray * binsize + i + 1]	×
116
UNCOV 117	if not notperhz:	×
UNCOV 118	power /= dfreq	×
119
UNCOV 120	logging.info("dfreq=" + str(dfreq))	×
121
UNCOV 122	return np.array(freq), np.array(power)	×

spedas / pyspedas / 25194239086

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