Elegir el filtro correcto para los datos del acelerómetro

28

Soy bastante nuevo en DSP y he realizado algunas investigaciones sobre posibles filtros para suavizar los datos del acelerómetro en Python. Un ejemplo del tipo de datos que estaré experimentando se puede ver en la siguiente imagen:

Datos del acelerómetro

Esencialmente, estoy buscando consejos para suavizar estos datos para eventualmente convertirlos en velocidad y desplazamiento. Entiendo que los acelerómetros de los teléfonos móviles son extremadamente ruidosos.

No creo que pueda usar un filtro de Kalman en este momento porque no puedo controlar el dispositivo para hacer referencia al ruido producido por los datos (¿leí que es esencial colocar el dispositivo plano y encontrar la cantidad de ruido de esas lecturas?)

FFT ha producido algunos resultados interesantes. Uno de mis intentos fue FFT la señal de aceleración, luego renderizar las frecuencias bajas para tener un valor FFT absoluto de 0. Luego utilicé la aritmética omega y la FFT inversa para obtener un gráfico de velocidad. Los resultados fueron los siguientes:

Señal filtrada

¿Es esta una buena manera de hacer las cosas? Estoy tratando de eliminar la naturaleza ruidosa general de la señal, pero es necesario identificar picos obvios, como alrededor de 80 segundos.

También me he cansado de usar un filtro de paso bajo en los datos del acelerómetro original, lo que ha hecho un gran trabajo alisándolo, pero no estoy realmente seguro de a dónde ir desde aquí. ¡Cualquier guía sobre dónde ir desde aquí sería realmente útil!

EDITAR: Un poco de código:

for i in range(len(fz)): 
    testing = (abs(Sz[i]))/Nz

    if fz[i] < 0.05:
        Sz[i]=0

Velfreq = []
Velfreqa = array(Velfreq)
Velfreqa = Sz/(2*pi*fz*1j)
Veltimed = ifft(Velfreqa)
real = Veltimed.real

Así que, esencialmente, he realizado una FFT en mis datos de acelerómetro, dando a Sz, filtrando altas frecuencias usando un simple filtro de pared de ladrillos (sé que no es lo ideal). Luego uso aritmética omega en la FFT de los datos. También muchas gracias a datageist por agregar mis imágenes a mi publicación :)

Michael M
fuente
¡Bienvenido a DSP! ¿Es la curva roja en su segunda imagen una versión "suavizada" de los datos originales (verdes)?
Phonon el
La curva roja es (¡ojalá!) Una curva de velocidad generada a partir de fft seguida de filtrado, seguida de aritmética omega (dividida por 2 * pi f j), seguida de inv. fft
Michael M
1
Quizás si incluye una expresión matemática más precisa o un seudocódigo para lo que hizo, aclararía un poco las cosas.
Phonon
Agregó algunos ahora, esa es la sensación general del código ..
Michael M
1
Mi pregunta sería: ¿qué esperas ver en los datos? No sabrá si tiene un buen enfoque a menos que sepa algo sobre la señal subyacente que espera ver después del filtrado. Además, el código que mostró es confuso. Aunque no muestra la inicialización de la fzmatriz, parece que está aplicando un filtro de paso alto en su lugar.
Jason R

Respuestas:

13

Como lo señaló @JohnRobertson en Bag of Tricks for Denoising Signals Mientras mantiene las transiciones nítidas, la eliminación de ruido de Variaton total (TV) es otra buena alternativa si su señal es constante por partes. Este puede ser el caso de los datos del acelerómetro, si su señal sigue variando entre diferentes mesetas.

μρ

yXμX-y2+reX1re

function denoise()

f = [-1*ones(1000,1);3*ones(100,1);1*ones(500,1);-2*ones(800,1);0*ones(900,1)];
plot(f);
axis([1 length(f) -4 4]);
title('Original');
g = f + .25*randn(length(f),1);
figure;
plot(g,'r');
title('Noisy');
axis([1 length(f) -4 4]);
fc = denoisetv(g,.5);
figure;
plot(fc,'g');
title('De-noised');
axis([1 length(f) -4 4]);

function f = denoisetv(g,mu)
I = length(g);
u = zeros(I,1);
y = zeros(I,1);
rho = 10;

eigD = abs(fftn([-1;1],[I 1])).^2;
for k=1:100
    f = real(ifft(fft(mu*g+rho*Dt(u)-Dt(y))./(mu+rho*eigD)));
    v = D(f)+(1/rho)*y;
    u = max(abs(v)-1/rho,0).*sign(v);
    y = y - rho*(u-D(f));
end

function y = D(x)
y = [diff(x);x(1)-x(end)];

function y = Dt(x)
y = [x(end)-x(1);-diff(x)];

Resultados:

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

