PaperlayerView exact revision on arXiv

Source-linked AI summary

Modulation Analysis with Higher-Order Spectra

Christopher K. Kovach, Sukhbinder Kumar

arXiv:2609.03172v1eess.SP

TL;DR

HOS capture cross-frequency dependencies, but their real-world meaning has been perceived as obscure. The paper interprets the trispectrum as a cross spectrum of the Wigner-Ville distribution and reports beta-range modulation patterns in recordings. It also relates time-invariant fourth-order statistics to conditional second-order nonstationarities.

  • Problem

    HOS describe cross-frequency dependencies, but their perceived obscurity of real-world meaning has hampered adoption.

    Source: Introduction · pp. 1–2

  • Method

    The paper treats the Wigner-Ville distribution as a multivariate time series and uses the trispectrum to represent linear dependence of power across frequencies.

    Source: Trispectrum as a measure of linear… · p. 5

  • Results

    The reported modulation band spans 19 to 25 Hz, with amplitude-modulation energy from below 0.5 Hz to about 1.5 Hz, consistent with aperiodic burst-like envelope modulation.

    Source: Data set · p. 15

  • Takeaways & Limitations

    For ergodic conditionally nonstationary systems, time-invariant fourth-order statistics can recover the spectrum of conditional second-order nonstationarities.

    Source: Aside on stationarity · p. 3

    The analysis assumes the process has a finite trispectrum and is therefore fourth-order stationary, with nonstationarity conditioned on an unknown random time delay.

    Source: Boundedness · p. 2, Conditional nonstationarity · p. 2

Abstract

from arXiv · show

A need to identify modulation of spectral power arises frequently in the analysis of physiological signals. Estimation of power in the relevant bands through filtering and envelope extraction has several limitations: the choice of filter may bias any resulting estimate, while additive Gaussian noise becomes non-Gaussian due to the nonlinearity of envelope computation. The present work considers how spectral decompositions of higher-order cumulants (higher-order spectra, HOS) avoid these limitations, with an emphasis on the use of the trispectrum to identify modulated oscillations. Specifically, it is shown: 1) The trispectrum may be interpreted as a measure of linear dependencies of power across frequencies by viewing it as the cross spectrum of the Wigner-Ville distribution. 2) A particular two-dimensional subdomain of the trispectrum is useful for identifying modulated carriers, recovering essential spectral properties of both the modulating and carrier signals while avoiding the cubic complexity of full trispectrum estimation. A representation of this subdomain, the modulogram, is demonstrated as a tool for identifying and distinguishing different forms of modulation. 3) As a cumulant-derived measure, the modulogram is not biased by additive Gaussian noise. 4) Modulogram phase retains information by which temporal patterns of modulation may be identified. 5) A recently described additive decomposition of HOS (HOSD) further aids identification when applied to the trispectrum. These developments are illustrated with the blind detection of beta bursts in rodent and human local field potential recordings. Finally, the relationship between the present approach and prior techniques of blind identification (BI) through moment maximization, including blind deconvolution and independent component analysis, is considered.

Introduction

Higher-order spectra (HOS) capture cross-frequency dependencies, yet their real-world meaning has been seen as obscure. This work frames modulation identification through the trispectrum and explains how time-invariant higher-order statistics can describe aspects of nonstationary signals.

Introduction

  • HOS are frequency-domain representations of time-shift-invariant higher-order cumulants; they are additive for independent processes, vanish for stationary Gaussian processes, and retain phase information.Source: Introduction · p. 1
  • HOS describe cross-frequency dependencies, but their perceived obscurity of real-world meaning has hampered adoption; this work develops trispectral blind identification of characteristic modulation patterns.Source: Introduction · pp. 1–2

Statement and justification of the problem

Aside on stationarity

  • For ergodic conditionally nonstationary systems, time-invariant moments describe bounded ensemble statistics, allowing conditional second-order nonstationarities to be recovered from time-invariant fourth-order statistics.Source: Aside on stationarity · p. 3

A Statistical View of Amplitude Modulation

