Menentukan simpangan rata-rata dan standar dalam waktu nyata

Apa yang akan menjadi cara ideal untuk menemukan mean dan standar deviasi sinyal untuk aplikasi waktu nyata. Saya ingin dapat memicu controller ketika sinyal lebih dari 3 standar deviasi dari rata-rata untuk jangka waktu tertentu.

Saya mengasumsikan DSP khusus akan melakukan ini dengan mudah, tetapi apakah ada "jalan pintas" yang mungkin tidak memerlukan sesuatu yang begitu rumit?

Ada kekurangan dalam jawaban Jason R, yang dibahas dalam Knuth's "Art of Computer Programming" vol. 2. Masalah muncul jika Anda memiliki standar deviasi yang merupakan sebagian kecil dari rata-rata: perhitungan E (x ^ 2) - (E (x) ^ 2) menderita sensitivitas parah hingga kesalahan pembulatan titik mengambang.

Anda bahkan dapat mencobanya sendiri dalam skrip Python:

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]] 
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

Saya mendapatkan -128.0 sebagai jawaban, yang jelas tidak valid secara komputasi, karena matematika memperkirakan bahwa hasilnya harus non-negatif.

Knuth mengutip pendekatan (saya tidak ingat nama penemunya) untuk menghitung rata-rata berjalan dan standar deviasi yang kira-kira seperti ini:

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

dan kemudian setelah setiap langkah, nilai dari mmean, dan standar deviasi dapat dihitung sebagai sqrt(S/n)atau sqrt(S/n-1)tergantung pada definisi favorit Anda tentang standar deviasi.

Persamaan yang saya tulis di atas sedikit berbeda dari yang ada di Knuth, tetapi ini setara secara komputasi.

Ketika saya memiliki beberapa menit lagi, saya akan memberi kode rumus di atas dalam Python dan menunjukkan bahwa Anda akan mendapatkan jawaban yang tidak negatif (yang diharapkan dekat dengan nilai yang benar).

pembaruan: ini dia.

test1.py:

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5] 
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

hasil:

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

Anda akan mencatat bahwa masih ada beberapa kesalahan pembulatan, tetapi itu tidak buruk, sedangkan naive_statshanya muntah.

sunting: Hanya memperhatikan komentar Belisarius yang mengutip Wikipedia yang menyebutkan algoritma Knuth.

Apa yang akan menjadi cara ideal untuk menemukan mean dan standar deviasi sinyal untuk aplikasi waktu nyata. Saya ingin dapat memicu controller ketika sinyal lebih dari 3 standar deviasi dari rata-rata untuk jangka waktu tertentu.

Pendekatan yang tepat dalam situasi seperti ini biasanya untuk menghitung rata-rata berlari tertimbang secara eksponensial dan standar deviasi. Dalam eksponensial tertimbang rata-rata, perkiraan mean dan varians bias terhadap sampel terbaru memberikan Anda memperkirakan mean dan varians selama terakhir detik$\tau$ , yang mungkin apa yang Anda inginkan, daripada biasa aritmatika rata-rata semua sampel yang pernah dilihat.

Dalam domain frekuensi, "rata-rata berlari tertimbang secara eksponensial" hanyalah kutub nyata. Mudah diterapkan di domain waktu.

Implementasi domain waktu

Biarkan meandan meansqjadilah estimasi saat ini dari mean dan mean dari kuadrat dari sinyal. Pada setiap siklus, perbarui estimasi ini dengan sampel baru x:

% update the estimate of the mean and the mean square:
mean = (1-a)*mean + a*x
meansq = (1-a)*meansq + a*(x^2)

% calculate the estimate of the variance:
var = meansq - mean^2;

% and, if you want standard deviation:
std = sqrt(var);

Di sini adalah konstanta yang menentukan panjang efektif rata-rata berjalan. Bagaimana memilih dijelaskan di bawah ini di "analisis".$0 < a < 1$$a$

Apa yang dinyatakan di atas sebagai program imperatif juga dapat digambarkan sebagai diagram aliran sinyal:

Analisis

Algoritma di atas menghitung mana adalah input pada sampel , dan adalah output (yaitu estimasi rata-rata). Ini adalah filter IIR tiang tunggal yang sederhana. Mengambil transformasi , kami menemukan fungsi transfer .$y_i = a x_i + (1-a) y_{i-1}$$x_i$$i$$y_i$$z$

Mengondensasi filter IIR ke dalam blok mereka sendiri, diagram sekarang terlihat seperti ini:

Untuk menuju ke domain kontinu, kami membuat substitusi mana adalah waktu sampel dan adalah laju sampel. Memecahkan , kami menemukan bahwa sistem kontinu memiliki kutub pada . $z = e^{s T}$$T$$f_s = 1/T$$1-(1-a)e^{-sT}=0$$s = \frac{1}{T} \log (1-a)$

Pilih :$a$

Referensi

• Sumber diagram Simulink dapat diunduh dari https://gist.github.com/1942771

