Escolhendo o filtro correto para dados do acelerômetro

Eu sou bastante novo no DSP e fiz algumas pesquisas sobre possíveis filtros para suavizar dados do acelerômetro em python. Um exemplo do tipo de dados que estou enfrentando pode ser visto na imagem a seguir:

Essencialmente, estou procurando conselhos para suavizar esses dados e eventualmente convertê-los em velocidade e deslocamento. Entendo que os acelerômetros dos telefones celulares são extremamente barulhentos.

Eu não acho que posso usar um filtro Kalman no momento porque não consigo me apossar do dispositivo para fazer referência ao ruído produzido pelos dados (eu li que é essencial colocar o dispositivo na horizontal e encontrar a quantidade de ruído dessas leituras?)

A FFT produziu alguns resultados interessantes. Uma das minhas tentativas foi realizar a FFT do sinal de aceleração e, em seguida, renderizar as frequências baixas para ter um valor absoluto da FFT igual a 0. Então usei a FFT aritmética e inversa ômega para obter um gráfico de velocidade. Os resultados foram os seguintes:

Essa é uma boa maneira de fazer as coisas? Estou tentando remover a natureza geral barulhenta do sinal, mas picos óbvios, como em cerca de 80 segundos, precisam ser identificados.

Também me cansei de usar um filtro passa-baixo nos dados originais do acelerômetro, o que fez um ótimo trabalho em suavizá-lo, mas não tenho muita certeza de onde ir a partir daqui. Qualquer orientação sobre onde ir a partir daqui seria realmente útil!

EDIT: Um pouco 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

Então, essencialmente, eu fiz uma FFT nos dados do meu acelerômetro, fornecendo Sz, filtrando altas frequências usando um simples filtro de parede de tijolos (eu sei que não é o ideal). Em seguida, use aritmética ômega na FFT dos dados. Também muito obrigado ao datageist por adicionar minhas imagens no meu post :)

Conforme apontado por @JohnRobertson no Pacote de truques para diminuir o ruído de sinais durante a manutenção de transições nítidas, o denoising de Variaton total (TV) é outra boa alternativa se o sinal for constante em partes. Pode ser o caso dos dados do acelerômetro, se o seu sinal continuar variando entre diferentes platôs.

Abaixo está um código do Matlab que executa o denoising da TV nesse sinal. O código é baseado no artigo Um método lagrangiano aumentado para restauração de vídeo com variação total . Os parâmetros e devem ser ajustados de acordo com o nível de ruído e as características do sinal.$\mu$$\rho$

Se é o sinal ruidoso e é o sinal a ser avaliado, a função a ser minimizada é , onde é o operador das diferenças finitas.$y$$x$$\mu\|{x-y}\|^2+\|{Dx}\|_1$$D$

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:

O problema é que seu ruído tem um espectro plano. Se você assume um ruído gaussiano branco (o que acaba sendo uma boa suposição), sua densidade de espectro de potência é constante. Grosso modo, isso significa que seu ruído contém todas as frequências. É por isso que qualquer abordagem de frequência, por exemplo, DFT ou filtros passa-baixo, não é boa. Quais seriam suas frequências de corte, já que seu ruído está em todo o espectro?

Uma resposta para essa pergunta é o filtro Wiener, que requer conhecimento das estatísticas do seu ruído e do sinal desejado. Basicamente, o sinal barulhento (sinal + ruído) é atenuado nas frequências em que se espera que o ruído seja maior do que o seu sinal e é amplificado onde seu sinal é maior do que o seu ruído.

No entanto, eu sugeriria abordagens mais modernas que usam processamento não linear, por exemplo, waenoet denoising. Esses métodos fornecem excelentes resultados. Basicamente, o sinal ruidoso é primeiro decomposto em wavelets e depois pequenos coeficientes são zerados. Essa abordagem funciona (e o DFT não) devido à natureza de multi-resolução das wavelets. Ou seja, o sinal é processado separadamente nas bandas de frequência definidas pela transformada wavelet.

No MATLAB, digite 'wavemenu' e depois 'SWT denoising 1-D'. Em seguida, 'Arquivo', 'Análise de exemplo', 'Sinais ruidosos', 'com Haar no nível 5, Blocos ruidosos'. Este exemplo usa a wavelet Haar, que deve funcionar bem para o seu problema.

Eu não sou bom em Python, mas acredito que você pode encontrar alguns pacotes NumPy que executam denoising Haar wavelet.

Conforme sugestão de Daniel Pipa, dei uma olhada no wavelet denoising e encontrei este excelente artigo de Francisco Blanco-Silva.

Aqui, modifiquei seu código Python para que o processamento de imagens funcione com dados 2D (acelerômetro) em vez de dados 3D (imagem).

Observe que o limite é "compensado" pelo limiar suave no exemplo de Francisco. Considere isso e modifique para seu aplicativo.

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)

Onde:

• wavelet- nome da string da forma de wavelet a ser usada (veja pywt.wavelist(), por exemplo 'haar')
• noise_sigma - desvio padrão de ruído dos dados
• data - matriz de valores para filtrar (por exemplo, dados do eixo x, y ou z)