Data qualification: stationarity and peaks

By José M. Requena Plens · Published Jul 21, 2026, updated Jul 26, 2026

Key references: Bendat & Piersol 2010 Wald & Wolfowitz 1940 Rice 1945

Every average in this documentation - a PSD, a Leq, a GUM budget - assumes the record is stationary: that the process generating it did not drift while it was being measured. Bendat & Piersol devote Section 10.3 of Random Data to qualifying records before analysis, and phonometry.metrology implements its quantitative core: distribution-free trend and stationarity tests with the book’s own acceptance regions, and the Rice statistics - level crossings, apparent frequency, peak rates and heights - that summarize what a qualified Gaussian record looks like and flag one that is not.

1. The reverse arrangement test

Given a sequence of $N$ observations $x_{1}, \dots, x_{N}$ - parameter estimates, segment levels, anything - count the pairs $i < j$ with $x_{i} > x_{j}$ . Each such pair is a reverse arrangement, and for independent observations of one random variable their total $A$ has (B&P Eqs. (4.54)-(4.55))

μ_{A} = \frac{N ( N - 1 )}{4}, σ_{A}^{2} = \frac{N ( 2 N + 5 ) ( N - 1 )}{72},

with no assumption about the distribution of the $x_{i}$ . A monotonic trend pushes $A$ to an extreme (0 for a rising sequence, $N (N - 1) /2$ for a falling one), so the hypothesis of no trend is accepted at significance $α$ when $A$ falls inside a two-sided region - B&P Table A.6, whose $α = 0.05$ rows trend_test reproduces exactly and the conformance suite pins. The book’s Example 4.4 (twenty observations, $A = 86$ , accepted between 64 and 125) runs verbatim:

from phonometry import trend_test

values = [5.2, 6.2, 3.7, 6.4, 3.9, 4.0, 3.9, 5.3, 4.0, 4.6,
          5.9, 6.5, 4.3, 5.7, 3.1, 5.6, 5.2, 3.9, 6.2, 5.0]

res = trend_test(values)                  # B&P Example 4.4
print(res.statistic, res.bounds)          # 86, (64, 125)
print(res.trend_free, round(res.p_value, 3))   # True, 0.586
res.plot()

Two sequences of twenty observations plotted against the sample index: Bendat and Piersol Example 4.4, which fluctuates around five with eighty-six reverse arrangements and is accepted as trend free, and the same fluctuations with an added rising drift climbing from about five to nine, whose thirty-eight reverse arrangements fall below the lower acceptance bound of sixty-four and are rejected as trended at the five percent level

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import trend_test

example = np.array([5.2, 6.2, 3.7, 6.4, 3.9, 4.0, 3.9, 5.3, 4.0, 4.6,
                    5.9, 6.5, 4.3, 5.7, 3.1, 5.6, 5.2, 3.9, 6.2, 5.0])
drifting = example + np.linspace(0.0, 4.0, example.size)   # rising drift

res_flat = trend_test(example)
res_drift = trend_test(drifting)

fig, ax = plt.subplots(figsize=(10, 6))
index = np.arange(1, example.size + 1)
ax.plot(index, res_flat.values, "o-",
        label=f"Example 4.4: A = {res_flat.statistic}, accepted")
ax.plot(index, res_drift.values, "s-",
        label=f"Rising drift: A = {res_drift.statistic}, rejected")
ax.set_xlabel("Sample index")
ax.set_ylabel("Sequence value")
ax.legend()
plt.show()

The p-value comes from the exact null distribution of $A$ (the inversion counts of a random permutation), computed up to $N = 100$ - the range of Table A.6 - and from the book’s normal approximation beyond; the verdict follows the book’s tabulated region. .plot() draws the tested sequence against its sample index with the count, the acceptance region and the verdict in the legend (for method="runs" it also marks the sequence median that classifies each value).

2. Stationarity of a record

