Source: https://jmrplens.github.io/phonometry/es/signals/metrology/data-qualification/
Todo promedio de esta documentación - una [PSD](/phonometry/es/signals/spectra/spectral-analysis/), un
[Leq](/phonometry/es/signals/levels/levels/), un
[balance GUM](/phonometry/es/signals/metrology/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 cualificar 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 cualificado
y delatan al que no lo es.
La cualificación es una decisión, no un estadístico, así que se lee mejor
como un flujo: segmentar el registro, contar inversiones, y dejar que la
región de la Tabla A.6 diga si se puede promediar.
## 0. Antes de los tests: validar el registro
La Sección 10.3 de Bendat y Piersol cubre tres actividades y este módulo
implementa una de ellas. La *cualificación*, los tests de más abajo, pregunta si
un registro es estacionario. La *validación* pregunta si es siquiera un registro
de aquello, y va primero, porque cada defecto de la lista siguiente o suspende el
test de estacionariedad por el motivo equivocado o, peor, lo aprueba. La
validación sigue siendo manual, tal como la describe el libro, pero son seis
líneas de numpy y una decisión por cada una:
```python
rng = np.random.default_rng(0)
x = 0.15 * rng.standard_normal(1 << 16) # tu registro, en unidades digitales
full_scale = 1.0 # 1.0 en audio float, 32768 en int16
clipped = int(np.count_nonzero(np.abs(x) >= 0.999 * full_scale))
dead = int(np.max(np.diff(np.flatnonzero(np.diff(x) != 0.0), prepend=0)))
offset = float(np.mean(x) / np.std(x))
spikes = int(np.count_nonzero(np.abs(x - np.mean(x)) > 5 * np.std(x)))
kurtosis = float(np.mean((x - np.mean(x)) ** 4) / np.std(x) ** 4)
print(clipped, dead, round(offset, 4), spikes, round(kurtosis, 2))
```
| Defecto | La comprobación | La decisión |
| :--- | :--- | :--- |
| Recorte | muestras iguales o superiores a 0,999 del fondo de escala | más de un puñado y el registro se edita o se rechaza: el recorte desinfla los cruces de nivel alto del §4 e infla todas las medias cuadráticas |
| Caídas, canal muerto | la racha más larga de muestras consecutivas idénticas | cualquier racha de más de unos pocos milisegundos es un fallo de transmisión, no una señal |
| Desplazamiento de continua, deriva infrasónica | la media del registro frente a su desviación típica, y la parte de la PSD por debajo de 20 Hz frente a la banda de interés | elimina la media y filtra en paso alto antes del §4, cuyas fórmulas suponen un proceso de media nula |
| Picos, contaminación transitoria | muestras más allá de 5σ, o una curtosis muy por encima de 3 en un registro supuestamente gaussiano | edítalos y dilo, o cualifica el registro por trozos |
| Zumbido de red | picos a 50 o 60 Hz y sus armónicos en una [PSD de Welch](/phonometry/es/signals/spectra/spectral-analysis/) | un fallo de la cadena, no de la fuente; arregla las masas, no los datos |
El registro sintético de arriba imprime `0 1 0.0024 0 3.02`, que es el aspecto de
lo limpio: ninguna muestra recortada, ninguna muestra repetida, un desplazamiento
de 0,002σ, ni una sola excursión más allá de 5σ en 65 536 muestras gaussianas y
una curtosis de 3,0, como debe tener un registro gaussiano.
Junto a los números, los metadatos de adquisición tienen que viajar con el
fichero o nada de esto se puede reconstruir después: frecuencia de muestreo,
profundidad de bits, ajuste de ganancia, el
[factor de calibración](/phonometry/es/signals/metrology/calibration/) y la toma de
la que salió, micrófono y pantalla antiviento, y, en exteriores, velocidad del
viento y temperatura. La regla es sencillamente que la validación precede a la
cualificación: primero editar o rechazar, después probar la estacionariedad.
## 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 - cuenta 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 batería de conformidad. El Ejemplo 4.4 del libro
(veinte observaciones, $A = 86$, aceptado en $(64, 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]: semiabierto
print(res.trend_free, round(res.p_value, 3)) # True, 0.586
res.plot(language="es")
```
*El test no mira nunca la pendiente ajustada, solo el orden de los valores. A la
izquierda, el propio Ejemplo 4.4 de Bendat y Piersol: $A = 86$, holgadamente
dentro de la región de aceptación $(64, 125]$, $p = 0{,}586$. A la derecha, las
mismas veinte fluctuaciones con una deriva añadida que las lleva de 5,2 a 9,0:
las fluctuaciones no cambian, pero el orden sí, y $A$ cae a 38, por debajo del
límite inferior, con $p = 0{,}0001$. Lo que empuja $A$ hacia cero es una deriva
monótona, que es por lo que el recuento funciona sin conocer la distribución.*
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()
```
`bounds` es **semiabierto**: la hipótesis de ausencia de tendencia se acepta
cuando `lower < A <= upper`, así que $A = 125$ pasa y $A = 64$ no. El test de
rachas del §3 usa el mismo convenio, y el veredicto de estacionariedad del §2
también.
El valor p 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).
**Cómo leer el veredicto.** Es un contraste de hipótesis, no una medida.
`trend_free = True` significa que los datos *no dan indicio* de tendencia al
nivel elegido, que es mucho más débil que decir «el registro es estacionario»;
$p = 0{,}586$ dice que el conteo observado es de lo más corriente en una
secuencia sin tendencia, y una $p$ apenas por encima de 0,05 no dice casi nada en
ninguna dirección. La consecuencia que muerde en la práctica es la
multiplicidad: con $\alpha = 0{,}05$, uno de cada veinte registros estacionarios
se rechaza por construcción, así que cribar las 31 bandas de tercio de octava de
un registro perfectamente estacionario produce uno o dos rechazos como resultado
*esperado*. Una sola banda que falla no es, por tanto, indicio de nada; una racha
de bandas contiguas que fallan, o un fallo del registro de banda ancha, sí lo es.
El mando es `alpha=`: apriétalo cuando se criben muchas bandas, aflójalo cuando
haya que cualificar un registro de forma conservadora, y las cotas de aceptación
se mueven con él, porque salen de la distribución nula exacta y no de una
consulta a una tabla.
## 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: divide el
registro en $N$ intervalos iguales lo bastante largos para ser
independientes, calcula una **media cuadrática por intervalo** y prueba
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")
```
*La misma secuencia de veinte segmentos cuenta la estacionariedad con la misma
aritmética. Las medias cuadráticas por segmento del registro estable se quedan
entre 0,982 y 1,034 y dan $A = 91$, dentro de $(64, 125]$. La rampa de ganancia
del 20 % las sube de 1,016 a 1,476, un cambio de un factor 1,45 que nadie
llamaría dramático a ojo sobre una traza temporal, y $A$ se desploma a 7. Bastan
veinte números y un recuento de pares comparado con una tabla; no entra nada del
ancho de banda ni de las unidades del registro.*
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. Conviene fijarse en que el mismo conteo de inversiones
responde a dos nombres de atributo distintos en las dos clases de resultado,
porque solo una de ellas tiene una secuencia por segmento que llevar:
| | `trend_test` → `TrendTestResult` | `stationarity_test` → `StationarityTestResult` |
| :--- | :--- | :--- |
| el conteo $A$ | `statistic` | `count` |
| el veredicto | `trend_free` | `stationary` |
| región de aceptación | `bounds` (semiabierta) | `bounds` (semiabierta) |
| valor p exacto | `p_value` | `p_value` |
| la secuencia | `values` | `segment_values`, `segment_times`, `n_segments` |
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.
**¿Cómo de pequeña es la deriva que llega a ver?** Cada media cuadrática de
segmento es ella misma una variable aleatoria: para una banda de anchura $B$
observada durante $T_\text{seg}$ segundos, su error aleatorio normalizado es
$1/\sqrt{B\,T_\text{seg}}$. Toma los 20 segmentos por defecto sobre un registro
de 10 s de la banda de tercio de octava de 1 kHz ($B = 231$ Hz,
$T_\text{seg} = 0{,}5$ s): $BT = 116$, así que cada valor de segmento se dispersa
alrededor de un 9 %, que son 0,4 dB. Una deriva mucho menor que eso no se puede
separar de la dispersión, mientras que una deriva monótona de un par de
desviaciones típicas de segmento a lo largo del registro se rechaza casi siempre.
De ahí salen dos consecuencias de diseño, y apuntan en sentido contrario a la
maniobra evidente. Para ver una deriva *más pequeña*, alarga el **registro**, no
el número de segmentos: la dispersión cae como $1/\sqrt{T_\text{seg}}$ mientras
la tendencia se queda donde está, así que trocear el mismo registro en más
segmentos hace más ruidoso cada uno. Y mantén $T_\text{seg}$ muy por encima de
$1/B$ (al menos $BT = 5$ por intervalo, y no menos de diez ciclos de la
frecuencia más baja de interés), o los valores de segmento dejan de ser
independientes y la región de aceptación tabulada ya no vale. Un test de 20
segmentos sobre un registro limitado en banda a 20 Hz y por encima necesita, por
tanto, intervalos de al menos 0,5 s, y con ello un registro de al menos 10 s.
Hay una forma de no estacionariedad que se escapa de este test por construcción:
una fuente que alterna entre dos regímenes estables, una máquina que carga y
descarga, tráfico en pelotones, produce una secuencia de segmentos no monótona
que el test de inversiones puede aceptar tan tranquilo mientras el registro es a
todas luces no estacionario. Para eso está `method="runs"` en el §3, y la
respuesta honesta suele ser trocear el registro por régimen (`segment_times`
enseña dónde) y cualificar cada parte por separado.
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 pasa `statistic="mean"` o prueba 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.
La primera salvedad es la que le cuesta un resultado equivocado a quien lee,
porque hace que el test *acepte* un registro que debería rechazar. Merece la pena
verla:
*Un barrido de amplitud constante de 200 Hz a 2 kHz es todo lo no estacionario
que puede ser un registro, y la media cuadrática de banda completa no lo ve: los
veinte valores de segmento van de 0,4998 a 0,5002 y el conteo $A = 96$ se queda
cómodamente dentro de $(64, 125]$. Limita antes el mismo registro a 200-400 Hz y
los valores de segmento se desploman en cuanto el barrido abandona la banda,
dando $A = 173$, por encima de la cota superior, y rechazado. El remedio de la
salvedad de arriba no es una formalidad: lo que el test ve depende por completo
del estadístico y de la banda que se le entreguen.*
Mostrar el código de esta figura
```python
from scipy import signal
# stationarity_test es el que se importó en el §2 de arriba.
fs = 8192.0
t = np.arange(1 << 16) / fs
glide = np.sin(2 * np.pi * (200 * t + (2000 - 200) / (2 * t[-1]) * t ** 2))
sos = signal.butter(6, [200 / (fs / 2), 400 / (fs / 2)], btype="band",
output="sos")
full = stationarity_test(glide, fs)
banded = stationarity_test(signal.sosfiltfilt(sos, glide), fs)
print(full.count, full.stationary) # 96 True -- ciego
print(banded.count, banded.stationary) # 173 False -- lo pilla
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 7))
for ax, res, label in ((ax1, full, "banda completa"),
(ax2, banded, "200-400 Hz")):
ax.plot(np.arange(1, res.n_segments + 1), res.segment_values, "o-")
ax.set_ylim(-0.03, 0.58)
ax.set_title(f"{label}: A = {res.count}, aceptación {res.bounds}")
ax.set_xlabel("Índice de segmento")
ax.set_ylabel("Media cuadrática por segmento")
plt.show()
```
## 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="Por encima de la mediana")
ax.scatter(idx[~above], seq[~above], color="r",
label="Por debajo de 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: una sinusoide 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_\mathrm{c}$ da $N_0 = 2\sqrt{f_\mathrm{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/signals/spectra/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) # una sinusoide 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")
```
*Dos momentos del espectro predicen una curva que abarca casi tres décadas de
tasa. La tasa de cruces por cero medida en este registro de 800-1200 Hz es de
2012,7 cruces por segundo frente a la forma cerrada del Ejemplo 5.13,
$2\sqrt{f_\mathrm{c}^2 + B^2/12} = 2013{,}3$, una frecuencia aparente de 1007 Hz; a
partir de ahí la exponencial de Rice $N_a = N_0 e^{-a^2/2\sigma_x^2}$ sigue los
recuentos con menos de un 2 % de diferencia hasta $2\sigma$, donde la tasa ya ha
caído a 268/s. Pasado $3\sigma$ el recuento es de un puñado de sucesos por
segundo y la dispersión es la estadística del registro, no un fallo del modelo.*
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 de banda limitada:
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.
*Izquierda: el aspecto real del cribado. Tres registros con la misma varianza y
la misma banda de 800-1200 Hz, cada uno dibujado con su propia curva de Rice para
que la desviación sea un desajuste del modelo y no un artefacto de escala. La
referencia gaussiana cae sobre la curva en todos los niveles. El recorte duro a
2,5σ deja intactas las tasas por debajo del nivel de recorte y las elimina por
completo por encima: la firma es un acantilado, no una pendiente. Los picos
dispersos de 6σ levantan las dos colas por encima de la curva a partir de unos
2,5σ sin apenas mover el centro. Derecha: el aviso complementario del §5, medido.
Un registro que es ruido paso bajo ideal hasta 2 kHz más un suelo plano 50 dB por
debajo mantiene $r$ en su valor de forma cerrada 0,745 mientras la banda de
análisis se mantenga cerca de la física, y lo pierde conforme la banda se abre:
0,74 a 6 kHz, 0,71 a 13 kHz y 0,47 a 22 kHz. Un suelo de ruido que nunca se
notaría en una gráfica de espectro cuesta un tercio del factor de irregularidad
en cuanto la banda de análisis es diez veces más ancha de la cuenta.*
Por último, la mitad de baja frecuencia del mismo argumento, con la que abrió el
§4 y a la que no volvió: **las fórmulas de Rice están escritas para un proceso de
media nula**. Resta la media del registro y filtra en paso alto todo lo que quede
por debajo de la banda con significado físico antes de contar: 20 Hz para un
registro acústico de banda de audio, la propia frecuencia de corte inferior del
transductor para un registro de vibración. El viento y el retumbo del edificio
están casi por entero en $m_0$ y no aportan casi nada a $m_2$, así que inflan
$\sigma_x$ dejando $\sigma_v$ intacta y sesgan a la baja $N_0$ y todas las tasas
de cruce por nivel; una deriva lenta de continua puede suprimir por completo los
cruces de los niveles bajos, y un registro que nunca vuelve a cero no tiene
ningún cruce por cero. La regla es simétrica: limita la banda por los dos
extremos, el bajo por $m_0$ y el alto por $m_4$, y declara la banda que usaste
siempre que informes de una frecuencia aparente o de un factor de irregularidad,
porque ninguno de los dos números significa nada sin ella.
## 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)
```
*Un solo número decide qué curva siguen los picos. El ruido paso bajo ideal
tiene la forma cerrada $r = \sqrt5/3 = 0{,}745$, y el registro mide 0,746, que
es también por lo que el 12,6 % de sus 39 534 máximos son **negativos**, algo
imposible bajo Rayleigh. La escalera empírica se apoya en la mezcla de Rice para
ese $r$ y queda en todas partes entre los dos límites: la probabilidad de un
pico más allá de $4\sigma$ es $2{,}5\times10^{-4}$, una cuarta parte por debajo
del $e^{-8} = 3{,}4\times10^{-4}$ de Rayleigh. Suponer el límite de Rayleigh en
un registro que no es de banda estrecha sobrestima los extremos.*
Mostrar el código de esta figura
```python
from phonometry import peak_statistics
from phonometry.metrology.data_qualification 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 - limita antes el registro a la
banda físicamente significativa.
## Dónde encaja esto
La cualificació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/signals/spectra/spectral-analysis/) y las
fórmulas de error
aleatorio de los [estimadores de
correlación](/phonometry/es/signals/spectra/correlation-delay/)
suponen que el registro es estacionario, igual que la idea misma de *el*
$L_\mathrm{eq}$ de una medición en [Niveles](/phonometry/es/signals/levels/levels/). Cuando un
registro no pasa el test, divídelo en el cambio (la secuencia
`segment_values` muestra dónde), analiza los trozos, o pasa a las
vistas de tiempo corto - el
[espectrograma calibrado](/phonometry/es/signals/spectra/time-frequency/) -
que no suponen
estacionariedad. Y cuando pasa, la
[maquinaria GUM](/phonometry/es/signals/metrology/gum-uncertainty/) puede
propagar el error
aleatorio *restante* de cada estimación promediada con la conciencia
tranquila.
## Qué cubre esta guía
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).
La sección 10.3 del libro se titula «cualificació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
- [Análisis espectral](/phonometry/es/signals/spectra/spectral-analysis/): el intervalo de confianza chi-cuadrado que da por supuesta la estacionariedad que esta página comprueba.
- [Análisis tiempo-frecuencia](/phonometry/es/signals/spectra/time-frequency/): la vista de corto plazo para los registros que suspenden el test, y para los barridos que su media cuadrática no ve.
- [Incertidumbre GUM](/phonometry/es/signals/metrology/gum-uncertainty/): propagar el error aleatorio que le queda a un promedio ya cualificado.
- [Niveles](/phonometry/es/signals/levels/levels/): el $L_\mathrm{eq}$ cuya definición supone una única medición estacionaria.
- Referencia de la API: [`metrology.data_qualification`](/phonometry/es/reference/api/metrology/data-qualification/).
- Teoría: [Incertidumbre de medida (GUM)](/phonometry/es/reference/theory/signal-analysis/#incertidumbre-de-medida-isoiec-guide-98-3-gum-y-suplemento-1): el marco de incertidumbre que alimentan los criterios de cualificación de esta página.