The paper models amplitude modulation as a time-varying envelope applied to a stationary carrier and interprets the trispectrum as a measure of linear power dependencies across frequencies. It develops a tractable trispectral subdomain for identifying modulation, while noting assumptions and limitations that affect interpretation and estimation.

A Statistical View of Amplitude Modulation

  • A 19- and 41-Hz sinusoid sum also admits an equivalent reading as a 30-Hz carrier modulated by a 5.5-Hz sinusoidal envelope.The WVD includes the midpoint energy, although it does not appear in the signal’s power spectral density.Source: Figure 1 · p. 4
  • AM identification distinguishes an unmodulated carrier spectrum from a non-constant envelope applied to that carrier, under a model in which the envelope is independent of the carrier.The envelope is modeled as the result of a separate stochastic process.Source: A Statistical View of Amplitude… · pp. 3–4

Trispectrum as a measure of linear cross-frequency dependencies of power

  • The trispectrum is the cross spectrum of the WVD, quantifying linear relationships between power across frequency bands.Treating WVD frequency bands as channels makes this interpretation intuitive; the trispectrum can be computed directly without retaining the WVD intermediate.Source: Trispectrum as a measure of linear… · p. 5

Identifying Amplitude Modulation

  • Real-valued AM induces positive correlations between frequencies equally offset from the carrier, whereas FM induces negative correlations that are reflected in trispectrum phase.For a narrowband carrier, the trispectrum contains a frequency-shifted copy of the envelope trispectrum; distortion becomes negligible when carrier bandwidth is narrow relative to the envelope’s autocorrelation timescale.Source: Identifying Amplitude Modulation · p. 7

The Trispectral Modulogram

  • Estimates restricted to degenerate subdomains can have Gaussian-noise-dependent bias and inflated error, while broadband bias confines resulting power correlations to a lowpass frequency range.The cited criticism concerns estimates based on degenerate domains, including the 0-Hz band of the WVD cross spectrum.Source: The Trispectral Modulogram · pp. 7–8
  • Reparameterizing the trispectrum by carrier frequency and sideband distance identifies a diagonal slice that can detect AM when carrier-band energy is retained, avoiding full trispectrum’s cubic complexity.The model assumes a lowpass modulator so energy remains within the original carrier band; the slice captures dependence among the upper and lower sidebands and carrier band.Source: The Diagonal Slice · pp. 8–9

Additive Decomposition

The modulogram resolves amplitude-modulation spectra better than the kurtogram in tested settings and its phase distinguishes FM from AM; HOSD can be applied to the windowed trispectrum.

Additive Decomposition

  • HOSD applies to windowed higher-order spectra such as the modulogram’s trispectral diagonal slice.The algorithm’s motivation and implementation are described in Appendix C.Source: Additive Decomposition · p. 9
  • The modulogram shows amplitude modulation except when its timescale falls below the estimator’s 0.67 Hz high-pass threshold; the kurtogram misses narrowband cases.Kurtogram performance improves as modulation timescale and fractional bandwidth increase.Source: Additive Decomposition · p. 10
  • Magnitude modulograms are nearly identical for simulated FM and AM, but phase is π rad for FM and 0 rad for AM.The first HOSD feature’s time-frequency pattern also distinguishes the two modulation types.Source: Figure 5 · p. 11

Materials and Methods

The methods apply the modulogram and HOSD to simulated and recorded signals, using selected frequency ranges and windowing strategies to examine modulated oscillations and bursts.

Data segmentation

  • Event-related window delays center segments on events, whereas a common alternative maintains fixed overlap; event-related analyses here use moment spectra for simplicity.The passage notes that event-related windows can be viewed as a degenerate subdomain of a fifth-order cross-moment spectrum, for which a cumulant estimate would include lower-order terms.Source: Window selection · p. 12

Modulogram estimation

  • The modulogram was estimated by averaging Fourier-coefficient products across tapered segments within the trispectral diagonal slice; tricoherence normalizes this estimate.The segments use a DFT, and the modulogram is equivalent to the normalized direct trispectrum estimator within that slice.Source: Modulogram estimation · p. 12

Kurtogram

  • The kurtogram estimates spectral kurtosis across window sizes from windowed spectral power and represents a normalized estimate of a degenerate trispectrum region.The described region has ω1 = ω2.Source: Kurtogram · pp. 12–13