The B&P Sec. 10.3.1.1 procedure turns the trend test into a stationarity test for a single time history: divide the record into $N$ equal intervals long enough to be independent, compute a mean square value per interval, and test that sequence. Nothing needs to be known about the record’s bandwidth, averaging distribution or units, and the test works on mean values, rms values or variances just as well (statistic=). A running 20 % gain drift - the book’s Example 10.3 scenario - is caught immediately, while the same noise without the drift passes:

import numpy as np
from phonometry import stationarity_test

fs = 8192.0
n = 1 << 16
noise = np.random.default_rng(42).standard_normal(n)

res = stationarity_test(noise, fs)        # 20 segments, mean squares
print(res.stationary, res.count, res.bounds)   # True, 91, (64, 125)

drifting = noise * np.linspace(1.0, 1.2, n)    # +20 % gain ramp
res = stationarity_test(drifting, fs)
print(res.stationary, res.count)          # False, 7  (upward trend -> low A)
res.plot()

Twenty segment mean square values for two noise records: a steady record whose values fluctuate around one with a reverse arrangement count of ninety-one, accepted as stationary, and the same noise with a twenty percent gain ramp whose segment mean squares climb steadily to one point five, giving only seven reverse arrangements, rejected as nonstationary at the five percent level

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import stationarity_test

fs = 8192.0
n = 1 << 16
steady = np.random.default_rng(42).standard_normal(n)
ramp = np.random.default_rng(42).standard_normal(n) * np.linspace(1.0, 1.2, n)

res_steady = stationarity_test(steady, fs)
res_ramp = stationarity_test(ramp, fs)

fig, ax = plt.subplots(figsize=(10, 6))
index = np.arange(1, res_steady.n_segments + 1)
ax.plot(index, res_steady.segment_values, "o-",
        label=f"Steady noise: A = {res_steady.count}, accepted")
ax.plot(index, res_ramp.segment_values, "s-",
        label=f"+20 % gain ramp: A = {res_ramp.count}, rejected")
ax.set_xlabel("Segment index")
ax.set_ylabel("Segment mean square")
ax.legend()
plt.show()

The result records the per-segment sequence (segment_values, segment_times), the count, the Table A.6 bounds, the exact p_value and the stationary verdict; .plot() draws the sequence with the verdict in the legend. The default of 20 segments matches the book’s worked examples - more segments resolve faster drifts but each interval must stay long against the record’s lowest frequencies for the values to be independent.

Two caveats from the book are worth repeating. A record can be nonstationary with a stationary mean square (a frequency glide, for instance), so pass statistic="mean" or test band-filtered versions when that matters; and the test needs the trend to be slow against the segment length, or it dissolves into the random fluctuation of the segment values.

3. The runs test

The classical companion (tabulated in the third edition of Random Data; the exact distribution is Wald & Wolfowitz 1940) classifies each value as above or below the sequence median and counts runs of like classification. A trend or slow drift produces few long runs; rapid alternation produces too many. method="runs" applies it with the exact conditional distribution at any $N$ - for twenty values split 10/10 the acceptance region at $α = 0.05$ is the classical $(6, 15]$ :

import numpy as np
from phonometry import trend_test

rng = np.random.default_rng(3)
res = trend_test(rng.standard_normal(40), method="runs")
print(res.statistic, res.bounds, res.trend_free)
res.plot()   # the sequence about its median with the runs verdict

alternating = np.tile([1.0, -1.0], 10)   # 20 runs: rejected the other way
print(trend_test(alternating, method="runs").trend_free)   # False

Two panels of sequences classified about their median, points above in blue and below in red. Left: forty Gaussian observations with twenty runs inside the acceptance region from fourteen to twenty-seven, verdict trend-free. Right: a strictly alternating twenty-value sequence also giving twenty runs, but against the acceptance region from six to fifteen for twenty values, verdict rejected for too-rapid alternation

