Dekonvolusi Sinyal 1D Dikaburkan oleh Kernel Gaussian

Saya telah menggabungkan sinyal acak dengan aa Gaussian dan menambahkan noise (Poisson noise dalam kasus ini) untuk menghasilkan sinyal berisik. Sekarang saya ingin mendekonvolusi sinyal berisik ini untuk mengekstrak sinyal asli menggunakan Gaussian yang sama.

Masalahnya adalah bahwa saya memerlukan kode yang melakukan pekerjaan dekonvolusi dalam 1D. (Saya sudah menemukan beberapa di 2D tapi tujuan utama saya adalah 1D).

Bisakah Anda menyarankan saya beberapa paket atau program yang dapat melakukannya? (Lebih disukai dalam MATLAB)

Terima kasih sebelumnya atas bantuannya.

Saya sudah menjelaskannya sekali di StackOverflow .

---

Sinyal Anda dapat direpresentasikan sebagai vektor, dan konvolusi adalah perkalian dengan matriks N-diagonal (di mana N adalah panjang filter). Saya mengasumsikan demi jawaban bahwa filter jauh lebih kecil daripada sinyal

Sebagai contoh:

Vektor / sinyal Anda adalah:

V1
V2
...
Vn

Filter Anda (elemen convolving) adalah:

  [b1 b2 b3];

Jadi matriksnya adalah nxn: (Biarkan disebut A):

[b2 b3 0  0  0  0.... 0]
[b1 b2 b3 0  0  0.... 0]
[0  b1 b2 b3 0  0.... 0]
.....
[0  0  0  0  0  0...b2 b3]

Konvolusi adalah:

A*v;

Dan de-convolution adalah

A^(-1) * ( A) * v;

Jelas, dalam beberapa kasus de-konvolusi tidak dimungkinkan. Ini adalah kasus-kasus ketika Anda memiliki A. tunggal Bahkan matriks yang tidak tunggal, tetapi hampir tunggal, dapat menjadi masalah, karena mereka akan memiliki kesalahan numerik yang besar. Anda dapat memperkirakannya dengan menghitung nomor kondisi dari matriks.

Jika A memiliki kondisi rendah, Anda dapat menghitung kebalikannya, dan menerapkannya pada hasilnya.

---

Sekarang, mari kita lihat beberapa contoh di Matlab:

Pertama, saya telah membuat fungsi yang menghitung matriks konvolusi.

function A = GetConvolutionMatrix(b,numA)
    A = zeros(numA,numA);
    vec = [b  zeros(1,numA-numel(b))];
    for i=1:size(A,1)
        A(i,:) = circshift(vec,[1 i]);
    end
end

Sekarang, mari kita coba melihat apa yang terjadi dengan kernel yang berbeda:

    b = [1 1 1];
    A = GetConvolutionMatrix(b,10);
    disp(cond(A));

Nomor kondisinya adalah:

 7.8541

Yang ini bermasalah, seperti yang diharapkan. Setelah rata-rata, sulit untuk mendapatkan kembali sinyal aslinya.

Sekarang, mari kita coba rata-rata lebih ringan:

b = [0.1 0.8 0.1];
A = GetConvolutionMatrix(b,10);
disp(cond(A));

Hasilnya adalah:

1.6667

Itu sesuai dengan intuisi kita, rata-rata sinyal asli yang ringan jauh lebih mudah untuk dibalik.

Kita juga dapat melihat bagaimana matriks invers terlihat:

 figure;imagesc(inv(A));

Berikut ini satu baris dari matriks:

  0.0003   -0.0026    0.0208   -0.1640    1.2910   -0.1640    0.0208   -0.0026    0.0003   -0.0001

Kita dapat melihat bahwa sebagian besar energi di setiap garis terkonsentrasi dalam 3-5 koefisien di sekitar pusat. Oleh karena itu, untuk melakukan dekonvolusi, kita dapat dengan mudah menggabungkan sinyal lagi dengan perkiraan ini:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