Continuous wavelet analysis

Simulations

  • The simulated signals were formed by multiplying a Gaussian carrier by a non-negative modulating signal produced with a logistic sigmoid nonlinearity.The carrier was band-filtered Gaussian white noise; the modulator was filtered, standardized Gaussian noise passed through a sigmoid with σ = 2.Source: Simulations · p. 13

Beta Burst Detection in Rodent LFP Data

  • The rodent data came from three rats recorded over nine 30-minute sessions each during a reward-conditioning paradigm targeting beta power.The selected beta band was 12–40 Hz; reward delivery followed power exceeding an adaptive threshold.Source: Data set · p. 13
  • For rodent beta-burst analysis, recordings were segmented using both reward-centered and event-independent windows, and HOSD was applied to concatenated sessions.HOSD used carrier frequencies of 10–40 Hz and envelope frequencies of 0.2–5 Hz; burst times were refined against peaks in the reconstructed signal’s analytic envelope.Source: Modulogram estimation · p. 13, Event identification · p. 14, Event characterization · p. 14

Spindle Detection

  • The spindle analysis used thalamic LFP from one propofol-sedated human subject and applied HOSD to the trispectral diagonal slice.The analysis used 8–20 Hz carrier frequencies and 0.1–2 Hz envelope modulation frequencies; the recordings were downsampled to 250 Hz and divided into 10-second Hann-tapered epochs.Source: Spindle Detection · p. 14, Data set · p. 18

Results

The results show that modulograms identify simulated modulation and distinguish AM from FM through phase. In rodent recordings, HOSD recovered beta components and burst patterns across rats; it also recovered spindle-like components in human thalamic data.

AM Detection in Simulated Data

  • In simulated signals, the modulogram resolved modulating-signal spectra better than the kurtogram, which lost resolution and noise robustness with short windows.The kurtogram failed to clearly reveal low-fractional-bandwidth AM in noise; modulogram detection outside DC was limited to modulation timescales shorter than the analysis window.Source: AM Detection in Simulated Data · p. 18

Distinguishing AM and FM

  • AM and FM produced opposite-sign upper- and lower-sideband correlations, with modulogram phase at 0 rad for AM and π rad for FM.Intermediate phases were described as conveying the degree to which upper-sideband amplitude leads lower-sideband amplitude.Source: Distinguishing AM and FM · p. 18

Beta Burst Detection

  • For rat 1, modulograms showed beta modulation in all nine sessions with both windowing methods, and the first HOSD component captured the central modulation band.Subsequent components captured energy at the carrier-band edges, likely reflecting burst-frequency variability; the event-centered modulogram showed somewhat greater relative enhancement.Source: Beta Burst Detection · pp. 18–19
  • HOSD recovered oscillatory components in rats 2 and 3 despite faint or absent modulogram evidence, and burst rates increased over training in all three rats.The authors report that these behavioral patterns agree with Karvat et al..Source: Beta Burst Detection · p. 19
  • The first HOSD component had a median estimated SNR of -4.8 dB in rat 1; the highest-SNR beta components had medians of -9.1 dB in rat 2 and -17 dB in rat 3.Estimated SNR generally decreased with component number; the overall rat 1 median across 10–40 Hz was -9.4 dB.Source: Beta Burst Detection · p. 19

Spindle detection

  • In approximately seven minutes of human thalamic LFP during propofol anesthesia, HOSD recovered spindle-like components around 14 Hz and 12 Hz.The modulogram showed the bands’ characteristic frequencies and slower modulation timescale relative to beta bursts.Source: Spindle detection · p. 19

Discussion

The discussion frames the trispectrum as an interpretable, Gaussian-noise-robust way to study cross-frequency power relationships, while its useful two-dimensional slice and HOSD address practical complexity. The examples and comparisons also clarify where the approach adds value and where its interpretations remain bounded.