The two-sided verdict of the runs test: 20 runs is unremarkable for 40 Gaussian observations (left, accepted), but the same count from a strictly alternating 20-value sequence exceeds its upper bound (right, rejected): too many runs is as non-random as too few.

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import trend_test

rng = np.random.default_rng(3)
sequences = [rng.standard_normal(40), np.tile([1.0, -1.0], 10)]

# One line per sequence — the tested values with the verdict:
trend_test(sequences[0], method="runs").plot()
plt.show()

# Both verdicts side by side, classified about the median:
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
for ax, seq in zip(axes, sequences):
    res = trend_test(seq, method="runs")
    idx = np.arange(1, seq.size + 1)
    median = np.median(seq)
    above = seq > median
    ax.plot(idx, seq, "0.7", lw=0.7)
    ax.scatter(idx[above], seq[above], label="above the median")
    ax.scatter(idx[~above], seq[~above], color="r",
               label="below the median")
    ax.axhline(median, color="k", linestyle="--")
    lo, hi = res.bounds
    verdict = "trend-free" if res.trend_free else "rejected"
    ax.set_title(f"r = {res.statistic} runs, accept ({lo}, {hi}]: {verdict}")
    ax.set(xlabel="Sample index", ylabel="Sequence value")
    ax.legend()
plt.show()

The reverse arrangement test is the more powerful of the two against the monotonic trends that dominate practice (B&P Sec. 4.5.2), which is why it is the default everywhere; the runs test adds sensitivity to non-monotonic clustering.

4. Level crossings and the apparent frequency

For a zero-mean Gaussian record with one-sided autospectrum $G (f)$ , Rice’s classical results give the expected rate of zero crossings (both slopes, B&P Eq. (5.195)) from the plain frequency moments $m_{k} = \int f^{k} G (f) d f$ :

N_{0} = \frac{1}{π} \frac{σ _{v}}{σ _{x}} = 2 \frac{m _{2}}{m _{0}}, N_{a} = N_{0} e^{- a^{2} /2 σ_{x}^{2}},

where $N_{a}$ is the crossing rate of level $a$ (Eq. (5.196)). $N_{0} /2$ is the record’s apparent frequency: a 60 Hz sine crosses zero 120 times per second, low-pass noise of bandwidth $B$ gives $N_{0} = 2 B / 3$ (an apparent frequency of $0.58 B$ , Example 5.12) and a band centred on $f_{c}$ gives $N_{0} = 2 f_{c}^{2} + B^{2} /12$ (Example 5.13). level_crossing_rate counts the actual crossings of each level and puts the Rice curve next to them, taking the moments from the record’s own Welch autospectrum:

import numpy as np
from phonometry import level_crossing_rate

fs = 8192.0
t = np.arange(1 << 16) / fs
x = np.sin(2 * np.pi * 60.0 * t)          # a 60 Hz sine ...

res = level_crossing_rate(x, fs)
print(round(res.zero_crossing_rate, 1))   # ... has 120 zeros per second
print(round(res.apparent_frequency, 1))   # 60.0
res.plot()

Measured level-crossing rates of a bandlimited Gaussian noise record between eight hundred and twelve hundred hertz, plotted as dots against the crossing level from minus three point five to plus three point five signal units on a logarithmic rate axis, falling from about two thousand crossings per second at level zero to a few per second at three sigma, with the Rice exponential curve passing through every measured point

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import level_crossing_rate

fs = 20480.0
n = 1 << 19
rng = np.random.default_rng(0)
freqs = np.fft.rfftfreq(n, 1 / fs)        # bandlimited Gaussian noise:
spec = rng.standard_normal(freqs.size) + 1j * rng.standard_normal(freqs.size)
spec[(freqs < 800.0) | (freqs > 1200.0)] = 0.0
x = np.fft.irfft(spec, n)

res = level_crossing_rate(x, fs, levels=np.linspace(-3.5, 3.5, 29) * np.std(x))

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(res.levels, res.rice_rates, label="Rice (Eq. 5.196)")
ax.plot(res.levels, res.rates, "o", label="Measured")
ax.set_yscale("log")
ax.set_xlabel("Level a [signal units]")
ax.set_ylabel("Crossings per second [1/s]")
ax.legend()
plt.show()