Kernel ini terlihat menarik! Ini adalah operator penajaman. Intuisi kita benar, mengasah membatalkan kabur.

Jika Anda telah menambahkan noise acak, Anda tidak dapat memperoleh sinyal asli ... Anda dapat mencoba memisahkan sinyal dalam domain frekuensi (jika noise dan sinyal dari frekuensi yang berbeda). Tetapi tampaknya apa yang Anda cari adalah filter Wiener .

Saya pikir ini masih merupakan masalah terbuka.

Ada banyak makalah penelitian yang mencoba memulihkan sinyal asli sebaik mungkin.

Salah satu pendekatan klasik adalah melalui Metode berbasis Wavelet .

Ada juga pendekatan kamus seperti ini .

Anda bisa mendapatkan pandangan yang lebih mendalam tentang masalah dengan mengikuti penelitian yang dilakukan oleh David L. Donho, Michael Elad, Alfred M. Bruckstein dll.

Jika saya memahami masalahnya dengan benar, kami dapat memformalkan masalahnya sebagai berikut:

Kami memiliki model sinyal,

$y = Hx + \eta$

$y$$H$$\eta$$x$

$\eta$

Saya tidak bekerja pada pemulihan sinyal di bawah kebisingan Poisson, tetapi saya googled dan menemukan makalah ini mungkin berguna. Pendekatan serupa dalam konteks itu dapat bermanfaat untuk masalah ini.

Dekonvolusi data yang berisik dikenal sebagai masalah yang buruk, karena kebisingan diperbesar secara sewenang-wenang dalam sinyal yang direkonstruksi. Oleh karena itu, diperlukan metode regularisasi untuk menstabilkan solusi. Di sini, Anda dapat menemukan paket MATLAB yang membahas masalah ini dengan menerapkan algoritma regularisasi Tikhonov:

https://github.com/soheil-soltani/TranKin .

Saya akan pergi ke awal pertanyaan. Ada fungsi dekonvolusi dalam MATLAB yang digunakan untuk aplikasi pemrosesan gambar. Namun, Anda juga dapat menggunakan fungsi ini untuk sinyal 1D. Sebagai contoh,

% a random signal
sig_clean = zeros(1,200); 
sig_clean(80:100)=100;

figure
subplot(1,3,1)
plot(sig_clean,'b-.','LineWidth',2)
legend('Clean Signal')

% convolve it with a gaussian
x=1:30;
h = exp(-(x-15).^2/20); h=h/sum(h);
sig_noisy = conv(sig_clean,h,'same');

% and add noise
sig_noisy = awgn(sig_noisy,0,'measured');

