Esta página aún no se ha traducido para esta versión. Puede ver la versión más reciente de esta página en inglés.

Suavizado de señal

Este ejemplo muestra cómo utilizar filtros de media móvil y remuestreo para aislar el efecto de los componentes periódicos de la hora del día en las lecturas de temperatura por hora, así como eliminar el ruido de línea no deseado de una medición de voltaje de bucle abierto. El ejemplo también muestra cómo suavizar los niveles de una señal de reloj preservando los bordes mediante un filtro mediano. El ejemplo también muestra cómo utilizar un filtro de Hampel para eliminar los valores atípicos grandes.

Motivación

El suavizado es la forma en que descubrimos patrones importantes en nuestros datos, dejando fuera cosas que no son importantes (es decir, ruido). Utilizamos el filtrado para realizar este suavizado. El objetivo del suavizado es producir cambios lentos en el valor para que sea más fácil ver las tendencias en nuestros datos.

A veces, cuando examina los datos de entrada, es posible que desee suavizar los datos para ver una tendencia en la señal. En nuestro ejemplo tenemos un conjunto de lecturas de temperatura en Celsius tomadas cada hora en el aeropuerto Logan para todo el mes de enero, 2011.

load bostemp days = (1:31*24)/24; plot(days, tempC) axis tight ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Tenga en cuenta que podemos ver visualmente el efecto que la hora del día tiene sobre las lecturas de temperatura. Si solo está interesado en la variación diaria de temperatura durante el mes, las fluctuaciones horarias solo contribuyen al ruido, lo que puede dificultar la discernir las variaciones diarias. Para eliminar el efecto de la hora del día, ahora nos gustaría suavizar nuestros datos mediante el uso de un filtro de media móvil.

Un filtro de media móvil

En su forma más simple, un filtro de media móvil de longitud N toma el promedio de cada N muestras consecutivas de la onda.

Para aplicar un filtro de media móvil a cada punto de datos, construimos nuestros coeficientes de nuestro filtro para que cada punto esté igualmente ponderado y aporte 1/24 al promedio total. Esto nos da la temperatura promedio sobre cada período de 24 horas.

hoursPerDay = 24; coeff24hMA = ones(1, hoursPerDay)/hoursPerDay;  avg24hTempC = filter(coeff24hMA, 1, tempC); plot(days,[tempC avg24hTempC]) legend('Hourly Temp','24 Hour Average (delayed)','location','best') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Retardo de filtro

Tenga en cuenta que la salida filtrada se retrasa aproximadamente doce horas. Esto es debido al hecho de que nuestro filtro de media móvil tiene un retraso.

Cualquier filtro simétrico de longitud N tendrá un retardo de (N-1)/2 muestras. Podemos tener en cuenta este retraso manualmente.

fDelay = (length(coeff24hMA)-1)/2; plot(days,tempC, ...      days-fDelay/24,avg24hTempC) axis tight legend('Hourly Temp','24 Hour Average','location','best') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Extracción de diferencias medias

Alternativamente, también podemos utilizar el filtro de media móvil para obtener una mejor estimación de cómo la hora del día afecta a la temperatura general. Para ello, primero, reste los datos suavizados de las mediciones de temperatura por hora. A continuación, segmenta los datos diferenciados en días y toma el promedio de los 31 días del mes.

figure deltaTempC = tempC - avg24hTempC; deltaTempC = reshape(deltaTempC, 24, 31).';  plot(1:24, mean(deltaTempC)) axis tight title('Mean temperature differential from 24 hour average') xlabel('Hour of day (since midnight)') ylabel('Temperature difference (\circC)')

Extracción del sobre pico

A veces también nos gustaría tener una estimación suavemente variable de cómo los máximos y mínimos de nuestra señal de temperatura cambian diariamente. Para ello podemos utilizar la función para conectar máximos y mínimos extremos detectados durante un subconjunto del período de 24 horas.envelope En este ejemplo, nos aseguramos de que haya al menos 16 horas entre cada extremo alto y extremo bajo. También podemos tener una idea de cómo los máximos y mínimos están en tendencia tomando el promedio entre los dos extremos.

[envHigh, envLow] = envelope(tempC,16,'peak'); envMean = (envHigh+envLow)/2;  plot(days,tempC, ...      days,envHigh, ...      days,envMean, ...      days,envLow)     axis tight legend('Hourly Temp','High','Mean','Low','location','best') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Filtros de media móvil ponderada