The Rice curve holds for Gaussian records, and that is precisely its second use: measured rates that fall systematically off the curve are a quick non-Gaussianity screen (B&P Sec. 5.5.1.1) - clipping shows up as missing high-level crossings, impulsive contamination as an excess. Count-based rates need the record comfortably oversampled: a crossing between two samples of the same sign goes uncounted.

5. Peaks: rates, the irregularity factor and heights

The same moments fix the expected rate of local maxima $M = (1/2 π) (σ_{a} / σ_{v}) = m_{4} / m_{2}$ (Eq. (5.211)) and with it the dimensionless irregularity factor

r = \frac{N _{0}}{2 M} = \frac{m _{2}}{m _{0} m _{4}} \in (0, 1],

the single number that fixes the distribution of peak heights (B&P Sec. 5.5.4). At $r = 1$ - narrow bandwidth data, one maximum per zero-crossing cycle - peaks are Rayleigh distributed: $Prob [peak > a] = e^{- a^{2} /2 σ_{x}^{2}}$ (Eq. (5.206)), which is the one-in-3000 chance of a peak beyond $4 σ$ of B&P Example 5.14. As $r \to 0$ ever more ripples ride on each cycle, negative maxima appear, and the peak heights approach the plain Gaussian amplitude distribution. In between, Rice’s mixture (Eqs. (5.217)/(5.223)) interpolates the two, available as peak_exceedance() / peak_density() on the result:

import numpy as np
from phonometry import peak_statistics

fs = 8192.0
t = np.arange(1 << 16) / fs
x = np.sin(2 * np.pi * 60.0 * t)          # narrowband: r is essentially 1

res = peak_statistics(x, fs)
print(round(res.irregularity_factor, 3))       # 0.997
print(res.peak_exceedance(4.0))                # exp(-8): about 1 in 3000
res.plot()   # empirical exceedance against the Rice mixture (figure below)

$Peak-height exceedance probability of low-pass Gaussian noise on a logarithmic axis against the standardized peak height from minus two point five to four point five: the empirical staircase from half a million samples follows the Rice mixture curve for irregularity factor zero point seven four six, clearly below the dashed Rayleigh limit and above the dotted Gaussian limit, with a visible fraction of negative maxima$ $Peak-height exceedance probability of low-pass Gaussian noise on a logarithmic axis against the standardized peak height from minus two point five to four point five: the empirical staircase from half a million samples follows the Rice mixture curve for irregularity factor zero point seven four six, clearly below the dashed Rayleigh limit and above the dotted Gaussian limit, with a visible fraction of negative maxima$

Show the code for this figure

import matplotlib.pyplot as plt
import numpy as np
from phonometry import peak_statistics
from phonometry.metrology.random_data import _rice_peak_exceedance

fs = 20480.0
n = 1 << 19
rng = np.random.default_rng(3)
freqs = np.fft.rfftfreq(n, 1 / fs)        # low-pass noise: r = sqrt(5)/3
spec = rng.standard_normal(freqs.size) + 1j * rng.standard_normal(freqs.size)
spec[freqs > 2000.0] = 0.0
x = np.fft.irfft(spec, n)

res = peak_statistics(x, fs)
peaks = res.peak_values
empirical = 1.0 - np.arange(1, peaks.size + 1) / peaks.size
z = np.linspace(-2.5, 4.5, 400)

fig, ax = plt.subplots(figsize=(10, 6))
ax.plot(z, _rice_peak_exceedance(z, 1.0), "--", label="Rayleigh (r = 1)")
ax.plot(z, _rice_peak_exceedance(z, 0.0), ":", label="Gaussian (r = 0)")
ax.plot(z, res.peak_exceedance(z),
        label=f"Rice (r = {res.irregularity_factor:.3f})")