subplot(1,3,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',3)
legend('Blurred and noise added signal','Clean Signal')

( sig_noisy = sig_clean * h + noise) Lalu mengapa tidak mendekonvolusikan sinyal output dengan hfungsi dan mendapatkan (hampir) sinyal input. Saya menggunakan dekonvolusi Wiener di sini

sig_deconvolved=deconvwnr(sig_noisy,h,1);

subplot(1,3,3)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_deconvolved,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Deconvolved Signal')

Atau, jika Anda tidak tahu hfungsi, tetapi tahu input dan output, kali ini mengapa tidak mendekonvolusikan sinyal input dengan output yang akan memberikan h^-1fungsi. Kemudian Anda dapat menggunakannya sebagai filter untuk menyaring sinyal berisik. ( sig_clean = sig_noisy * h^-1)

h_inv=deconvwnr(sig_clean,sig_noisy,1);

figure;
subplot(1,2,1)
plot(h_inv)
legend('h^-^1')


sig_filtered=conv(sig_noisy,h_inv,'same');
subplot(1,2,2)
plot(sig_noisy,'r')
hold on, plot(sig_clean,'b-.','LineWidth',2)
hold on, plot(sig_filtered,'k--','LineWidth',2)
legend('Blurred and noise added signal','Clean Signal','Filtered Signal')

Saya harap ini membantu.

Konvolusi adalah perkalian dan penjumlahan dari dua sinyal satu ke yang lain. Saya berbicara tentang dua sinyal deterministik. Jika Anda ingin memisahkan satu dari yang lain maka ini sesuai dengan solusi sistem persamaan. Seperti yang Anda ketahui, sistem persamaan tidak selalu dapat dipecahkan. Sistem persamaan dapat ditentukan secara berlebihan, tidak ditentukan atau dapat dipecahkan dengan tepat.

Jika Anda menambahkan beberapa kebisingan, maka Anda kehilangan beberapa informasi dan Anda tidak bisa mendapatkan informasi ini kembali. Apa yang dapat Anda lakukan adalah menyelesaikan lagi sistem persamaan linear dengan mempertimbangkan fakta bahwa setiap koefisien ditambahkan dengan istilah noise. Atau seperti yang Anda lihat di jawaban lain untuk pertanyaan Anda, Anda mungkin ingin terlebih dahulu memperkirakan sinyal asli dari sinyal berisik dan kemudian mencoba menyelesaikan sistem persamaan.

Penting untuk dicatat bahwa noise ditambahkan ke koefisien yang dikalikan dan disimpulkan. Oleh karena itu, bahkan mungkin kasus bahwa sistem persamaan Anda pada akhirnya tidak dapat dipecahkan secara unik. Untuk memastikan bahwa ini dapat dipecahkan secara unik, matriks koefisien Anda harus kuadrat dan peringkat penuh.

Ini akan sulit dilakukan. Konvolusi dengan Gaussian setara dengan perkalian dengan Transformasi Fourier Gaussian dalam domain frekuensi. Ini kebetulan juga merupakan Gaussian pada dasarnya ini adalah low pass filter dan yang benar-benar efektif. Setelah Anda menambahkan noise, semua informasi yang ada di "stop band" Gaussian dihancurkan. Tidak ada cara untuk memulihkannya.

De-convolution pada dasarnya mengalikan dengan kebalikan dari respon frekuensi. Inilah masalahnya: Kebalikan dari respons frekuensi menjadi sangat, sangat besar di mana Gaussian asli sangat kecil. Pada frekuensi ini pada dasarnya Anda akan memperkuat kebisingan dengan jumlah besar. Bahkan jika semuanya benar-benar bebas noise, kemungkinan besar Anda akan mengalami masalah numerik.

Pendekatan

Ada banyak metode untuk Dekonvolusi (Yaitu operator degradasi adalah linear dan Time / Space Invariant) di luar sana.
Semua dari mereka mencoba untuk berurusan dengan fakta bahwa masalahnya adalah Poised Poised dalam banyak kasus.

Metode yang lebih baik adalah metode yang menambahkan beberapa regularisasi pada model data yang akan dipulihkan.
Ini bisa berupa model statistik (Priors) atau pengetahuan apa pun.
Untuk gambar, model yang bagus adalah potongan yang bijak, halus atau jarangnya gradien.

Tetapi demi jawabannya, pendekatan parametrik sederhana akan dilakukan - -Minimalkan Least Squares Error antara data yang dipulihkan dalam model dengan pengukuran.

Model

Model kuadrat terkecil sederhana.
Fungsi objektif sebagai fungsi data diberikan oleh:

$$f \left( x \right) = \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$$

Masalah optimasi diberikan oleh:

$$\arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| h \ast x - y \right\|_{2}^{2}$$

$x$$h$$y$
$x \in \mathbb{R}^{n}$$h \in \mathbb{R}^{k}$$y \in \mathbb{R}^{m}$$m = n - k + 1$

Ini adalah operasi linear dalam ruang terbatas sehingga dapat ditulis menggunakan Formulir Matriks:

$$\arg \min_{x} f \left( x \right) = \arg \min_{x} \frac{1}{2} \left\| H x - y \right\|_{2}^{2}$$

$H \in \mathbb{R}^{m \times n}$

Larutan

Solusi Least Squares diberikan oleh:

$$\hat{x} = \left( {H}^{T} H \right)^{-1} {H}^{T} y$$

${H}^{T} H$$\operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)}$

