Source: https://jmrplens.github.io/phonometry/es/guides/data-qualification/
Todo promedio de esta documentación - una [PSD](/phonometry/es/guides/spectral-analysis/), un
[Leq](/phonometry/es/guides/levels/), un
[balance GUM](/phonometry/es/guides/gum-uncertainty/) - supone que el registro es
**estacionario**: que el proceso que lo genera no derivó mientras se medía.
Bendat y Piersol dedican la Sección 10.3 de *Random Data* a calificar los
registros antes del análisis, y `phonometry.metrology` implementa su núcleo
cuantitativo: **tests de tendencia y de estacionariedad** libres de
distribución con las regiones de aceptación del propio libro, y las
**estadísticas de Rice** - cruces por nivel, frecuencia aparente, tasas y
alturas de pico - que resumen el aspecto de un registro gaussiano calificado
y delatan al que no lo es.
## 1. El test de inversiones de orden
Dada una secuencia de $N$ observaciones $x_1, \dots, x_N$ - estimaciones de
parámetros, niveles por segmento, lo que sea - cuéntense los pares $i < j$
con $x_i > x_j$. Cada par es una **inversión de orden** (*reverse
arrangement*), y para observaciones independientes de una misma variable
aleatoria su total $A$ tiene (B&P Ecs. (4.54)-(4.55))
$$
\mu_A = \frac{N(N-1)}{4}, \qquad
\sigma_A^2 = \frac{N(2N+5)(N-1)}{72},
$$
sin ninguna suposición sobre la distribución de las $x_i$. Una tendencia
monótona empuja $A$ a un extremo (0 si la secuencia sube, $N(N-1)/2$ si
baja), así que la hipótesis de *ausencia de tendencia* se acepta al nivel de
significación $\alpha$ cuando $A$ cae dentro de una región bilateral - la
Tabla A.6 de B&P, cuyas filas de $\alpha = 0{,}05$ reproduce exactamente
`trend_test` y fija la suite de conformidad. El Ejemplo 4.4 del libro
(veinte observaciones, $A = 86$, aceptado entre 64 y 125) corre tal cual:
```python
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 Ejemplo 4.4
print(res.statistic, res.bounds) # 86, (64, 125)
print(res.trend_free, round(res.p_value, 3)) # True, 0.586
res.plot(language="es")
```
Mostrar el código de esta figura
```python
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) # deriva ascendente
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"Ejemplo 4.4: A = {res_flat.statistic}, aceptado")
ax.plot(index, res_drift.values, "s-",
label=f"Deriva ascendente: A = {res_drift.statistic}, rechazado")
ax.set_xlabel("Índice de muestra")
ax.set_ylabel("Valor de la secuencia")
ax.legend()
plt.show()
```
El p-valor sale de la distribución nula exacta de $A$ (los conteos de
inversiones de una permutación aleatoria), calculada hasta $N = 100$ - el
alcance de la Tabla A.6 - y de la aproximación normal del libro más allá;
el veredicto sigue la región tabulada del libro. `.plot()` dibuja la
secuencia analizada frente a su índice de muestra con el conteo, la región
de aceptación y el veredicto en la leyenda (con `method="runs"` marca
además la mediana de la secuencia que clasifica cada valor).
## 2. Estacionariedad de un registro
El procedimiento de la Sec. 10.3.1.1 de B&P convierte el test de tendencia
en un test de estacionariedad para una única historia temporal: divídase el
registro en $N$ intervalos iguales lo bastante largos para ser
independientes, calcúlese una **media cuadrática por intervalo** y pruébese
esa secuencia. No hace falta saber nada del ancho de banda del registro, de
su distribución ni de sus unidades, y el test funciona igual de bien sobre
medias, valores eficaces o varianzas (`statistic=`). Una deriva de ganancia
del 20 % - el escenario del Ejemplo 10.3 del libro - se detecta de
inmediato, mientras que el mismo ruido sin deriva pasa:
```python
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 segmentos, medias cuadráticas
print(res.stationary, res.count, res.bounds) # True, 91, (64, 125)
drifting = noise * np.linspace(1.0, 1.2, n) # rampa de ganancia del +20 %
res = stationarity_test(drifting, fs)
print(res.stationary, res.count) # False, 7 (tendencia al alza -> A bajo)
res.plot(language="es")
```
Mostrar el código de esta figura
```python
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"Ruido estable: A = {res_steady.count}, aceptado")
ax.plot(index, res_ramp.segment_values, "s-",
label=f"Rampa de ganancia del +20 %: A = {res_ramp.count}, rechazado")
ax.set_xlabel("Índice de segmento")
ax.set_ylabel("Media cuadrática por segmento")
ax.legend()
plt.show()
```
El resultado registra la secuencia por segmento (`segment_values`,
`segment_times`), el conteo, los `bounds` de la Tabla A.6, el `p_value`
exacto y el veredicto `stationary`; `.plot()` dibuja la secuencia con el
veredicto en la leyenda. El valor por defecto de 20 segmentos coincide con
los ejemplos resueltos del libro - más segmentos resuelven derivas más
rápidas, pero cada intervalo debe seguir siendo largo frente a las
frecuencias más bajas del registro para que los valores sean
independientes.
Merece la pena repetir dos salvedades del libro. Un registro puede ser no
estacionario con media cuadrática estacionaria (un barrido de frecuencia,
por ejemplo), así que pásese `statistic="mean"` o pruébense versiones
filtradas por bandas cuando importe; y el test necesita que la *tendencia*
sea lenta frente a la longitud del segmento, o se disuelve en la
fluctuación aleatoria de los valores de segmento.
## 3. El test de rachas
El complemento clásico (tabulado en la tercera edición de *Random Data*; la
distribución exacta es de Wald y Wolfowitz, 1940) clasifica cada valor como
por encima o por debajo de la mediana de la secuencia y cuenta **rachas**
de clasificación igual. Una tendencia o deriva lenta produce pocas rachas
largas; la alternancia rápida produce demasiadas. `method="runs"` lo aplica
con la distribución condicional exacta para cualquier $N$ - para veinte
valores repartidos 10/10 la región de aceptación con $\alpha = 0{,}05$ es la
clásica $(6, 15]$:
```python
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(language="es") # la secuencia sobre su mediana con el veredicto
alternating = np.tile([1.0, -1.0], 10) # 20 rachas: rechazado por el otro lado
print(trend_test(alternating, method="runs").trend_free) # False
```
*El veredicto bilateral del test de rachas: 20 rachas no llaman la atención
en 40 observaciones gaussianas (izquierda, aceptado), pero el mismo
recuento en una secuencia alternante de 20 valores supera su cota superior
(derecha, rechazado): demasiadas rachas es tan poco aleatorio como
demasiado pocas.*
Mostrar el código de esta figura
```python
from phonometry import trend_test
rng = np.random.default_rng(3)
sequences = [rng.standard_normal(40), np.tile([1.0, -1.0], 10)]
# Una línea por secuencia: los valores analizados con el veredicto:
trend_test(sequences[0], method="runs").plot(language="es")
plt.show()
# Los dos veredictos lado a lado, clasificados respecto a la mediana:
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="sobre la mediana")
ax.scatter(idx[~above], seq[~above], color="r",
label="bajo la mediana")
ax.axhline(median, color="k", linestyle="--")
lo, hi = res.bounds
verdict = "sin tendencia" if res.trend_free else "rechazada"
ax.set_title(f"r = {res.statistic} rachas, aceptación ({lo}, {hi}]: {verdict}")
ax.set(xlabel="Índice de muestra", ylabel="Valor de la secuencia")
ax.legend()
plt.show()
```
El test de inversiones de orden es el más potente de los dos frente a las
tendencias monótonas que dominan la práctica (B&P Sec. 4.5.2), y por eso es
el valor por defecto en todas partes; el de rachas añade sensibilidad al
agrupamiento no monótono.
## 4. Cruces por nivel y frecuencia aparente
Para un registro gaussiano de media cero con autoespectro unilateral
$G(f)$, los resultados clásicos de Rice dan la tasa esperada de cruces por
cero (ambas pendientes, B&P Ec. (5.195)) a partir de los momentos
frecuenciales $m_k = \int f^k\,G(f)\,\mathrm{d}f$:
$$
N_0 = \frac{1}{\pi}\frac{\sigma_v}{\sigma_x} = 2\sqrt{\frac{m_2}{m_0}},
\qquad
N_a = N_0\, e^{-a^2/2\sigma_x^2},
$$
donde $N_a$ es la tasa de cruces del nivel $a$ (Ec. (5.196)). $N_0/2$ es la
**frecuencia aparente** del registro: un seno de 60 Hz cruza el cero 120
veces por segundo, el ruido paso bajo de ancho de banda $B$ da
$N_0 = 2B/\sqrt{3}$ (una frecuencia aparente de $0{,}58B$, Ejemplo 5.12) y
una banda centrada en $f_c$ da $N_0 = 2\sqrt{f_c^2 + B^2/12}$
(Ejemplo 5.13). `level_crossing_rate` cuenta los cruces reales de cada
nivel y pone al lado la curva de Rice, tomando los momentos del propio
[autoespectro de Welch](/phonometry/es/guides/spectral-analysis/) del registro:
```python
from phonometry import level_crossing_rate
fs = 8192.0
t = np.arange(1 << 16) / fs
x = np.sin(2 * np.pi * 60.0 * t) # un seno de 60 Hz...
res = level_crossing_rate(x, fs)
print(round(res.zero_crossing_rate, 1)) # ...tiene 120 ceros por segundo
print(round(res.apparent_frequency, 1)) # 60.0
res.plot(language="es")
```
Mostrar el código de esta figura
```python
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) # ruido gaussiano limitado en banda:
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 (Ec. 5.196)")
ax.plot(res.levels, res.rates, "o", label="Medida")
ax.set_yscale("log")
ax.set_xlabel("Nivel a [unidades de la señal]")
ax.set_ylabel("Cruces por segundo [1/s]")
ax.legend()
plt.show()
```
La curva de Rice vale para registros gaussianos, y ese es precisamente su
segundo uso: unas tasas medidas que se apartan sistemáticamente de la curva
son un cribado rápido de no gaussianidad (B&P Sec. 5.5.1.1) - el recorte
aparece como cruces de nivel alto que faltan, la contaminación impulsiva
como un exceso. Las tasas basadas en conteo necesitan el registro
holgadamente sobremuestreado: un cruce entre dos muestras del mismo signo
queda sin contar.
## 5. Picos: tasas, factor de irregularidad y alturas
Los mismos momentos fijan la tasa esperada de máximos locales
$M = (1/2\pi)(\sigma_a/\sigma_v) = \sqrt{m_4/m_2}$ (Ec. (5.211)) y con ella
el **factor de irregularidad** adimensional
$$
r = \frac{N_0}{2M} = \frac{m_2}{\sqrt{m_0\, m_4}} \in (0, 1],
$$
el número único que fija la distribución de las alturas de pico
(B&P Sec. 5.5.4). Con $r = 1$ - datos de banda estrecha, un máximo por
ciclo de cruce por cero - los picos siguen una distribución de
**Rayleigh**: $\mathrm{Prob}[\text{pico} > a] = e^{-a^2/2\sigma_x^2}$
(Ec. (5.206)), la probabilidad de uno entre 3000 de un pico más allá de
$4\sigma$ del Ejemplo 5.14 de B&P. Cuando $r \to 0$ cada ciclo lleva cada
vez más rizos encima, aparecen máximos negativos y las alturas de pico se
acercan a la distribución de amplitud **gaussiana** simple. Entre medias,
la mezcla de Rice (Ecs. (5.217)/(5.223)) interpola ambas, disponible como
`peak_exceedance()` / `peak_density()` en el resultado:
```python
from phonometry import peak_statistics
fs = 8192.0
t = np.arange(1 << 16) / fs
x = np.sin(2 * np.pi * 60.0 * t) # banda estrecha: r es en esencia 1
res = peak_statistics(x, fs)
print(round(res.irregularity_factor, 3)) # 0.997
print(res.peak_exceedance(4.0)) # exp(-8): aprox. 1 entre 3000
res.plot(language="es") # excedencia empírica contra la mezcla de Rice (figura de abajo)
```
Mostrar el código de esta figura
```python
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) # ruido paso bajo: 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="Gaussiana (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="Empírica")
ax.set_yscale("log")
ax.set_ylim(1e-5, 1.5)
ax.set_xlabel(r"Altura de pico estandarizada $z = a/\sigma_x$")
ax.set_ylabel("Prob[pico > z]")
ax.legend()
plt.show()
```
Para una banda paso bajo ideal las formas cerradas dan
$r = \sqrt{5}/3 = 0{,}745$, y el valor medido aterriza encima. El factor de
irregularidad es el puente estándar hacia la estimación de fatiga y de daño
vibroacústico, donde selecciona la corrección de conteo de ciclos entre la
hipótesis de Rayleigh de banda estrecha y las correcciones rainflow de
banda ancha. Un aviso práctico: $m_4$ pondera el espectro por $f^4$, así
que el ruido de instrumentación de banda ancha muy por encima de la banda
física infla $M$ y deprime $r$ sin avisar - limítese antes el registro a la
banda físicamente significativa.
## Dónde encaja esto
La calificación va *antes* de las estadísticas que calculan las demás
páginas de esta sección: el intervalo de confianza chi-cuadrado de una
[PSD de Welch](/phonometry/es/guides/spectral-analysis/) y las fórmulas de error
aleatorio de los [estimadores de correlación](/phonometry/es/guides/correlation-delay/)
suponen que el registro es estacionario, igual que la idea misma de *el*
Leq de una medición en [Niveles](/phonometry/es/guides/levels/). Cuando un
registro no pasa el test, divídase en el cambio (la secuencia
`segment_values` muestra dónde), analícense los trozos, o pásese a las
vistas de tiempo corto - el
[espectrograma calibrado](/phonometry/es/guides/time-frequency/) - que no suponen
estacionariedad. Y cuando pasa, la
[maquinaria GUM](/phonometry/es/guides/gum-uncertainty/) puede propagar el error
aleatorio *restante* de cada estimación promediada con la conciencia
tranquila.
## Qué cubre esta guía
**Cubierto.** La sección 4.5.2 con la Tabla A.6 de Bendat y Piersol (el
test de tendencia por inversiones de orden, `trend_test`, su valor p exacto
y por aproximación normal, y el Ejemplo 4.4), el test de rachas de Wald y
Wolfowitz (1940) (`trend_test(method="runs")`), el procedimiento de
estacionariedad por medias cuadráticas de segmento de la sección 10.3.1.1
con el Ejemplo 10.3 (`stationarity_test`), y las estadísticas de cruces por
nivel y de picos de Rice (1945) de la sección 5.5 (`level_crossing_rate`,
`peak_statistics`, el factor de irregularidad y sus límites Rayleigh y
gaussiano para la altura de pico).
**No cubierto.** La sección 10.3 del libro se titula "calificación de
datos" en un sentido más amplio y también cubre la clasificación, la
validación y la edición de registros; `phonometry.metrology` implementa
solo su núcleo cuantitativo, el test de estacionariedad por medias
cuadráticas. Clasificar el tipo de un registro, validarlo frente a límites
físicos o editar transitorios espurios siguen siendo pasos manuales tal
como los describe el libro, fuera de este módulo.
## Véase también
- Referencia de la API: [`metrology.random_data`](/phonometry/es/reference/api/metrology/random-data/).