Daniel R. Pipa
fuente
Realmente me gusta esta respuesta, voy a seguir adelante y probarlo. Lo siento, me tomó tanto tiempo responder.
Michael M
Excelente respuesta Gracias por los detalles Estoy buscando la versión C de este código. ¿Alguien ha portado este código matlab a C que les gustaría compartir? Gracias.
Rene Limberger el
¿Qué significa constante por piezas?
tilaprimera 01 de
6

El problema es que su ruido tiene un espectro plano. Si supone un ruido gaussiano blanco (que resulta ser una buena suposición), la densidad del espectro de potencia es constante. En términos generales, significa que su ruido contiene todas las frecuencias. Es por eso que cualquier enfoque de frecuencia, por ejemplo, DFT o filtros de paso bajo, no es bueno. ¿Cuáles serían sus frecuencias de corte ya que su ruido está en todo el espectro?

Una respuesta a esta pregunta es el filtro Wiener, que requiere el conocimiento de las estadísticas de su ruido y su señal deseada. Básicamente, la señal ruidosa (señal + ruido) se atenúa en las frecuencias en las que se espera que el ruido sea mayor que su señal, y se amplifica donde su señal se espera que sea mayor que su ruido.

Sin embargo, sugeriría enfoques más modernos que usan procesamiento no lineal, por ejemplo, eliminación de ondas. Estos métodos proporcionan excelentes resultados. Básicamente, la señal ruidosa se descompone primero en wavelets y luego se ponen a cero los coeficientes pequeños. Este enfoque funciona (y DFT no) debido a la naturaleza de múltiples resoluciones de las wavelets. Es decir, la señal se procesa por separado en las bandas de frecuencia definidas por la transformada wavelet.

En MATLAB, escriba 'wavemenu' y luego 'SWT denoising 1-D'. Luego 'Archivo', 'Análisis de ejemplo', 'Señales ruidosas', 'con Haar en el nivel 5, Bloques ruidosos'. Este ejemplo usa wavelet de Haar, que debería funcionar bien para su problema.

No soy bueno en Python, pero creo que puedes encontrar algunos paquetes de NumPy que realizan la eliminación de ruido de ondas Haar.

Daniel R. Pipa
fuente
1
No estaría de acuerdo con tu primera declaración. Está asumiendo que la señal de interés cubre el ancho de banda completo de la secuencia de entrada, lo cual es poco probable. Todavía es posible obtener una relación señal / ruido mejorada utilizando un filtro lineal en este caso, eliminando el ruido fuera de banda. Si la señal está muy sobremuestreada, puede obtener una gran mejora con un enfoque tan simple.
Jason R
Es cierto, y esto se logra mediante el filtro Wiener, cuando conoce las estadísticas de su señal y su ruido.
Daniel R. Pipa
Aunque la teoría detrás de la eliminación de ruido wavelet es complicada, la implementación es tan simple como el enfoque que describió. Implica solo bancos de filtros y umbrales.
Daniel R. Pipa el
1
Estoy investigando esto ahora, publicaré mi progreso arriba, ¡gracias a ustedes y a Phonon por toda su ayuda hasta ahora!
Michael M
@DanielPipa No tengo acceso a los paquetes de matlab en cuestión. ¿Puede proporcionar un documento u otra referencia que describa el método que corresponde a su código matlab.
John Robertson
0

Según la sugerencia de Daniel Pipa, eché un vistazo a la eliminación de ondas y encontré este excelente artículo de Francisco Blanco-Silva.

Aquí he modificado su código Python para que el procesamiento de imágenes funcione con datos 2D (acelerómetro) en lugar de datos 3D (imagen).

Tenga en cuenta que el umbral está "compuesto" para el umbral suave en el ejemplo de Francisco. Considere esto y modifique para su aplicación.

def wavelet_denoise(data, wavelet, noise_sigma):
    '''Filter accelerometer data using wavelet denoising

    Modification of F. Blanco-Silva's code at: https://goo.gl/gOQwy5
    '''
    import numpy
    import scipy
    import pywt

    wavelet = pywt.Wavelet(wavelet)
    levels  = min(15, (numpy.floor(numpy.log2(data.shape[0]))).astype(int))

    # Francisco's code used wavedec2 for image data
    wavelet_coeffs = pywt.wavedec(data, wavelet, level=levels)
    threshold = noise_sigma*numpy.sqrt(2*numpy.log2(data.size))

    new_wavelet_coeffs = map(lambda x: pywt.threshold(x, threshold, mode='soft'),
                             wavelet_coeffs)

    return pywt.waverec(list(new_wavelet_coeffs), wavelet)

Dónde:

  • wavelet- nombre de cadena de la forma wavelet que se utilizará (ver pywt.wavelist(), por ejemplo 'haar')
  • noise_sigma - desviación estándar del ruido de los datos
  • data - matriz de valores para filtrar (por ejemplo, datos de eje x, y o z)
ryanjdillon
fuente