Dekonwolucja sygnałów 1D zamazana przez jądro Gaussa

Przekształciłem losowy sygnał z aa Gaussa i dodałem szum (w tym przypadku szum Poissona), aby wygenerować głośny sygnał. Teraz chciałbym dekonwolować ten głośny sygnał, aby wyodrębnić oryginalny sygnał przy użyciu tego samego gaussowskiego.

Problem polega na tym, że potrzebuję kodu, który wykonuje dekonwolucję w 1D. (Znalazłem już trochę w 2D, ale moim głównym celem jest 1D).

Czy możesz mi zasugerować kilka pakietów lub programów, które są w stanie to zrobić? (Najlepiej w MATLAB)

Z góry dziękuję za pomoc.

Wyjaśniłem to raz na StackOverflow .

---

Twój sygnał może być reprezentowany jako wektor, a splot jest mnożony przez macierz N-diagonalną (gdzie N jest długością filtra). Zakładam, ze względu na odpowiedź, że filtr jest znacznie mniejszy niż sygnał

Na przykład:

Twój wektor / sygnał to:

V1
V2
...
Vn

Twój filtr (element zwijający) to:

  [b1 b2 b3];

Zatem macierz to nxn: (Niech to się nazywa 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]

Konwolucja to:

A*v;

I dekonwolucja jest

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

Oczywiście w niektórych przypadkach dekonwolucja nie jest możliwa. Są to przypadki, gdy masz liczbę pojedynczą A. Nawet macierze, które nie są liczbą pojedynczą, ale są bliskie byciu liczbą pojedynczą, mogą być problematyczne, ponieważ będą miały duży błąd numeryczny. Możesz to oszacować, obliczając numer warunku macierzy.

Jeśli A ma niski warunek, możesz obliczyć odwrotność i zastosować ją do wyniku.

---

Zobaczmy teraz kilka przykładów w Matlabie:

Po pierwsze, stworzyłem funkcję, która oblicza macierz splotu.

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

Teraz spróbujmy zobaczyć, co dzieje się z różnymi jąderami:

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

Numer warunku to:

 7.8541

Ten jest problematyczny, zgodnie z oczekiwaniami. Po uśrednieniu trudno jest odzyskać oryginalny sygnał.

A teraz spróbujmy łagodniejszego uśrednienia:

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

Wynik to:

1.6667

Jest to zgodne z naszą intuicją, łagodne uśrednianie oryginalnego sygnału jest znacznie łatwiejsze do odwrócenia.

Możemy również zobaczyć, jak wygląda macierz odwrotna:

 figure;imagesc(inv(A));

Oto jedna linia z matrycy:

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

Widzimy, że większość energii w każdej linii koncentruje się w 3-5 współczynnikach wokół centrum. Dlatego, aby dekonwolować, możemy po prostu ponownie zwołać sygnał za pomocą tego przybliżenia:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

To jądro wygląda interesująco! Jest operatorem ostrzenia. Nasza intuicja jest prawidłowa, wyostrzanie eliminuje rozmycie.

Jeśli dodałeś losowy szum, nie możesz uzyskać oryginalnego sygnału ... Możesz spróbować oddzielić sygnały w dziedzinie częstotliwości (jeśli szum i sygnał mają różne częstotliwości). Wygląda jednak na to, że szukasz filtru Wienera .

Myślę, że to wciąż otwarty problem.

Istnieje wiele prac badawczych, które próbują odzyskać oryginalny sygnał najlepiej, jak potrafią.

Klasycznym podejściem są metody oparte na falkach .

Istnieją również podejścia słownikowe takie jak ten .

Bardziej szczegółowe spojrzenie na problem można uzyskać, badając badania przeprowadzone przez Davida L. Donho, Michaela Elada, Alfreda M. Brucksteina itp.

Jeśli poprawnie zrozumiałem problem, możemy go sformalizować w następujący sposób:

Mamy model sygnału,

$y = Hx + \eta$

gdzie $y$ jest obserwacją, $H$ jest operatorem splotu, oraz $\eta$to hałas. Chcemy oszacować$x$ poprzez wykorzystanie obserwacji i znajomości charakterystyk hałasu.

W tym przypadku, $\eta$jest symulowany z rozkładu Poissona. Jednak wyżej wymienione podejścia słownikowe leżą u podstaw założenia szumu Gaussa. W tym przypadku Gaussian jest operatorem splotu, a nie hałasem.

