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 las lecturas de temperatura horaria, 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 preservan los bordes utilizando un filtro mediano. En el ejemplo también se muestra cómo utilizar un filtro Hampel para eliminar los valores atípicos grandes.

Motivación

Suavizar es cómo descubrimos patrones importantes en nuestros datos mientras dejamos fuera cosas que no son importantes (es decir, ruido). Utilizamos el filtrado para realizar este alisado. El objetivo de suavizar es producir cambios lentos de valor para que sea más fácil ver las tendencias en nuestros datos.

A veces, cuando se examinan los datos de entrada es posible que desee suavizar los datos con el fin de 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, 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)')

Observe que podemos ver visualmente el efecto que la hora del día tiene sobre las lecturas de la temperatura. Si sólo está interesado en la variación diaria de la temperatura durante el mes, las fluctuaciones por hora sólo contribuyen con el ruido, lo que puede hacer que las variaciones diarias difíciles de discernir. Para eliminar el efecto de la hora del día, ahora quisiéramos alisar nuestros datos usando un filtro medio móvil.

Un filtro medio móvil

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

Para aplicar un filtro medio móvil a cada punto de datos, construimos nuestros coeficientes de nuestro filtro de modo que cada punto sea igualmente ponderado y contribuya 1/24 al promedio 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 aproximadamente doce horas. Esto se debe al hecho de que nuestro filtro medio móvil tiene un retraso.

Cualquier filtro simétrico de longitud n tendrá un retardo 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 la temperatura horaria. Luego, segmenta los datos diferenciados en días y toma 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)')

Extrayendo el sobre máximo

A veces también quisiéramos tener una estimación suavemente diversa de cómo los altos y los mínimos de nuestra señal de la temperatura cambian diariamente. Para ello, podemos utilizar la función envelope para conectar las alturas y mínimos extremos detectados en un subconjunto del período de 24 horas. En este ejemplo, nos aseguramos de que haya por lo menos 16 horas entre cada extremo alto y extremo bajo. También podemos tener una idea de cómo las altas y bajas están en la tendencia al tomar 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 medios móviles ponderados

Otras clases de filtros medios móviles no pesan cada muestra igualmente.

Otro filtro común sigue la expansión binomial de . Este tipo de filtro aproxima una curva normal para valores grandes de n. Es útil para filtrar el ruido de alta frecuencia para la pequeña n. Para encontrar los coeficientes para el filtro binomial, convolve con sí mismo y luego iterativamente convolverá la salida con un número de veces prescrito. 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 de Gauss es el filtro medio móvil exponencial. Este tipo de filtro medio móvil ponderado 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 mayor 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 por un día.

axis([3 4 -5 2])

Savitzky-filtros Golay

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

Para rastrear 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 conveniencia, puede utilizar la función sgolayfilt para implementar un filtro de suavizado Savitzky-Golay. Para utilizar sgolayfilt, se especifica un segmento de longitud impar de los datos y un orden polinómico estrictamente inferior a la longitud del segmento. La función sgolayfilt 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.

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, muestreamos el voltaje de bucle abierto a través de la entrada de un instrumento analógico en presencia de interferencias de 60 Hz de ruido de línea de alimentación de CA. Muestreamos 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')

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

Si se construye un filtro medio móvil ponderado uniformemente, se 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. Vamos a intentar "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 si bien 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 capturar un ciclo completo de la señal de 60 Hz por nuestro filtro medio móvil.

Si remuestreamos la señal a 17 * 60 Hz = 1020 Hz, podemos utilizar nuestro filtro medio 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 busca. 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 los filtros de Savitzky-Golay respectivamente debajo-corrigen y sobre-corrigen cerca de los bordes de la señal del reloj.

Una manera simple de preservar los bordes, pero aún 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 mediante 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 añadir algunos picos de ruido artificial:

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

Dado que cada espiga introducida tiene una duración de sólo una muestra, podemos utilizar un filtro mediano de sólo tres elementos para eliminar las espigas.

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

El filtro eliminó las espigas, 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 sólo los valores que equivalen a algunas desviaciones estándar alejadas del valor medio 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 filtrado y remuestreo, consulte la caja 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. 4th Ed. London: Macmillan, 1983.