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ñales

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 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 mientras se conservan los bordes mediante el uso de un filtro mediano. En el ejemplo también se muestra cómo utilizar un filtro Hampel para eliminar valores atípicos grandes.

Motivación

El suavizado es la forma en que descubrimos patrones importantes en nuestros datos mientras omitimos cosas que no son importantes (es decir, ruido). Usamos 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 durante todo el mes de enero de 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 en las lecturas de temperatura. Si solo está interesado en la variación de temperatura diaria durante el mes, las fluctuaciones por hora sólo contribuyen al ruido, lo que puede hacer que las variaciones diarias sean difíciles de discernir. 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 forma de 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 contribuya con 1/24 a la media total. Esto nos da la temperatura media 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 del filtro

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

Cualquier filtro simétrico de longitud N tendrá un retraso de (N-1)/2 muestras. Podemos explicar 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, en primer lugar, reste los datos suavizados de las mediciones de temperatura por hora. A continuación, segmente los datos diferenciados en días y tome el promedio durante 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 de sobre pico

A veces también nos gustaría tener una estimación suave mente variable de cómo los máximos y bajos 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 altículos y bajos 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 ponderados

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ños. Para encontrar los coeficientes para el filtro binomial, convolve
<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>
con sí mismo y luego iterativamente convolve 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.

Se ajusta un filtro de media móvil ponderado exponencialmente mediante un parámetro alfa entre cero y uno. Un valor más alto de alfa 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)')

Amplíe las lecturas durante un día.

axis([3 4 -5 2])

Filtros Savitzky-Golay

Observará que al suavizar los datos, los valores extremos se recortaron un poco.

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

Como comodidad, puede utilizar la función para implementar un filtro de suavizado Savitzky-Golay.sgolayfilt Para utilizar , especifique un segmento de longitud impar de los datos y un orden polinómial 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 principio 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 remuestrear una señal para aplicar correctamente una media móvil.

En nuestro siguiente ejemplo, hemos muestreado el voltaje de bucle abierto a través de la entrada de un instrumento analógico en presencia de interferencia de 60 Hz de ruido de línea de alimentación de CA. Hemos muestreado la tensión 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')

Vamos a intentar eliminar el efecto del ruido de línea mediante el uso de un filtro de media móvil.

Si construye un filtro de media móvil ponderado uniformemente, 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 muestrea 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 a 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 suaviza 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 completo de la señal de 60 Hz por nuestro filtro de media móvil.

Si remuestreamos la señal a 17 * 60 Hz a 1020 Hz, podemos utilizar nuestro filtro de media móvil de 17 puntos para eliminar el ruido de 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 tienen bordes afilados que no queremos 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 los filtros Savitzky-Golay respectivamente sub-correctos y sobre-correctos cerca de los bordes de la señal del reloj.

Una forma sencilla de conservar los bordes, pero 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 valores atípicos a través del filtro Hampel

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 artificial:

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

Dado que cada pico que introdujimos tiene una duración de una sola muestra, podemos usar un filtro mediano de solo 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 Hampel funciona de forma similar a un filtro mediano, sin embargo, reemplaza solo los valores que son equivalentes a unas pocas desviaciones estándar lejos del valor mediano local.

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

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

Lectura adicional

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

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