Cualificación de datos: estacionariedad y picos

Referencias: Bendat y Piersol 2010 Wald y Wolfowitz 1940 Rice 1945

Todo promedio de esta documentación - una PSD, un Leq, un balance GUM - 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.

Diagrama de flujo de decisión de la cualificación de datos: un registro temporal, antes de confiar en cualquier promedio PSD, Leq o GUM, se divide en 20 segmentos iguales cuyas medias cuadráticas forman una secuencia, el recuento de inversiones A con media sin tendencia 95 se compara con la región de aceptación de la Tabla A.6, mayor que 64 y como máximo 125, al nivel del 5 por ciento, y el rombo se bifurca hacia una caja verde de estacionario, ruido estable con A igual a 91 aceptado de modo que valen los intervalos chi-cuadrado y las fórmulas de error, o hacia una caja roja de no estacionario, una rampa de ganancia del 20 por ciento con A igual a 7 rechazada, aconsejando dividir en el cambio o pasar a análisis de corto plazo; las notas al pie nombran el test de rachas como compañero y el punto ciego de la media cuadrática ante deslizamientos de frecuencia

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:

import numpy as np

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	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 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))

μ_{A} = \frac{N ( N - 1 )}{4}, σ_{A}^{2} = \frac{N ( 2 N + 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 $α$ cuando $A$ cae dentro de una región bilateral - la Tabla A.6 de B&P, cuyas filas de $α = 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:

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")

Dos secuencias de veinte observaciones dibujadas frente al índice de muestra: el Ejemplo 4.4 de Bendat y Piersol, que fluctúa en torno a cinco con ochenta y seis inversiones de orden y se acepta como sin tendencia, y las mismas fluctuaciones con una deriva ascendente añadida que sube de unos cinco a nueve, cuyas treinta y ocho inversiones de orden caen por debajo del límite inferior de aceptación de sesenta y cuatro y se rechazan como con tendencia al nivel del cinco por ciento

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

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)   # 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 $α = 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:

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 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")

Veinte medias cuadráticas por segmento de dos registros de ruido: un registro estable cuyos valores fluctúan en torno a uno con un conteo de inversiones de orden de noventa y uno, aceptado como estacionario, y el mismo ruido con una rampa de ganancia del veinte por ciento cuyas medias cuadráticas por segmento suben de forma sostenida hasta uno coma cinco, dando solo siete inversiones de orden, rechazado como no estacionario al nivel del cinco por ciento

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

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"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_{seg}$ segundos, su error aleatorio normalizado es $1/ B T_{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_{seg} = 0, 5$ s): $B T = 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/ T_{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_{seg}$ muy por encima de $1/ B$ (al menos $B T = 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:

Tres paneles apilados. Arriba: la frecuencia instantánea de un barrido de amplitud constante que sube linealmente de 200 Hz a 2 kHz en ocho segundos, con la banda de 200 a 400 Hz sombreada. En medio: las veinte medias cuadráticas de segmento del registro de banda completa, planas en 0,5 con un recuento de inversiones de 96 dentro de la región de aceptación, aceptado. Abajo: el mismo registro limitado a la banda de 200 a 400 Hz, cuyas medias cuadráticas de segmento se desploman a cero a partir del tercer segmento, dando un recuento de 173 por encima de la cota superior, rechazado

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

import matplotlib.pyplot as plt
import numpy as np
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 $α = 0, 05$ es la clásica $(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(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

Dos paneles de secuencias clasificadas respecto a su mediana, puntos por encima en azul y por debajo en rojo. Izquierda: cuarenta observaciones gaussianas con veinte rachas dentro de la región de aceptación de catorce a veintisiete, veredicto sin tendencia. Derecha: una secuencia de veinte valores estrictamente alternante que también da veinte rachas, pero contra la región de aceptación de seis a quince para veinte valores, veredicto rechazado por alternancia demasiado rápida

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

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)]

# 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) 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}},

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} = 2 B / 3$ (una frecuencia aparente de $0, 58 B$ , Ejemplo 5.12) y una banda centrada en $f_{c}$ da $N_{0} = 2 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 del registro:

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)          # 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")

