Memilih filter yang benar untuk data accelerometer

Saya cukup baru untuk DSP, dan telah melakukan beberapa penelitian tentang kemungkinan filter untuk menghaluskan data accelerometer dengan python. Contoh tipe data yang akan saya alami dapat dilihat pada gambar berikut:

Pada dasarnya, saya mencari saran untuk memuluskan data ini untuk akhirnya mengubahnya menjadi kecepatan dan perpindahan. Saya mengerti bahwa akselerometer dari ponsel sangat bising.

Saya tidak berpikir saya dapat menggunakan filter Kalman saat ini karena saya tidak bisa mendapatkan perangkat untuk referensi kebisingan yang dihasilkan oleh data (saya membaca bahwa penting untuk menempatkan perangkat datar dan menemukan jumlah suara dari pembacaan itu?)

FFT telah menghasilkan beberapa hasil menarik. Salah satu upaya saya adalah untuk FFT sinyal akselerasi, kemudian membuat frekuensi rendah untuk memiliki nilai FFT absolut dari 0. Kemudian saya menggunakan aritmatika omega dan invers FFT untuk mendapatkan plot kecepatan. Hasilnya adalah sebagai berikut:

Apakah ini cara yang baik untuk melakukan sesuatu? Saya mencoba untuk menghilangkan sifat bising keseluruhan sinyal tetapi puncak yang jelas seperti sekitar 80 detik perlu diidentifikasi.

Saya juga lelah menggunakan low pass filter pada data accelerometer asli, yang telah melakukan pekerjaan yang sangat baik untuk melicinkannya, tetapi saya tidak begitu yakin ke mana harus pergi dari sini. Bimbingan tentang ke mana harus pergi dari sini akan sangat membantu!

EDIT: Sedikit kode:

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

Jadi pada dasarnya, saya telah melakukan FFT pada data accelerometer saya, memberikan Sz, menyaring frekuensi tinggi menggunakan filter dinding bata sederhana (saya tahu itu tidak ideal). Kemudian ive gunakan omega aritmatika pada FFT data. Juga terima kasih banyak kepada ahli data untuk menambahkan gambar saya ke dalam posting saya :)

Seperti yang ditunjukkan oleh @JohnRobertson dalam Bag of Tricks for Denoising Signals Sementara Mempertahankan Transisi Tajam , Total Variaton (TV) denoising adalah alternatif lain yang baik jika sinyal Anda konstan konstan. Ini mungkin menjadi kasus untuk data accelerometer, jika sinyal Anda terus bervariasi antara plateaux berbeda.

Di bawah ini adalah kode Matlab yang melakukan denoising TV dalam sinyal seperti itu. Kode ini didasarkan pada kertas Metode Augmented Lagrangian untuk Pemulihan Video Variasi Total . Parameter dan harus disesuaikan sesuai dengan tingkat kebisingan dan karakteristik sinyal.$\mu$$\rho$

Jika adalah sinyal berisik dan adalah sinyal yang diperkirakan, fungsi yang akan diminimalkan adalah , di mana adalah operator perbedaan hingga.$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)];

Hasil:

Masalahnya adalah bahwa suara Anda memiliki spektrum datar. Jika Anda menganggap noise Gaussian putih (yang ternyata menjadi asumsi yang baik) kerapatan spektrum dayanya konstan. Secara kasar, ini berarti bahwa suara Anda mengandung semua frekuensi. Itu sebabnya pendekatan frekuensi apa pun, misalnya DFT atau filter low-pass, tidak bagus. Apa yang akan menjadi frekuensi cut-off Anda karena kebisingan Anda ada di seluruh spektrum?

Satu jawaban untuk pertanyaan ini adalah filter Wiener, yang membutuhkan pengetahuan tentang statistik kebisingan Anda dan sinyal yang Anda inginkan. Pada dasarnya, sinyal berisik (sinyal + noise) dilemahkan pada frekuensi di mana noise diharapkan lebih besar daripada sinyal Anda, dan itu diperkuat di mana sinyal Anda diharapkan lebih besar daripada noise Anda.

Namun, saya akan menyarankan pendekatan yang lebih modern yang menggunakan pemrosesan non-linear, misalnya denoising wavelet. Metode-metode ini memberikan hasil yang sangat baik. Pada dasarnya, sinyal berisik pertama kali didekomposisi menjadi wavelet dan kemudian koefisien kecil menjadi nol. Pendekatan ini berfungsi (dan DFT tidak) karena sifat multi-resolusi wavelet. Artinya, sinyal diproses secara terpisah dalam pita frekuensi yang ditentukan oleh transformasi wavelet.

Dalam MATLAB, ketik 'wavemenu' dan kemudian 'SWT denoising 1-D'. Kemudian 'File', 'Analisis Contoh', 'Sinyal bising', 'dengan Haar di level 5, blok Bising'. Contoh ini menggunakan wavelet Haar, yang seharusnya berfungsi dengan baik untuk masalah Anda.

Saya tidak pandai Python, tapi saya percaya Anda dapat menemukan beberapa paket NumPy yang melakukan denoising Haar wavelet.

Sesuai saran Daniel Pipa, saya melihat denoising wavelet dan menemukan artikel yang bagus ini oleh Francisco Blanco-Silva.

Di sini saya telah memodifikasi kode Python-nya untuk pemrosesan gambar agar bekerja dengan data 2D (akselerometer) daripada 3D (gambar).

Perhatikan , ambang "dibuat" untuk soft-thresholding dalam contoh Francisco. Pertimbangkan ini dan modifikasi untuk aplikasi Anda.

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)

Dimana:

• wavelet- nama string dari bentuk wavelet yang akan digunakan (lihat pywt.wavelist(), misalnya 'haar')
• noise_sigma - standar deviasi kebisingan dari data
• data - array nilai yang akan difilter (mis. data sumbu x, y, atau z)