Analisis Nomor Kondisi

Apa yang ada di balik nomor kondisi ini?
Orang bisa menjawabnya menggunakan Aljabar Linier.
Namun, menurut saya, pendekatan yang lebih intuitif akan memikirkannya dalam Domain Frekuensi.

Pada dasarnya operator degradasi mengurangi energi, pada umumnya, frekuensi tinggi.
Sekarang, karena pada frekuensi ini pada dasarnya adalah perkalian elemen yang bijak, orang akan mengatakan cara mudah untuk membalikkannya adalah pembagian elemen bijak oleh filter terbalik.
Ya, itulah yang dilakukan di atas.
Masalah muncul dengan kasus filter melemahkan energi praktis menjadi nol. Maka kita memiliki masalah nyata ...
Ini pada dasarnya adalah apa yang dikatakan Nomor Kondisi, betapa sulitnya beberapa frekuensi dilemahkan relatif terhadap yang lain.

Di atas orang dapat melihat Nomor Kondisi (Menggunakan unit [dB]) sebagai fungsi dari parameter STD Filter Gaussian.
Seperti yang diharapkan, semakin tinggi STD, semakin buruk angka kondisinya karena STD yang lebih tinggi berarti LPF yang lebih kuat (Nilai yang turun pada akhirnya adalah masalah numerik).

Solusi Numerik

Ensemble Gaussian Blur Kernel telah dibuat.

$n = 300$$k = 31$$m = 270$

Dalam MATLAB, Sistem Linier diselesaikan dengan pinv()menggunakan Pseudo Inverse dan \operator berbasis SVD .

Seperti yang dapat kita lihat, menggunakan SVD solusinya jauh kurang sensitif seperti yang diharapkan.

Mengapa ada kesalahan?
Mencari solusi (Untuk STD tertinggi):

Karena orang dapat melihat sinyal dipulihkan dengan sangat baik kecuali untuk awal dan akhir.
Hal ini disebabkan oleh penggunaan Konvolusi Valid yang memberi tahu kita sedikit tentang sampel tersebut.

Kebisingan

Jika kami menambahkan noise, semuanya akan terlihat berbeda!
Alasan hasil yang baik sebelumnya adalah karena fakta bahwa MATLAB dapat menangani DR dari data dan menyelesaikan persamaan meskipun mereka memiliki nomor kondisi yang besar.

Tetapi angka kondisi besar berarti filter terbalik menguatkan kuat (Untuk membalikkan redaman kuat) beberapa frekuensi.
Ketika itu mengandung noise itu berarti noise akan diperkuat dan restorasi akan menjadi buruk.

Seperti yang bisa dilihat di atas, sekarang rekonstruksi tidak akan berhasil.

Ringkasan

Jika seseorang mengetahui Operator Degradasi dengan tepat dan SNR sangat bagus, metode dekonvolusi sederhana akan berfungsi.
Masalah utama dari dekonvolusi adalah seberapa keras Operator Degradasi melemahkan frekuensi.
Semakin ia mengurangi SNR yang dibutuhkan untuk memulihkan (Ini pada dasarnya adalah ide di balik Wiener Filter ).
Frekuensi yang ditetapkan ke nol tidak dapat dipulihkan!

Dalam praktiknya, agar hasil yang stabil harus ditambahkan beberapa prior.

Kode ini tersedia di Pemrosesan Sinyal StackExchange saya, Q2969 GitHub Repository .