Discussion

  • The trispectrum links power across frequencies, resists additive Gaussian-noise bias, and retains phase information that distinguishes modulation patterns.Its two-dimensional diagonal slice reduces the dimensionality of trispectrum analysis when the envelope is lowpass or has a non-zero DC offset; HOSD further projects the estimate into interpretable waveform components.Source: Discussion · pp. 19–20

Relationship to ICA

  • HOSD resembles a spatiotemporal ICA applicable when sources outnumber channels, but its components may reflect one source’s nonstationarity rather than independent sources.Direct cumulant maximization suits super-Gaussian, not sub-Gaussian, sources; distinguishing these interpretations remains open.Source: Relationship to ICA · p. 21

Relationship to Blind deconvolution

  • Blind deconvolution uses static cumulant optimization with a Toeplitz or banded-Toeplitz mixing constraint, and classical diagonal maximization assumes super-Gaussian sources.Amplitude modulation of a Gaussian carrier produces super-Gaussian statistics in the setting considered.Source: Relationship to Blind deconvolution · p. 21

Relationship to Cyclostationarity

  • Unlike peak-based cyclostationary detection, the trispectrum conveys aperiodic as well as periodic modulation information without additive Gaussian-noise bias.Peak detection in a transformed power spectrum may suffice for cyclostationary signals, but not for the general aperiodic or weakly periodic case discussed here.Source: Relationship to Cyclostationarity · p. 22

Relationship to Spectral Kurtosis and the Kurtogram

  • The kurtogram represents carrier frequency and modulation time scale, but its broadband smoothing limits modulating-frequency resolution and prevents frequency-modulation detection.The sideband amplitudes of an FM signal are anticorrelated and cancel when blended by broadband bias.Source: Relationship to Spectral Kurtosis and… · p. 22

On the modulation spectrum

  • The expected squared modulation spectrum is biased by signal power and can suggest a modulation-frequency feature even in phase-randomized stationary Gaussian data.By itself, the modulation spectrum does not account for this bias because it is not defined as an expectation.Source: On the modulation spectrum · p. 23

Identifying oscillatory bursts with the modulogram and HOSD

  • HOSD recovered expected beta-burst features across the three rats despite weak or absent modulogram evidence in rats 2 and 3.For rat 1, the highest median SNR across recovered components was -4.8 dB; maximum estimated SNRs for rats 2 and 3 were -9 dB and -17 dB, respectively. The prominent conditioning-related change was burst rate, consistent with prior reports.Source: Identifying oscillatory bursts with… · p. 23
  • HOSD’s greater apparent sensitivity over raw trispectral estimates may reflect projection onto a lower-dimensional manifold with deterministic spectral structure.The modulogram evidence was weaker for rat 2 and absent for rat 3, while HOSD recovered oscillatory components in the target band.Source: Identifying oscillatory bursts with… · pp. 23–24

Conclusion

The paper connects blind identification to interpretable representations of higher-order statistics. It argues that viewing cumulants through the dynamic properties of lower-order statistics may clarify recovered components and support scientific interpretation.

Conclusion

  • The paper links blind identification to specific, interpretable representations of the underlying higher-order statistics.It suggests that relating higher cumulants of time series to dynamic properties of lower-order statistics may build intuition and offer a priori insight into components recovered through blind identification.Source: Conclusion · p. 24
  • Higher-order-moment techniques can separate non-Gaussian features, but their black-box character limits a priori insight into why an algorithm converges on a solution.The paper notes that ICA and related methods are often used when success is easy to verify or when removing nuisance signals is the practical goal.Source: Conclusion · p. 24
  • Questions of interpretation and reliability continue to constrain the use of ICA and related methods for scientific insight.The paper proposes interpreting higher cumulants through lower-order statistics as one route to clarifying recovered components.Source: Conclusion · p. 24

Main properties

Higher-order spectra preserve time-shift-invariant information about temporal structure and have useful cumulant properties, including Gaussian-noise robustness and additivity. Their estimation choices affect consistency and the time scales of nonstationarity captured; polycoherence normalization is intended to quantify dependence.