Otros tipos de filtros de media móvil no ponderan cada muestra por igual.

Otro filtro común sigue la expansión binomial de

<math display="inline">
<mrow>
<msup>
<mrow>
<mrow>
<mo>[</mo>
<mrow>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
<mo stretchy="false">,</mo>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
</mrow>
<mo>]</mo>
</mrow>
</mrow>
<mrow>
<mi mathvariant="italic">n</mi>
</mrow>
</msup>
</mrow>
</math>
. Este tipo de filtro se aproxima a una curva normal para valores grandes de n. Es útil para filtrar el ruido de alta frecuencia para n pequeño. Para encontrar los coeficientes para el filtro binomial, herramienta
<math display="inline">
<mrow>
<mrow>
<mo>[</mo>
<mrow>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
<mo stretchy="false">,</mo>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
</mrow>
<mo>]</mo>
</mrow>
</mrow>
</math>
consigo mismo y luego iterativamente herramienta la salida con
<math display="inline">
<mrow>
<mrow>
<mo>[</mo>
<mrow>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
<mo stretchy="false">,</mo>
<mrow>
<mrow>
<mn>1</mn>
</mrow>
<mo stretchy="false">/</mo>
<mrow>
<mn>2</mn>
</mrow>
</mrow>
</mrow>
<mo>]</mo>
</mrow>
</mrow>
</math>
un número prescrito de veces. En este ejemplo, utilice cinco iteraciones totales.

h = [1/2 1/2]; binomialCoeff = conv(h,h); for n = 1:4     binomialCoeff = conv(binomialCoeff,h); end  figure fDelay = (length(binomialCoeff)-1)/2; binomialMA = filter(binomialCoeff, 1, tempC); plot(days,tempC, ...      days-fDelay/24,binomialMA) axis tight legend('Hourly Temp','Binomial Weighted Average','location','best') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Otro filtro algo similar al filtro de expansión gaussiana es el filtro de media móvil exponencial. Este tipo de filtro de media móvil ponderada es fácil de construir y no requiere un tamaño de ventana grande.

Ajustar un filtro de media móvil exponencialmente ponderado por un parámetro alfa entre cero y uno. Un valor más alto de Alpha tendrá menos suavizado.

alpha = 0.45; exponentialMA = filter(alpha, [1 alpha-1], tempC); plot(days,tempC, ...      days-fDelay/24,binomialMA, ...      days-1/24,exponentialMA)  axis tight legend('Hourly Temp', ...        'Binomial Weighted Average', ...        'Exponential Weighted Average','location','best') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)')

Acerque las lecturas durante un día.

axis([3 4 -5 2])

Filtros Savitzky-Golay

Tendrá en cuenta que al suavizar los datos, los valores extremos se recortan un poco.

Para rastrear la señal un poco más de cerca, puede utilizar un filtro de media móvil ponderada que intente ajustar un polinomio de un orden especificado sobre un número especificado de muestras en un sentido de mínimos cuadrados.

Como conveniencia, puede utilizar la función para implementar un filtro de suavizado Savitzky-Golay.sgolayfilt Para utilizar, se especifica un segmento de longitud impar de los datos y un orden polinómico estrictamente inferior a la longitud del segmento.sgolayfilt La función calcula internamente los coeficientes polinómicos de suavizado, realiza la alineación de retardo y se encarga de los efectos transitorios al inicio y al final del registro de datos.sgolayfilt

cubicMA   = sgolayfilt(tempC, 3, 7); quarticMA = sgolayfilt(tempC, 4, 7); quinticMA = sgolayfilt(tempC, 5, 9); plot(days,[tempC cubicMA quarticMA quinticMA]) legend('Hourly Temp','Cubic-Weighted MA', 'Quartic-Weighted MA', ...        'Quintic-Weighted MA','location','southeast') ylabel('Temp (\circC)') xlabel('Time elapsed from Jan 1, 2011 (days)') title('Logan Airport Dry Bulb Temperature (source: NOAA)') axis([3 5 -5 2])

Remuestreo

A veces es beneficioso para remuestrear una señal con el fin de aplicar correctamente una media móvil.

En nuestro siguiente ejemplo, probamos la tensión de bucle abierto a través de la entrada de un instrumento analógico en presencia de interferencias de ruido de la línea de alimentación de CA de 60 Hz. Probamos el voltaje con una frecuencia de muestreo de 1 kHz.