Secara umum, satu metode untuk menangani masalah yang menggeneralisasi secara substansial ke masalah mengekstraksi dua atau lebih komponen adalah dengan mengambil spektra G¹, G² ⋯, Gⁿ sinyal # 1, # 2, ..., #n, mentabulasi total kuadrat Γ (ν) = | G¹ (ν) | ² + | G² (ν) | ² + ⋯ + | Gⁿ (ν) | ² pada setiap frekuensi ν, dan menormalkan G₁ (ν) ≡ G¹ (ν) * / Γ (ν), G₂ (ν) ≡ G² (ν) * / Γ (ν), ..., G_n (ν) ≡ Gⁿ (ν) * / Γ (ν). Masalah dengan ketelitian dan kebisingan tidak sesuai dengan kenyataan bahwa Γ (ν) ~ 0 dimungkinkan untuk beberapa frekuensi ν. Untuk menangani ini, tambahkan "sinyal" lain untuk mengekstrak G⁰ (ν) = konstan - sinyal "noise". Sekarang Γ (ν) akan dibatasi secara ketat di bawah ini. Ini hampir pasti terkait dengan regularisasi Tikhonov, tetapi saya tidak pernah menemukan atau menetapkan hasil kesetaraan atau korespondensi lainnya. Ini lebih sederhana dan lebih langsung dan intuitif.

Sebagai alternatif, Anda dapat memperlakukan G's sebagai vektor yang dilengkapi dengan produk dalam yang sesuai, misalnya «G, G '» ≡ ∫ G (ν) * G' (ν) dν, dan menganggap (G₀, G₁, ⋯, G_n) sebagai dual (mis. invers umum) dari (G⁰, G¹, ⋯, Gⁿ) - dengan asumsi, tentu saja, vektor komponen bebas linear.

Untuk dekonvolusi Gaussian, seseorang akan mengatur n = 1, G⁰ = sinyal "noise" dan G¹ = sinyal "Gaussian".

Jawaban yang diberikan oleh Andrey Rubshtein akan gagal total di hadapan kebisingan, karena masalah yang dijelaskan sangat sensitif terhadap kebisingan dan kesalahan pemodelan. Merupakan ide yang bagus untuk membangun matriks konvolusi, tetapi penggunaan regularisasi dalam inversi adalah mutlak harus dalam masalah seperti ini. Metode regularisasi yang sangat sederhana dan lurus ke depan (walaupun mahal secara komputasional) adalah Truncated Singular Value Decomposition (TSVD). Metode seperti Regulatorisasi Tikhonov dan Regulasi Variasi Totallayak dicoba. Regularisasi Tikhonov (dan bentuk umumnya) memiliki bentuk ditumpuk sangat elegan yang mudah diimplementasikan di Matlab. Lihat buku: Masalah invers linier dan nonlinier dengan aplikasi praktis oleh Samuli Siltanen dan Jennifer Mueller.

Sebenarnya pertanyaannya tidak jelas. Tetapi jawabannya menjelaskan apa yang Anda minta. Anda dapat membangun sistem persamaan aljabar linier seperti saran beberapa orang, itu benar, tetapi matriks yang dibangun berdasarkan sinyal yang dikenal disebut dengan kondisi buruk. Itu berarti ketika Anda mencoba untuk membalikkannya, kesalahan pemotongan membunuh solusi dan Anda menerima angka acak dalam hasilnya. Pendekatan umum dibatasi ekstrem. Anda meminimalkan norma solusi || x || dengan constrain || Axe - y || <delta. Jadi Anda mencari x seperti itu dengan norma terkecil yang tidak memungkinkan perbedaan antara Ax dan y menjadi besar. Sangat sederhana Anda perlu menambahkan apa yang disebut parameter regularisasi pada diagonal utama matriks yang diperoleh dalam penerapan kuadrat terkecil. Ini disebut regularisasi Tikhonov. Saya memiliki sampel pengkodean yang melakukan itu,