Main properties

  • Higher-order spectra retain temporal structure lost by the power spectrum and can recover transient waveforms up to time shift and, at even orders, sign inversion.The second-order spectrum discards phase and yields only energy as a function of frequency.Source: Time invariance · pp. 24–25, Phase preservation · p. 25
  • Cumulant HOS vanish for Gaussian noise and are additive for independent processes, supporting noise-robust estimation and blind decomposition.The paper distinguishes these cumulant properties from moments, which do not generally share them.Source: Indifference to Gaussian noise · p. 25, Additivity for additive mixtures · p. 25, Moment vs. cumulant HOS · p. 25
  • Discarding degenerate subdomains when estimating cumulant HOS loses asymptotically negligible information and yields statistically consistent estimators.For the AM application, this shortcut makes the analysis sensitive to second-order conditional nonstationarity at time scales shorter than the analysis window, unlike earlier approaches sensitive to longer scales.Source: Moment vs. cumulant HOS · p. 25
  • For a deterministic process, HOS factor into products of the signal spectrum; an LTI system driven by stationary white noise can also have deterministic HOS.Transient deterministic signals may be recovered from a limited set of HOS samples, but the text cautions that this does not extend to non-deterministic HOS.Source: Deterministic and non-deterministic HOS · p. 26, HOS inversion · p. 26
  • Under the paper’s ergodic-process view, finite observations act as randomly chosen intervals from a deterministic sequence, a conceptual tool for analyzing estimator properties.The text cautions that this view does not assert that the system itself has deterministic HOS over finite observation windows.Source: HOS of an ergodic process · p. 26

Polycoherence

  • Polycoherence normalization is intended to express statistical dependence on a 0-to-1 scale, though the normalization in Eq. 22 may exceed 1.The paper adopts an alternative normalization according to expected magnitude and notes that other unit-interval normalizations are possible.Source: Polycoherence · pp. 26–27

B Moments and cumulants of bounded ergodic processes

This appendix explains why moment and cumulant spectra agree outside degenerate frequency subdomains for stationary, including conditionally nonstationary, systems. It presents the difference as products of lower-order terms whose frequency partitions determine when those terms vanish.

B Moments and cumulants of bounded ergodic processes

  • For stationary, including conditionally nonstationary, systems, moment and cumulant spectra agree except in degenerate frequency subdomains.Degeneracy occurs when a subset of frequency arguments sums to zero.Source: Cumulant and moment spectra are equal… · p. 28
  • The moment-cumulant difference expands into products of lower-order moments and cumulants, with frequency partitions governing whether terms remain.The passage states that terms after the leading moment term include time-varying moments unless a partitioned frequency subset sums to zero.Source: Cumulant and moment spectra are equal… · p. 28
  • A constructive proof uses an expansion with at most two factors per term, simplifying the frequency partitioning to at most two subsets.The expansion is derived from the relation between the moment-generating and cumulant-generating functions.Source: Cumulant and moment spectra are equal… · p. 28

Cumulant spectra are proper continuous functions

Ergodicity supports estimating higher-order spectra from one time series, while consistency and exclusion of degenerate subdomains underpin cumulant estimates.

Cumulant spectra are proper continuous functions

  • Ergodicity allows sufficiently separated intervals from one time series to stand in for independent ensemble observations when cumulants decay at large lags.Ergodicity also implies cumulant HOS are proper continuous functions.Source: Cumulant spectra are proper continuous… · pp. 28–29
  • Window-overlapped segment-averaged HOS estimates average over N intervals whose shifts and windows adjust to the observation duration T.The estimator expectation converges on the time-invariant moment domain when the window spectrum approaches a unit impulse.Source: Cumulant spectra are proper continuous… · p. 29
  • Excluding degenerate subdomains from a moment estimator can yield a consistent cumulant estimator when its window support converges with the non-degenerate domain.The passage states that lower-order moments vanish for stationary processes in non-degenerate regions and that cumulant estimates are asymptotically independent.Source: Cumulant spectra are proper continuous… · p. 30

C HOSD Algorithm

HOSD uses higher-order spectra to estimate a filter matched to a target waveform, then permits windowing the HOS support to focus or simplify estimation.

