pyhctsa.operations.spectral.spectral_summaries

pyhctsa.operations.spectral.spectral_summaries(y, psd_meth='fft', window_type='hamming')

Statistics of the power spectrum of a time series.

Estimates the power spectrum (as a periodogram, a plain fast Fourier transform, or by Welch’s method) and returns many summary statistics: the location and width of its main peak, the number and prominence of peaks, the distribution of power values, the autocorrelation of power across frequency, the frequencies below which given fractions of the power lie, power-weighted moments of frequency, fits to the cumulative power, a spectral entropy and flatness, robust power-law fits to the log-log spectrum, the power in 2 and 5 equal frequency bands, and the number of crossings of the log spectrum at various levels. Many statistics have a log-domain version computed on log(S). The spectrum is floored at 1e-12 of its maximum before any statistic is computed, so that bins at the rounding level of the estimator do not determine the log-domain statistics.

Parameters:
y : array-like

The input time series.

psd_meth : {'periodogram', 'fft', 'welch'}, optional

The method for obtaining the spectrum from the signal:

  • ’periodogram’: periodogram (window across the whole series)

  • ’fft’: fast Fourier transform, single-sided power spectral density; no window is applied and the DC bin is dropped

  • ’welch’: Welch’s method, with windows of max(min(256, round(N/4)), 16) samples, 50% overlap, no detrending, and each segment zero-padded to nfft = max(256, next power of 2 >= window length) (also used by ‘periodogram’, with the window length N)

Default is 'fft'.

window_type : {'boxcar', 'rect', 'bartlett', 'hann', 'hamming', 'none'}, optional

The window to use (for ‘periodogram’ and ‘welch’; ignored by ‘fft’). Default is 'hamming'.

Returns:

The spectrum S is a power spectral density in angular frequency w (radians per sample, 0 to pi), normalized for all three estimators so that its area (the sum of S times the bin spacing dw) equals the variance of the series (~1 for a z-scored series). Statistics that accumulate over bins (cumulative-area fits, entropy) are therefore integrals over w, independent of the number of bins. Returns NaN for a constant series. Output fields:

Peaks

  • maxS, maxw: the maximum of S and the angular frequency at which it occurs.

  • maxWidth: half-power bandwidth of the dominant peak.

  • numPromPeaks_3, numPromPeaks_5, numPromPeaks_8: number of peaks of log(S) with prominence above 3, 5, 8 (natural-log units).

  • meanProm_5, meanPeakWidth_prom5: mean prominence and mean width of the peaks with prominence above 5.

  • width_weighted_prom, w_weighted_peak_prom: mean peak width and mean peak location, weighted by prominence.

  • peakPower_2, peakPower_5, peakPower_prom5: power (height x width) in the 2 and 5 tallest peaks and in the peaks with prominence above 5.

  • numPeaks_50power, peakpower_1: number of tallest peaks needed to hold half the peak power, and the fraction of peak power in the tallest peak.

Peaks are found on log(S), and the prominence thresholds are calibrated for (and hctsa registers these fields only with) psd_meth='welch'.

Distribution of power values (and of log power values, with prefix log)

iqr, logiqr, q25, median, q75, logq25, logmedian, logq75, std, stdlog (log of the std), logstd (std of log(S)), mom3 and logmom3 (skewness).

Autocorrelation of the spectrum across frequency

ac1, ac2, tau (first zero-crossing of the autocorrelation, in units of w), and logac1, logac2, logtau for log(S).

Cumulative power

wq_5, wq_10, wq_25, wq_50, wq_75, wq_90, wq_95, wq_99: the frequency below which 5%, 10%, …, 99% of the power lies (wq_50 is the median frequency).

Power-weighted moments of frequency

specCentroid, specSpread, specSkew, specKurt: mean, standard deviation, skewness and kurtosis of frequency weighted by power.

Fits to the cumulative area under S (a running integral over w)

fpoly2csS_p1, fpoly2csS_p2, fpoly2csS_p3, fpoly2_sse, fpoly2_r2, fpoly2_rmse: coefficients and goodness of a quadratic fit (fpoly2_sse is the integrated squared error); fpolysat_a, fpolysat_b, fpolysat_r2, fpolysat_rmse: parameters and goodness of the least-squares fit of a*w**2/(b+w**2), the global optimum over 1e-6 <= b/max(w)**2 <= 100 (see _polysat_fit; a curve that keeps rising like w**2 gives the largest b).

Entropy, flatness and areas

spect_shann_ent (differential Shannon entropy of the power distribution over frequency), spect_shann_ent_norm (exp(spect_shann_ent) as a fraction of the frequency range), sfm (spectral flatness measure, 10*log10(geometric mean / arithmetic mean)), areatopeak and ylogareatopeak (area under S, and under log(S), up to the peak).

Robust linear fits (a1 intercept, a2 gradient, sigrat OLS/robust sigma ratio, sigma, sea1 standard error of the intercept)

linfitloglog_all_{a1,a2,sigrat,sigma,sea1} (log(S) against log(w), all frequencies), linfitloglog_lf_a2 and linfitloglog_mf_a2 (lower half and middle half), linfitloglog_hf_{a1,a2,sigrat,sigma,sea1} (upper half) and linfitsemilog_all_{a1,sigrat,sigma,sea1} (log(S) against w).

Power in frequency bands (2 and 5 equal bands, from the lowest frequencies)

area_2_1, area_2_2, logarea_2_1, logarea_2_2 and area_5_1 … area_5_5, logarea_5_1 … logarea_5_5 (area under S and log(S) in each band); statav2_s, statav5_s (std across bands of the within-band std of S, relative to std(S)); logstatav2_m, logstatav2_s, logstatav5_m, logstatav5_s (the same for the band means and stds of log(S), relative to std(log(S))). When the number of bins is not divisible by the number of bands the last few bins are dropped.

Crossings of the log spectrum

ncross_log_f05, ncross_log_f10, ncross_log_f20, ncross_log_f50: crossings of a level 5%, 10%, 20%, 50% of the way from min(log S) to max(log S).

Return type:

dict