load openloop60hertz fs = 1000; t = (0:numel(openLoopVoltage)-1) / fs; plot(t,openLoopVoltage) ylabel('Voltage (V)') xlabel('Time (s)') title('Open-loop Voltage Measurement')

Intentemos eliminar el efecto del ruido de la línea utilizando un filtro de media móvil.

Si construye un filtro de media móvil uniformemente ponderado, eliminará cualquier componente que sea periódico con respecto a la duración del filtro.

Hay aproximadamente 1000/60 = 16,667 muestras en un ciclo completo de 60 Hz cuando se muestrean a 1000 Hz. Intentemos "redondear" y usar un filtro de 17 puntos. Esto nos dará un filtrado máximo a una frecuencia fundamental de 1000 Hz/17 = 58,82 Hz.

plot(t,sgolayfilt(openLoopVoltage,1,17)) ylabel('Voltage (V)') xlabel('Time (s)') title('Open-loop Voltage Measurement') legend('Moving average filter operating at 58.82 Hz', ...        'Location','southeast')

Tenga en cuenta que mientras que el voltaje se alisa significativamente, todavía contiene una pequeña ondulación de 60 Hz.

Podemos reducir significativamente la ondulación si remuestreamos la señal para que capturemos un ciclo completo completa de la señal de 60 Hz por nuestro filtro de media móvil.

Si remuestreamos la señal a 17 * 60 Hz = 1020 Hz, podemos utilizar nuestro filtro de 17 puntos de media móvil para eliminar el ruido de la línea de 60 Hz.

fsResamp = 1020; vResamp = resample(openLoopVoltage, fsResamp, fs); tResamp = (0:numel(vResamp)-1) / fsResamp; vAvgResamp = sgolayfilt(vResamp,1,17); plot(tResamp,vAvgResamp) ylabel('Voltage (V)') xlabel('Time (s)') title('Open-loop Voltage Measurement') legend('Moving average filter operating at 60 Hz', ...     'Location','southeast')

Filtro mediano

La media móvil, la media móvil ponderada y los filtros Savitzky-Golay suavizan todos los datos que filtran. Esto, sin embargo, puede no ser siempre lo que se quiere. Por ejemplo, ¿qué pasa si nuestros datos se toman de una señal de reloj y tiene bordes afilados que no deseamos suavizar? Los filtros discutidos hasta ahora no funcionan tan bien:

load clockex  yMovingAverage = conv(x,ones(5,1)/5,'same'); ySavitzkyGolay = sgolayfilt(x,3,5); plot(t,x, ...      t,yMovingAverage, ...      t,ySavitzkyGolay)  legend('original signal','moving average','Savitzky-Golay')

La media móvil y Savitzky-Golay filtros respectivamente bajo-corregir y sobre-corregir cerca de los bordes de la señal de reloj.

Una forma sencilla de preservar los bordes, pero aún así suavizar los niveles es utilizar un filtro mediano:

yMedFilt = medfilt1(x,5,'truncate'); plot(t,x, ...      t,yMedFilt) legend('original signal','median filter')

Eliminación de outlier a través de Hampel Filter

Muchos filtros son sensibles a los valores atípicos. Un filtro que está estrechamente relacionado con el filtro mediano es el filtro Hampel. Este filtro ayuda a eliminar los valores atípicos de una señal sin suavizar excesivamente los datos.

Para ver esto, cargue una grabación de audio de un silbato de tren y agregue algunos picos de ruido artificiales:

load train y(1:400:end) = 2.1; plot(y)

Dado que cada pico que introdujimos tiene una duración de sólo una muestra, podemos utilizar un filtro mediano de sólo tres elementos para eliminar los picos.

hold on plot(medfilt1(y,3)) hold off legend('original signal','filtered signal')

El filtro eliminó los picos, pero también eliminó un gran número de puntos de datos de la señal original. Un filtro de Hampel funciona de forma similar a un filtro mediano, sin embargo, reemplaza sólo los valores que son equivalentes a algunas desviaciones estándar lejos del valor de la mediana local.

hampel(y,13) legend('location','best')

Sólo se eliminan los valores atípicos de la señal original.

Leer más

Para obtener más información sobre el filtrado y el remuestreo, consulte la caja de herramientas de procesamiento de señal.

Referencia: Kendall, Maurice G., Alan Stuart, y J. Keith Ord. .The Advanced Theory of Statistics, Vol. 3: Design and Analysis, and Time-Series 4º Ed. Londres: Macmillan, 1983.