Tasas de cruce por nivel medidas de un registro de ruido gaussiano de banda limitada entre ochocientos y mil doscientos hercios, dibujadas como puntos frente al nivel de cruce desde menos tres coma cinco hasta más tres coma cinco unidades de la señal en un eje de tasas logarítmico, cayendo desde unos dos mil cruces por segundo en el nivel cero hasta unos pocos por segundo a tres sigmas, con la curva exponencial de Rice pasando por todos los puntos medidos

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 f_{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 σ_{x}^{2}}$ sigue los recuentos con menos de un 2 % de diferencia hasta $2 σ$ , donde la tasa ya ha caído a 268/s. Pasado $3 σ$ 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

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)        # 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: tasas de cruce por nivel de tres registros de la misma varianza y el mismo ancho de banda frente al nivel de cruce en unidades de sigma, cada uno con su propia curva de Rice. La referencia gaussiana cae sobre su curva, el registro recortado duro en 2,5 sigma no tiene ningún cruce más allá del nivel de recorte, y el registro con picos dispersos de seis sigma queda por encima de su curva en las dos colas. Derecha: el factor de irregularidad de un registro paso bajo de dos kilohercios con un suelo de ruido 50 decibelios por debajo, representado frente a la frecuencia de corte superior de análisis, plano en el valor de forma cerrada 0,745 cerca de la banda física y cayendo a 0,47 en 22 kilohercios

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 $σ_{x}$ dejando $σ_{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 π) (σ_{a} / σ_{v}) = m_{4} / m_{2}$ (Ec. (5.211)) y con ella el factor de irregularidad adimensional

r = \frac{N _{0}}{2 M} = \frac{m _{2}}{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: $Prob [pico > a] = e^{- a^{2} /2 σ_{x}^{2}}$ (Ec. (5.206)), la probabilidad de uno entre 3000 de un pico más allá de $4 σ$ 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:

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)          # 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)

$Probabilidad de excedencia de la altura de pico de ruido gaussiano paso bajo en eje logarítmico frente a la altura de pico estandarizada desde menos dos coma cinco hasta cuatro coma cinco: la escalera empírica de medio millón de muestras sigue la curva de mezcla de Rice para un factor de irregularidad de cero coma siete cuatro seis, claramente por debajo del límite de Rayleigh discontinuo y por encima del límite gaussiano punteado, con una fracción visible de máximos negativos$

Un solo número decide qué curva siguen los picos. El ruido paso bajo ideal tiene la forma cerrada $r = 5 /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 σ$ es $2, 5 \times 1 0^{- 4}$ , una cuarta parte por debajo del $e^{- 8} = 3, 4 \times 1 0^{- 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

import matplotlib.pyplot as plt
import numpy as np
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 = 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 y las fórmulas de error aleatorio de los estimadores de correlación suponen que el registro es estacionario, igual que la idea misma de el $L_{eq}$ de una medición en Niveles. 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 - que no suponen estacionariedad. Y cuando pasa, la maquinaria GUM 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 «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: el intervalo de confianza chi-cuadrado que da por supuesta la estacionariedad que esta página comprueba.
Análisis tiempo-frecuencia: 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: propagar el error aleatorio que le queda a un promedio ya cualificado.
Niveles: el $L_{eq}$ cuya definición supone una única medición estacionaria.
Referencia de la API: metrology.data_qualification.
Teoría: Incertidumbre de medida (GUM): el marco de incertidumbre que alimentan los criterios de cualificación de esta página.

Referencias

Bendat, J. S. y Piersol, A. G. (2010). Random data: Analysis and measurement procedures (4.ª ed.). Wiley. https://doi.org/10.1002/9781118032428Sección 4.5.2 con la Tabla A.6 (el test no paramétrico de tendencia por inversiones de orden, su media y varianza, los puntos porcentuales tabulados y el Ejemplo 4.4), Sección 10.3 (cualificación de datos: clasificación, validación, edición; el procedimiento de estacionariedad por medias cuadráticas de segmento de 10.3.1.1 y el Ejemplo 10.3) y Sección 5.5 (cruces por nivel y valores de pico: Ecs. (5.195)-(5.196), (5.206), (5.211), (5.217)-(5.223) y Ejemplos 5.12-5.15). ISBN 978-0-470-24877-5. El test de rachas complementario aparecía en los puntos porcentuales de la distribución de rachas de la tercera edición.
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.xLas derivaciones originales de las tasas esperadas de cruces por nivel y de máximos y de la distribución de picos del ruido gaussiano que presenta la Sección 5.5 de Bendat y Piersol (su Ref. 6; las Partes I-II están en el volumen 23, 1944).
Wald, A. y 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/1177731909La distribución condicional exacta del número de rachas, de la que se calculan las regiones de aceptación de rachas respecto de la mediana.

Creado y mantenido por José Manuel Requena Plens ORCID Todos los proyectos