Bagaimana cara menggunakan filter Savitzky Golay untuk menemukan maxima lokal (di antara sampel) dalam sinyal 1D yang diambil secara terpisah?

9

Saya memiliki sinyal seismik y (i): masukkan deskripsi gambar di sini

Di sini saya telah menemukan satu maksimum: i = 152,54, y = 222,29 secara manual dan diplot dengan warna merah.

Saya ingin mencari semua maxima secara otomatis.

Saya membaca bahwa Savitzky Golay Filter (SGF) dapat digunakan untuk menemukan perkiraan halus dari sinyal dan turunannya, dan bahwa salah satu manfaat dari SGF adalah mempertahankan minima dan maxima jauh lebih baik daripada filter lainnya. Ini kedengarannya bagus untuk saya gunakan.

Saya menemukan skrip Matlab yang menghasilkan koefisien SGF. Dan menggunakan ini untuk menemukan bahwa urutan ke 4 koefisien SGF untuk turunannya. Saya membuat kode skrip Matlab kecil itu

  1. menemukan turunan dari sinyal dengan menggabungkan sinyal dengan koefisien SGF urutan ke-4 untuk turunannya
  2. menemukan pasangan sampel (i, i +1) di mana turunan berubah tanda
  3. menemukan zero crossing derivatif dengan interpolasi linier antara i dan i + 1

Naskah:

function [maxX,maxY] = findLocalMax(y)
% Kernel for 4th order Savitzky-Golay filter for finding derivative:
d4 = [0.0724 -0.1195 -0.1625 -0.1061 0 0.1061 0.1625 0.1195 -0.0724];

dy = conv(y,d4,'same'); % derivative

[m n] = size(dy);
maxX = [];
maxY = [];
for i = 1 : n - 1
  if dy(i) < 0 && dy(i+1) > 0 % max somewhere between i and i+1
    a = dy(i)/(dy(i) - dy(i+1)); % linear interpolation
    mx = i + a;
    maxX = [maxX mx];
    my = y(i)*(1-a) + y(i+1)*a; % linear interpolation
    maxY = [maxY my];
  end
end

Saya skrip saya harus menguji apakah turunan berubah dari negatif ke positif untuk mendapatkan fungsi untuk memberikan hasil yang diinginkan, namun ini membingungkan saya. Bukankah seharusnya turunan untuk maksimum berubah dari positif ke negatif? Apakah ada cara yang lebih baik untuk membedakan antara maxima dan minima?

Di bawah ini adalah hasil dari menggunakan fungsi ini untuk menemukan maksimum pada sinyal saya: masukkan deskripsi gambar di sini

Hasilnya terlihat bagus, tetapi saya perhatikan bahwa beberapa maxima tidak ditemukan: i = 143.13, 190.88, 256.97.

Apakah ini karena mereka harus dekat dengan maxima lain?

Bagaimana saya bisa mengendalikan dua maksima terdekat?

Terima kasih sebelumnya atas jawaban apa pun!

Andy
sumber
Bisakah Anda merencanakan output filter?
Jim Clay

Jawaban:

5

Meskipun saya tidak terbiasa dengan jenis filter khusus ini, berdasarkan plot yang telah Anda tunjukkan, saya akan menebak bahwa maxima yang tidak ditemukan oleh proses Anda hanya bertabrakan dengan resolusi waktu yang melekat dalam proses. Segala jenis "penghalusan" menyiratkan bahwa ada beberapa waktu-lokal mengolesi sinyal bunga, sehingga jika ada dua puncak terdekat, ada kemungkinan bahwa mereka akan digabung menjadi satu. Mungkin saja filter tingkat rendah mungkin menyajikan perilaku ini lebih sedikit, kemungkinan dengan mengorbankan jumlah pemulusan yang Anda dapatkan.

Jason R
sumber