C HOSD Algorithm

  • The integration in Eq. 38 produces a matched-filter output raised to K −1, while a clean copy of the target waveform helps the estimate approximate a filter matched to the noisy signal.The matched-filter peak marks the target delay, and exponentiation makes the output more impulse-like.Source: C HOSD Algorithm · p. 31
  • An estimated HOS can replace the unknown waveform HOS in Eq. 38, allowing a filter estimate when the waveform and its delays are initially unknown.Further algorithm details and applications to non-deterministic HOS are referenced in and Appendix C.Source: C HOSD Algorithm · p. 31

HOS windowing

  • A non-negative HOS-domain window can reduce coefficient-estimation cost, emphasize relevant regions, or exclude degenerate subdomains while retaining matched phase within its spectral range.The support need not span the full HOS domain.Source: HOS windowing · p. 31

HOS normalization

The method uses normalization and iterative HOS-based filtering to identify features, refine estimates, threshold detections, and reconstruct and successively remove components.

HOS normalization

  • HOS normalization optimizes the matched filter for Fourier-domain signal-to-noise ratio, but differing choices encode different assumptions about background-noise cross-frequency dependence.An improperly matched amplitude spectrum makes the filter suboptimal for noise suppression, although the expected output can still peak at the correct target delay.Source: HOS normalization · pp. 31–32

Iterative filter estimation

  • Iteratively averaging partial delay filters, applying the average to records, and realigning at output maxima converges efficiently on feature-waveform and matched-filter estimates.This procedure iterates across individual records.Source: Iterative filter estimation · p. 32

Power iteration

  • HOSD’s initial filter is only an approximation to a moment-maximizing filter, and repeated feature/filter re-estimation can refine it through a power-iteration-like procedure.For tensor decomposition with missing non-time-invariant moments, power iteration performs poorly from a random guess but can succeed with the preceding step’s starting approximation.Source: Power iteration · p. 32

Thresholding

  • A heuristic threshold makes the residual static Kth-order cumulant vanish or trades sensitivity against specificity according to the expected Gaussian-noise null distribution.For normalized HOS estimates, thresholding approximates a weighted probability distribution for feature delays.Source: Thresholding · p. 33, Component reconstruction · p. 33

Component reconstruction

  • Convolving the thresholded filter output with the estimated feature waveform and scaling to minimize squared error yields a reconstructed component signal.The thresholded output approximates a delay distribution weighted by amplitude.Source: Component reconstruction · p. 33

Decomposition through serial deflation

  • Subtracting a reconstructed component and repeating HOSD on the residual produces an additive decomposition, stopped by a component-count or cumulant threshold criterion.Deflation can continue until a predetermined number of components is recovered or the detection-filter output reaches the chosen static-cumulant threshold.Source: Decomposition through serial deflation · p. 33

Implementation

The implementation estimates fourth-order spectra from segmented data, restricts estimation to selected trispectrum regions, and applies HOSD to detect components, including weak oscillatory bursts.

Implementation

  • The implementation uses K = 4 and averages HOS coefficients from Fourier-coefficient products across segments.The data are segmented after line-noise removal and artifact rejection.Source: HOS estimation · pp. 33–34
  • The WOSA HOS estimates are asymptotically converging and unbiased as P and J increase with signal length, provided the system is bounded and ergodic.With segment spacing ∆P > 1, estimates are invariant under shifts modulo ∆P; denser spacing restores whole-series shift invariance but is computationally expensive and redundant.Source: HOS estimation · p. 34
  • HOS estimates are normalized to polycoherence to represent frequency relationships and optimize subsequent filter estimates, with mean-magnitude normalization used here.Several alternative normalizations may also be justified.Source: Normalization and HOS Windowing · p. 34
  • Restricting Q to selected carrier and modulation frequencies within the diagonal slice avoids the computational expense of estimating the full trispectrum.The diagonal slice conveys interaction between upper and lower sidebands and the carrier band.Source: Normalization and HOS Windowing · pp. 34–35, Figure 9 · p. 36
  • In rat 2, modulograms gave weak evidence of oscillatory bursting, while rat 3 showed no clear indication; HOSD recovered oscillatory features in both.The authors suggest these cases may represent weak bursting near the limit of detectability.Source: Figure 10 · p. 37
Loading 2609.03172v1…