Nie pracowałem nad odzyskiwaniem sygnału pod szumem Poissona, ale poszukałem go i stwierdziłem, że ten artykuł może być przydatny. Podobne podejścia w tym kontekście mogą być przydatne w przypadku tego problemu.

Dekonwolucja hałaśliwych danych jest znanym problemem, ponieważ szum jest arbitralnie powiększany w zrekonstruowanym sygnale. Dlatego do ustabilizowania rozwiązania wymagana jest metoda regularyzacji. Tutaj możesz znaleźć pakiet MATLAB, który rozwiązuje ten problem, implementując algorytm regularyzacji Tichonowa:

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

Przejdę do samego początku pytania. Istnieją funkcje dekonwolucji w MATLAB, które są używane w aplikacjach do przetwarzania obrazu. Tych funkcji można jednak użyć także do sygnałów 1D. Na przykład,

% 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) Dlaczego więc nie rozłożyć sygnału wyjściowego za pomocą hfunkcji i uzyskać (prawie) sygnał wejściowy. Używam tutaj dekonwolucji Wienera

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')

Ewentualnie, jeśli nie znasz hfunkcji, ale znasz wejście i wyjście, tym razem dlaczego nie rozłóż sygnału wejściowego za pomocą wyjścia, które da h^-1funkcję. Następnie możesz użyć go jako filtra do filtrowania szumu. ( 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')

Mam nadzieję, że to pomoże.

Konwolucja to zwielokrotnienie i zsumowanie dwóch sygnałów jeden do drugiego. Mówię o dwóch deterministycznych sygnałach. Jeśli chcesz dekonwolować jeden od drugiego, odpowiada to rozwiązaniu układu równań. Jak być może wiesz, układ równań nie zawsze jest rozwiązywalny. Układ równań może być zawyżony, niedookreślony lub dokładnie rozwiązany.

Jeśli dodasz trochę hałasu, stracisz część informacji i nie możesz ich odzyskać. To, co możesz zrobić, to ponownie rozwiązać liniowy układ równań, biorąc pod uwagę fakt, że do każdego współczynnika dodawany jest składnik szumowy. Lub, jak widać w innej odpowiedzi na twoje pytanie, możesz najpierw oszacować oryginalny sygnał z sygnału hałaśliwego, a następnie spróbować rozwiązać układ równań.

Należy zauważyć, że hałas jest dodawany do pomnożonych i zsumowanych współczynników. Dlatego może się zdarzyć, że twój układ równań ostatecznie nie będzie jednoznacznie rozwiązany. Aby mieć pewność, że jest to unikalne rozwiązanie, macierz współczynników powinna być kwadratowa i mieć pełną rangę.

To byłoby trudne do zrobienia. Konwolucja z gaussowskim jest równoważna pomnożeniu z transformatą Fouriera Gaussa w dziedzinie częstotliwości. Zdarza się, że jest to również gaussowski, w istocie jest to filtr dolnoprzepustowy i naprawdę skuteczny. Po dodaniu hałasu wszystkie informacje znajdujące się w „paśmie zatrzymania” Gaussa zostają zniszczone. Nie ma sposobu, aby to odzyskać.

Dekonwolucja zasadniczo zwielokrotnia się przez odwrotność odpowiedzi częstotliwościowej. Oto problem: odwrotność odpowiedzi częstotliwościowej staje się naprawdę duża, gdy oryginalny Gaussian jest bardzo mały. Przy tych częstotliwościach można zasadniczo zwiększyć hałas w ogromnych ilościach. Nawet jeśli wszystko byłoby całkowicie wolne od hałasu, najprawdopodobniej napotkalibyśmy problemy numeryczne.

Podejścia

Istnieje wiele metod dekonwolucji (mianowicie operator degradacji jest liniowy i niezmienny w czasie / przestrzeni).
Wszyscy próbują poradzić sobie z faktem, że problem jest źle przygotowany w wielu przypadkach.

Lepszymi metodami są te, które dodają pewną regularyzację do modelu danych do przywrócenia.
Mogą to być modele statystyczne (Priors) lub dowolna wiedza.
W przypadku obrazów dobry model jest gładki lub rzadki w gradiencie.

Jednak ze względu na odpowiedź zastosowane zostanie proste podejście parametryczne - Minimalizacja błędu najmniejszych kwadratów między przywróconymi danymi w modelu a pomiarami.

Model

Model najmniejszych kwadratów jest prosty.
Funkcja celu jako funkcja danych jest określona przez:

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

Problem optymalizacji wynika z:

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

Gdzie $x$ to dane do przywrócenia, $h$ to Jądro Rozmycia (w tym przypadku Gaussa) i $y$jest zbiorem podanych pomiarów.
Model zakłada, że ​​pomiary podano tylko dla ważnej części splotu. Mianowicie jeśli$x \in \mathbb{R}^{n}$ i $h \in \mathbb{R}^{k}$ następnie $y \in \mathbb{R}^{m}$ gdzie $m = n - k + 1$.

Jest to operacja liniowa w skończonej przestrzeni, dlatego można ją zapisać za pomocą formularza Matrix:

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

Gdzie $H \in \mathbb{R}^{m \times n}$ jest macierzą splotu.

Rozwiązanie

Rozwiązanie Least Squares daje:

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

Jak widać, wymaga inwersji macierzy.
Możliwość odpowiedniego rozwiązania tego problemu zależy od liczby warunków operatora${H}^{T} H$ który jest posłuszny $\operatorname{cond} \left( H \right) = \sqrt{\operatorname{cond} \left( {H}^{T} H \right)}$.

Analiza numeru stanu

Co kryje się za tym numerem warunku?
Można na nie odpowiedzieć za pomocą algebry liniowej.
Jednak bardziej intuicyjne, moim zdaniem, podejście byłoby myślenie o tym w dziedzinie częstotliwości.

Zasadniczo operator degradacji tłumi energię generalnie o wysokiej częstotliwości.
Ponieważ częstotliwość jest w zasadzie zwielokrotnieniem elementarnym, można powiedzieć, że najłatwiejszym sposobem na odwrócenie tego jest dzielenie elementów przez filtr odwrotny.
Cóż, to co zostało zrobione powyżej.
Problem pojawia się w przypadkach, gdy filtr tłumi energię praktycznie do zera. Mamy wtedy prawdziwe problemy ...
Właśnie to mówi nam numer stanu, jak mocno niektóre częstotliwości były tłumione w stosunku do innych.

Powyżej można zobaczyć numer warunku (za pomocą jednostek [dB]) jako funkcję parametru STD filtra Gaussa.
Zgodnie z oczekiwaniami, im wyższy STD, tym gorsza liczba warunków, ponieważ wyższy STD oznacza silniejszy LPF (wartości malejące na końcu są kwestiami numerycznymi).

Rozwiązanie numeryczne

Utworzono zespół Gaussian Blur Kernel.

Parametry są $n = 300$, $k = 31$ i $m = 270$.
Dane są losowe i nie dodano hałasu.

W MATLAB rozwiązano układ liniowy, w pinv()którym zastosowano Pseudo odwrotność opartą na SVD i \operator.

Jak widać, przy użyciu SVD rozwiązanie jest znacznie mniej czułe, zgodnie z oczekiwaniami.

Dlaczego wystąpił błąd?
Patrząc na rozwiązanie (dla najwyższego STD):

Jak widać, sygnał jest bardzo dobrze przywracany, z wyjątkiem początku i końca.
Wynika to z użycia Valid Convolution, która niewiele mówi nam o tych próbkach.

Hałas

Gdybyśmy dodali hałas, wszystko wyglądałoby inaczej!
Powodem, dla którego wyniki były dobre wcześniej, jest fakt, że MATLAB mógł obsłużyć DR danych i rozwiązać równania, mimo że miały dużą liczbę warunków.

Ale duża liczba warunków oznacza, że ​​filtr odwrotny silnie wzmacnia (Aby odwrócić silne tłumienie) niektóre częstotliwości.
Gdy zawierają szum, oznacza to, że hałas zostanie wzmocniony, a odbudowa będzie zła.

Jak widać powyżej, teraz rekonstrukcja nie będzie działać.

Podsumowanie

Jeśli ktoś dokładnie zna Operatora Degradacji, a SNR jest bardzo dobry, zadziałają proste metody dekonwolucji.
Głównym problemem dekonwolucji jest to, jak mocno Operator Degradacji tłumi częstotliwości.
Im bardziej tłumi, tym więcej SNR jest potrzebnych do przywrócenia (jest to w zasadzie idea filtra Wienera ).
Częstotliwości ustawione na zero nie mogą zostać przywrócone!

W praktyce, aby uzyskać stabilne wyniki, należy dodać kilka priorytetów.

Kod jest dostępny w moim repozytorium GitHub Stack29change Signal Processing Q2969 .

Ogólnie rzecz biorąc, jednym ze sposobów rozwiązania problemu, który zasadniczo uogólnia problem wyodrębnienia dwóch lub więcej składników, jest pobranie widm G¹, G² ⋯, G # sygnałów # 1, # 2, ..., #n, zestawienie sumaryczne kwadrat Γ (ν) = | G¹ (ν) | ² + | G² (ν) | ² + ⋯ + | Gⁿ (ν) | ² przy każdej częstotliwości ν i normalizuj G₁ (ν) ≡ G¹ (ν) * / Γ (ν), G₂ (ν) ≡ G² (ν) * / Γ (ν), ..., G_n (ν) ≡ Gⁿ (ν) * / Γ (ν). Problem niedokładności i szumu odpowiada temu, że Γ (ν) ~ 0 jest możliwe dla niektórych częstotliwości ν. Aby sobie z tym poradzić, dodaj kolejny „sygnał”, aby wyodrębnić G⁰ (ν) = stały - sygnał „szumu”. Teraz Γ (ν) będzie ściśle ograniczone poniżej. Jest to prawie na pewno związane z regularyzacją Tichonowa, ale nigdy nie znalazłem ani nie ustaliłem żadnego wyniku równoważności ani innej korespondencji. Jest prostszy, bardziej bezpośredni i intuicyjny.

Alternatywnie można traktować G jako wektory wyposażone w odpowiedni iloczyn wewnętrzny, np. «G, G '» ≡ ∫ G (ν) * G' (ν) dν, i przyjmować (G₀, G₁, ⋯, G_n) jako podwójny (np. uogólniona odwrotność) (G⁰, G¹, ⋯, Gⁿ) - zakładając oczywiście, że wektory składowe są liniowo niezależne.

Dla dekonwolucji Gaussa można ustawić n = 1, G⁰ = sygnał „szumu”, a G¹ = sygnał „Gaussa”.

Odpowiedź udzielona przez Andreya Rubshteina nie powiedzie się dobrze w obecności hałasu, ponieważ opisany problem jest bardzo wrażliwy na hałas i błędy modelowania. Dobrym pomysłem jest skonstruowanie matrycy splotowej, ale użycie regularyzacji w inwersji jest absolutną koniecznością w przypadku takiego problemu. Bardzo prostą i bezpośrednią metodą regularyzacji (choć kosztowną pod względem obliczeniowym) jest skrócona dekompozycja wartości pojedynczej (TSVD). Metody takie jak regularyzacja Tichonowa i regularyzacja całkowitej zmiennościwarto sprawdzić. Regulararyzacja Tichonowa (i jej ogólna postać) ma bardzo elegancką formę stosową, którą łatwo wdrożyć w Matlabie. Sprawdź książkę: Liniowe i nieliniowe problemy odwrotne z praktycznymi zastosowaniami Samuli Siltanen i Jennifer Mueller.

W rzeczywistości pytanie nie jest jasne. Ale odpowiedzi potwierdziły to, o co prosiłeś. Możesz zbudować układ liniowych równań algebraicznych, jak sugerują niektórzy ludzie, to prawda, ale macierz zbudowana na znanym sygnale jest tak zwana słabo uwarunkowana. Oznacza to, że gdy spróbujesz go odwrócić, błędy skracania zabijają rozwiązanie i otrzymujesz losowe liczby. Powszechnym podejściem jest ograniczenie ekstremum. Minimalizujesz normę rozwiązania || x || z ograniczeniem || Ax - y || <delta. Więc szukasz takiego x z najmniejszą normą, która nie pozwoli, aby różnica między Ax i Y była duża. Bardzo proste jest dodanie tak zwanego parametru regularyzacji na głównej przekątnej macierzy uzyskanej przy zastosowaniu najmniejszych kwadratów. Nazywa się to regularyzacją Tichonowa. Mam kodujące próbki, które to robią,