ax.plot(peaks, empirical, drawstyle="steps-post", label="Empirical")
ax.set_yscale("log")
ax.set_ylim(1e-5, 1.5)
ax.set_xlabel(r"Standardized peak height $z = a/\sigma_x$")
ax.set_ylabel("Prob[peak > z]")
ax.legend()
plt.show()

For an ideal low-pass band the closed forms give $r = 5 /3 = 0.745$ , and the measured value lands on it. The irregularity factor is the standard bridge to fatigue and vibro-acoustic damage estimation, where it selects the cycle-counting correction between the narrow-band Rayleigh assumption and broad-band rainflow corrections. One practical warning: $m_{4}$ weights the spectrum by $f^{4}$ , so wideband instrumentation noise far above the physical band silently inflates $M$ and deflates $r$ - band-limit the record to the physically meaningful range first.

Where this fits

Qualification comes before the statistics this section’s other pages compute: the chi-square confidence interval of a Welch PSD and the random-error formulas of the correlation estimators all assume the record is stationary, as does the very idea of the Leq of a measurement in Levels. When a record fails the test, split it at the change (the segment_values sequence shows where), analyse the pieces, or move to the short-time views - the calibrated spectrogram - that do not assume stationarity. And when it passes, the GUM machinery can propagate the remaining random error of every averaged estimate with a clean conscience.

What this guide covers

Covered. Bendat & Piersol Section 4.5.2 with Table A.6 (the reverse arrangement trend test, trend_test, its exact and normal-approximation p-value and Example 4.4), the Wald & Wolfowitz (1940) runs test (trend_test(method="runs")), the Section 10.3.1.1 segment mean-square stationarity procedure with Example 10.3 (stationarity_test), and the Rice (1945) level-crossing and peak statistics of Section 5.5 (level_crossing_rate, peak_statistics, the irregularity factor and its Rayleigh/Gaussian peak-height limits).

Not covered. Section 10.3 of the book is titled “data qualification” more broadly and also covers classification, validation and editing procedures; phonometry.metrology implements only its quantitative core, the segment mean-square stationarity test. Classifying a record’s type, validating it against physical limits, or editing out glitches stay manual steps the way the book describes them, outside this module.

References

Bendat, J. S., & Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4th ed.). Wiley. https://doi.org/10.1002/9781118032428Section 4.5.2 with Table A.6 (the nonparametric reverse arrangement trend test, its mean and variance, the tabulated percentage points and Example 4.4), Section 10.3 (data qualification: classification, validation, editing; the segment mean-square stationarity procedure of 10.3.1.1 and Example 10.3) and Section 5.5 (level crossings and peak values: Eqs. (5.195)-(5.196), (5.206), (5.211), (5.217)-(5.223) and Examples 5.12-5.15). ISBN 978-0-470-24877-5. The runs companion appeared in the third edition's percentage points of the run distribution.
Rice, S. O. (1945). Mathematical analysis of random noise. The Bell System Technical Journal, 24(1), 46-156. https://doi.org/10.1002/j.1538-7305.1945.tb00453.xThe original derivations of the expected level-crossing and maxima rates and the peak distribution of Gaussian noise that Bendat & Piersol Section 5.5 presents (their Ref. 6; Parts I-II are in volume 23, 1944).
Wald, A., & Wolfowitz, J. (1940). On a test whether two samples are from the same population. The Annals of Mathematical Statistics, 11(2), 147-162. https://doi.org/10.1214/aoms/1177731909The exact conditional distribution of the number of runs, from which the runs-about-the-median acceptance regions are computed.

Created and maintained byJosé Manuel Requena Plens· GitHub· PyPI· All projects

Data qualification: stationarity and peaks

1. The reverse arrangement test

2. Stationarity of a record

3. The runs test

4. Level crossings and the apparent frequency

5. Peaks: rates, the irregularity factor and heights

Where this fits

What this guide covers

See also

References