Metode yang saya gunakan sebelumnya dalam aplikasi pemrosesan tertanam adalah untuk mempertahankan akumulator jumlah dan jumlah kuadrat dari sinyal yang menarik:

Juga, pantau instan waktu saat ini dalam persamaan di atas (yaitu, perhatikan jumlah sampel yang telah Anda tambahkan ke akumulator). Kemudian, mean sampel dan simpangan baku pada saat adalah:$i$$i$

atau Anda dapat menggunakan:

tergantung pada metode estimasi standar deviasi yang Anda inginkan . Persamaan ini didasarkan pada definisi varians :

Saya telah menggunakan ini dengan sukses di masa lalu (meskipun saya hanya peduli dengan estimasi varians, bukan standar deviasi), meskipun Anda harus berhati-hati tentang tipe numerik yang Anda gunakan untuk memegang akumulator jika Anda akan menjumlahkan periode waktu yang lama; Anda tidak ingin meluap.

Sunting: Selain komentar di atas tentang overflow, harus dicatat bahwa ini bukan algoritma yang kuat secara numerik ketika diterapkan dalam aritmatika titik-mengambang, berpotensi menyebabkan kesalahan besar dalam perkiraan statistik. Lihatlah jawaban Jason S untuk pendekatan yang lebih baik dalam kasus itu.

Mirip dengan jawaban yang disukai di atas (Jason S.), dan juga berasal dari rumus yang diambil dari Knut (Vol.2, p 232), orang juga dapat memperoleh rumus untuk mengganti nilai, yaitu menghapus dan menambahkan nilai dalam satu langkah . Menurut pengujian saya, ganti memberikan presisi yang lebih baik daripada versi dua-langkah hapus / tambah.

Kode di bawah ini di Jawa, meandan sdapatkan pembaruan (variabel anggota "global"), sama seperti mdan sdi atas dalam posting Jason. Nilai countmengacu pada ukuran jendela n.

/**
 * Replaces the value {@code x} currently present in this sample with the
 * new value {@code y}. In a sliding window, {@code x} is the value that
 * drops out and {@code y} is the new value entering the window. The sample
 * count remains constant with this operation.
 * 
 * @param x
 *            the value to remove
 * @param y
 *            the value to add
 */
public void replace(double x, double y) {
    final double deltaYX = y - x;
    final double deltaX = x - mean;
    final double deltaY = y - mean;
    mean = mean + deltaYX / count;
    final double deltaYp = y - mean;
    final double countMinus1 = count - 1;
    s = s - count / countMinus1 * (deltaX * deltaX - deltaY * deltaYp) - deltaYX * deltaYp / countMinus1;
}

Jawaban Jason dan Nibot berbeda dalam satu aspek penting: Metode Jason menghitung std dev dan rata-rata untuk seluruh sinyal (karena y = 0), sedangkan Nibot adalah perhitungan "berjalan", yaitu menimbang sampel yang lebih baru lebih kuat daripada sampel dari masa lalu yang jauh.

Karena aplikasi tampaknya memerlukan std dev dan berarti sebagai fungsi waktu, metode Nibot mungkin yang lebih tepat (untuk aplikasi spesifik ini). Namun, bagian yang sulit sebenarnya adalah untuk mendapatkan bagian pembobotan waktu yang tepat. Contoh Nibot menggunakan filter tiang tunggal sederhana.

$E[x]$$x[n]$$E[x^2]$$x[n]^2$

Pilihan filter lowpass dapat dipandu oleh apa yang Anda ketahui tentang sinyal Anda dan resolusi waktu yang Anda butuhkan untuk estimasi Anda. Frekuensi cutoff yang lebih rendah dan urutan yang lebih tinggi akan menghasilkan akurasi yang lebih baik tetapi waktu respons lebih lambat.

Untuk menyulitkan lebih lanjut, satu filter diterapkan dalam domain linier dan lainnya dalam domain kuadrat. Mengkuadratkan secara signifikan mengubah konten spektral sinyal sehingga Anda mungkin ingin menggunakan filter yang berbeda dalam domain kuadrat.

Berikut adalah contoh tentang cara memperkirakan mean, rms dan std dev sebagai fungsi waktu.

%% example
fs = 44100; n = fs; % 44.1 kHz sample rate, 1 second
% signal: white noise plus a low frequency drift at 5 Hz)
x = randn(n,1) + sin(2*pi*(0:n-1)'*5/fs);
% mean estimation filter: since we are looking for effects in the 5 Hz range we use maybe a
% 25 Hz filter, 2nd order so it's not too sluggish
[b,a] = butter(2,25*2/fs);
xmeanEst = filter(b,a,x);
% now we estimate x^2, since most frequency double we use twice the bandwidth
[b,a] = butter(2,50*2/fs);
x2Est = filter(b,a,x.^2);
% std deviation estimate
xstd = sqrt(x2Est)-xmeanEst;
% and plot it
h = plot([x, xmeanEst sqrt(x2Est) xstd]);
grid on;
legend('x','E<x>','sqrt(E<x^2>)','Std dev');
set(h(2:4),'Linewidth',2);