Entfaltung von 1D-Signalen, die durch einen Gaußschen Kern verwischt wurden

Ich habe ein zufälliges Signal mit einem Gaußschen Signal gefaltet und Rauschen (in diesem Fall Poisson-Rauschen) hinzugefügt, um ein verrauschtes Signal zu erzeugen. Jetzt möchte ich dieses verrauschte Signal auflösen, um das ursprüngliche Signal mit demselben Gaußschen zu extrahieren.

Das Problem ist, dass ich einen Code brauche, der die Entfaltung in 1D erledigt. (Ich habe bereits einige in 2D gefunden, aber mein Hauptziel ist 1D).

Können Sie mir bitte einige Pakete oder Programme vorschlagen, die dazu in der Lage sind? (Vorzugsweise in MATLAB)

Vielen Dank im Voraus für die Hilfe.

Ich habe es einmal auf StackOverflow erklärt .

---

Ihr Signal kann als Vektor dargestellt werden, und die Faltung ist eine Multiplikation mit einer N-Diagonalmatrix (wobei N die Länge des Filters ist). Ich gehe aus Gründen der Antwort davon aus, dass der Filter viel kleiner als das Signal ist

Zum Beispiel:

Ihr Vektor / Signal ist:

V1
V2
...
Vn

Ihr Filter (Faltungselement) lautet:

  [b1 b2 b3];

Die Matrix lautet also nxn: (Sei 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]

Faltung ist:

A*v;

Und Entfaltung ist

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

Offensichtlich ist in einigen Fällen eine Entfaltung nicht möglich. Dies sind die Fälle, in denen Sie Singular A haben. Selbst Matrizen, die nicht Singular sind, aber fast Singular sind, können problematisch sein, da sie einen großen numerischen Fehler aufweisen. Sie können es schätzen, indem Sie die Bedingungsnummer der Matrix berechnen.

Wenn A eine niedrige Bedingung hat, können Sie die Umkehrung berechnen und auf das Ergebnis anwenden.

---

Schauen wir uns nun einige Beispiele in Matlab an:

Zuerst habe ich eine Funktion erstellt, die die Faltungsmatrix berechnet.

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

Versuchen wir nun zu sehen, was mit verschiedenen Kerneln passiert:

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

Die Bedingungsnummer lautet:

 7.8541

Dieser ist erwartungsgemäß problematisch. Nach der Mittelwertbildung ist es schwierig, das ursprüngliche Signal zurückzugewinnen.

Versuchen wir nun eine mildere Mittelung:

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

Das Ergebnis ist:

1.6667

Das passt gut zu unserer Intuition. Eine milde Mittelung des Originalsignals ist viel einfacher umzukehren.

Wir können auch sehen, wie die inverse Matrix aussieht:

 figure;imagesc(inv(A));

Hier ist eine Zeile aus der Matrix:

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

Wir können sehen, dass der größte Teil der Energie in jeder Linie in 3-5 Koeffizienten um das Zentrum konzentriert ist. Um sich zu entfalten, können wir das Signal daher einfach mit dieser Näherung wieder falten:

   [0.0208   -0.1640    1.2910   -0.1640    0.0208]

Dieser Kernel sieht interessant aus! Es ist ein Schärfoperator. Unsere Intuition ist richtig, das Schärfen hebt Unschärfe auf.

Wenn Sie zufälliges Rauschen hinzugefügt haben, können Sie das ursprüngliche Signal nicht erhalten ... Sie können versuchen, die Signale im Frequenzbereich zu trennen (wenn das Rauschen und das Signal unterschiedliche Frequenzen haben). Aber es scheint, dass Sie nach einem Wiener Filter suchen .

Ich denke, das ist immer noch ein offenes Problem.

Es gibt zahlreiche Forschungsarbeiten, die versuchen, das ursprüngliche Signal so gut wie möglich wiederherzustellen.

Ein klassischer Ansatz sind Wavelet-basierte Methoden .

Es gibt auch Wörterbuchansätze wie diesen .

Sie können sich einen tieferen Überblick über das Problem verschaffen, indem Sie den Untersuchungen von David L. Donho, Michael Elad, Alfred M. Bruckstein usw. folgen.

Wenn ich das Problem richtig verstanden habe, können wir das Problem wie folgt formalisieren:

Wir haben ein Signalmodell,

$y = Hx + \eta$

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

$\eta$

Ich habe nicht an der Signalwiederherstellung unter dem Poisson-Rauschen gearbeitet, aber ich habe gegoogelt und festgestellt, dass dieses Papier nützlich sein kann. Ähnliche Ansätze in diesem Zusammenhang können für dieses Problem nützlich sein.

Die Entfaltung von verrauschten Daten ist bekanntermaßen ein schlecht gestelltes Problem, da das Rauschen im rekonstruierten Signal willkürlich vergrößert wird. Daher ist ein Regularisierungsverfahren erforderlich, um die Lösung zu stabilisieren. Hier finden Sie ein MATLAB-Paket, das dieses Problem durch Implementierung des Regularisierungsalgorithmus von Tikhonov behebt:

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

Ich werde ganz am Anfang der Frage stehen. In MATLAB gibt es Entfaltungsfunktionen, die für Bildverarbeitungsanwendungen verwendet werden. Sie können diese Funktionen jedoch auch für 1D-Signale verwenden. Zum Beispiel,

% 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) Warum dann nicht das Ausgangssignal mit der hFunktion auflösen und das (fast) Eingangssignal erhalten? Ich benutze hier die Wiener Entfaltung

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

Wenn Sie die hFunktion nicht kennen , aber den Ein- und Ausgang kennen, können Sie dieses Mal das Eingangssignal mit dem Ausgang dekonvolvieren, der die h^-1Funktion liefert . Dann können Sie es als Filter verwenden, um das verrauschte Signal zu filtern. ( 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')

Ich hoffe, es hilft.

Faltung ist die Multiplikation und Summierung von zwei Signalen untereinander. Ich spreche von zwei deterministischen Signalen. Wenn Sie sich voneinander lösen wollen, entspricht dies der Lösung des Gleichungssystems. Wie Sie vielleicht wissen, sind Gleichungssysteme nicht immer lösbar. Das Gleichungssystem kann überbestimmt, unterbestimmt oder genau lösbar sein.

Wenn Sie etwas Rauschen hinzufügen, verlieren Sie einige Informationen und können diese Informationen nicht zurückerhalten. Sie können das lineare Gleichungssystem erneut lösen, indem Sie berücksichtigen, dass jedem Koeffizienten ein Rauschausdruck hinzugefügt wird. Wie Sie in einer anderen Antwort auf Ihre Frage sehen können, möchten Sie möglicherweise zuerst das ursprüngliche Signal aus dem verrauschten Signal abschätzen und dann versuchen, das Gleichungssystem zu lösen.

Es ist wichtig zu beachten, dass das Rauschen zu den multiplizierten und summierten Koeffizienten addiert wird. Daher kann es sogar vorkommen, dass Ihr Gleichungssystem letztendlich nicht eindeutig lösbar ist. Um sicherzugehen, dass es eindeutig lösbar ist, sollte Ihre Koeffizientenmatrix quadratisch und von vollem Rang sein.

Dies wäre schwierig zu tun. Die Faltung mit einem Gaußschen entspricht der Multiplikation mit einer Fourier-Transformation des Gaußschen im Frequenzbereich. Dies ist zufällig auch ein Gaußscher Wert, der im Wesentlichen ein Tiefpassfilter und ein wirklich effektiver Filter ist. Sobald Sie Rauschen hinzufügen, werden alle Informationen, die sich im "Stoppband" des Gaußschen befinden, zerstört. Es gibt keine Möglichkeit, das wiederherzustellen.

Die Entfaltung multipliziert sich im wesentlichen mit der Umkehrung des Frequenzgangs. Hier ist das Problem: Die Umkehrung des Frequenzgangs wird sehr, sehr groß, wenn der ursprüngliche Gaußsche Wert sehr klein ist. Bei diesen Frequenzen würden Sie das Rauschen grundsätzlich um große Beträge verstärken. Selbst wenn alles völlig rauschfrei wäre, würden Sie höchstwahrscheinlich auf numerische Probleme stoßen.

Ansätze

Es gibt viele Methoden zur Entfaltung (der Degradationsoperator ist nämlich linear und zeit- / rauminvariant).
Alle versuchen, mit der Tatsache umzugehen, dass das Problem in vielen Fällen schlecht ist.

Bessere Methoden sind solche, die dem Modell der wiederherzustellenden Daten eine gewisse Regularisierung hinzufügen.
Dies können statistische Modelle (Priors) oder beliebiges Wissen sein.
Für Bilder ist ein gutes Modell stückweise glatt oder spärlich in den Verläufen.

Um der Antwort willen wird jedoch ein einfacher parametrischer Ansatz gewählt: - Minimierung des Fehlers der kleinsten Quadrate zwischen den wiederhergestellten Daten im Modell und den Messungen.

Modell

Das Modell der kleinsten Quadrate ist einfach.
Die Zielfunktion als Funktion der Daten ist gegeben durch:

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

Das Optimierungsproblem ist gegeben durch:

$$\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$

Dies ist eine lineare Operation im endlichen Raum und kann daher unter Verwendung einer Matrixform geschrieben werden:

$$\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}$

Lösung

Die Lösung der kleinsten Quadrate ist gegeben durch:

$$\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)}$

Zustandsnummernanalyse

Was steckt hinter dieser Bedingungsnummer?
Man könnte es mit linearer Algebra beantworten.
Ein meiner Meinung nach intuitiverer Ansatz wäre es jedoch, im Frequenzbereich darüber nachzudenken.

Grundsätzlich dämpft der Degradationsoperator Energie mit im Allgemeinen hoher Frequenz.
Da dies in der Frequenz im Grunde genommen eine elementweise Multiplikation ist, würde man sagen, dass der einfache Weg, sie zu invertieren, die elementweise Division durch das inverse Filter ist.
Nun, es ist das, was oben gemacht wurde.
Das Problem tritt in Fällen auf, in denen der Filter die Energie praktisch auf Null abschwächt. Dann haben wir echte Probleme ...
Dies ist im Grunde das, was uns die Bedingungsnummer sagt, wie stark einige Frequenzen im Vergleich zu anderen gedämpft wurden.

Oben konnte man die Bedingungsnummer (unter Verwendung von [dB] -Einheiten) als Funktion des Gaussian Filter STD-Parameters sehen.
Wie erwartet, ist die Zustandszahl umso schlechter, je höher die STD ist, da eine höhere STD eine stärkere LPF bedeutet (Werte, die am Ende sinken, sind numerische Probleme).

Numerische Lösung

Das Ensemble des Gaußschen Weichzeichnungskerns wurde erstellt.

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

In MATLAB wurde das lineare System gelöst, pinv()das SVD-basiertes Pseudo Inverse und den \Operator verwendet.

Wie man sehen kann, ist die Lösung bei Verwendung der SVD erwartungsgemäß viel weniger empfindlich.

Warum liegt ein Fehler vor?
Suche nach einer Lösung (für die höchste sexuell übertragbare Krankheit):

Wie man sehen konnte, ist das Signal bis auf den Anfang und das Ende sehr gut wiederhergestellt.
Dies ist auf die Verwendung von Valid Convolution zurückzuführen, die wenig über diese Beispiele aussagt.

Lärm

Wenn wir Lärm hinzufügen würden, würden die Dinge anders aussehen!
Der Grund, warum die Ergebnisse zuvor gut waren, liegt in der Tatsache, dass MATLAB die DR der Daten verarbeiten und die Gleichungen lösen konnte, obwohl sie eine große Bedingungszahl hatten.

Eine große Bedingungszahl bedeutet jedoch, dass das inverse Filter einige Frequenzen stark verstärkt (um die starke Dämpfung umzukehren).
Wenn diese Geräusche enthalten, bedeutet dies, dass die Geräusche verstärkt werden und die Wiederherstellung schlecht ist.

Wie man oben sehen konnte, wird der Wiederaufbau jetzt nicht mehr funktionieren.

Zusammenfassung

Wenn man den Degradationsoperator genau kennt und das SNR sehr gut ist, funktionieren einfache Entfaltungsmethoden.
Das Hauptproblem der Entfaltung ist, wie stark der Degradationsoperator Frequenzen dämpft.
Je mehr es abschwächt, desto mehr SNR wird für die Wiederherstellung benötigt (dies ist im Grunde die Idee hinter Wiener Filter ).
Frequenzen, die auf Null gesetzt wurden, können nicht wiederhergestellt werden!

In der Praxis sollten einige Prioritäten hinzugefügt werden, um stabile Ergebnisse zu erzielen.

Der Code ist in meinem StackExchange Signal Processing Q2969 GitHub Repository verfügbar .

Im Allgemeinen besteht eine Methode zur Behandlung des Problems, das sich im Wesentlichen auf ein Problem der Extraktion von zwei oder mehr Komponenten verallgemeinert, darin, die Spektren G¹, G², G⋯ der Signale Nr. 1, Nr. 2, ..., Nr. N zu nehmen und die Gesamtsumme zu tabellieren Quadrat Γ (ν) = | G¹ (ν) | ² + | G² (ν) | ² + ⋯ + | Gⁿ (ν) | ² bei jeder Frequenz ν und Normalisierung von G₁ (ν) ≡ G¹ (ν) * / Γ (ν), G₂ (ν) ≡ G² (ν) * / Γ (ν), ..., G_n (ν) ≡ Gⁿ (ν) * / Γ (ν). Das Problem mit Unbestimmtheit und Rauschen entspricht der Tatsache, dass Γ (ν) ~ 0 für einige Frequenzen ν möglich ist. Fügen Sie dazu ein weiteres "Signal" hinzu, um G⁰ (ν) = Konstante zu extrahieren - das "Rausch" -Signal. Jetzt wird Γ (ν) unten streng begrenzt. Dies hängt mit ziemlicher Sicherheit mit der Regularisierung von Tichonow zusammen, aber ich habe nie ein Äquivalenzergebnis oder eine andere Entsprechung gefunden oder festgestellt. Es ist einfacher und direkter und intuitiver.

Alternativ können Sie die Gs als Vektoren behandeln, die mit einem geeigneten inneren Produkt ausgestattet sind, z. B. «G, G '» ≡ ≡ G (ν) * G' (ν) dν, und (G₀, G₁, ⋯, G_n) als Dual nehmen (zB die verallgemeinerte Umkehrung) von (G⁰, G¹, ⋯, Gⁿ) - vorausgesetzt natürlich, die Komponentenvektoren sind linear unabhängig.

Für die Gaußsche Entfaltung würde man n = 1, G⁰ = das "Rausch" -Signal und G¹ = das "Gaußsche" Signal einstellen.

Die Antwort von Andrey Rubshtein wird bei Vorhandensein von Rauschen kläglich scheitern, da das beschriebene Problem sehr empfindlich auf Rauschen und Modellierungsfehler reagiert. Es ist eine gute Idee, eine Faltungsmatrix zu konstruieren, aber die Verwendung der Regularisierung bei der Inversion ist bei einem solchen Problem ein absolutes Muss. Eine sehr einfache und unkomplizierte Regularisierungsmethode (obwohl rechenintensiv) ist die TSVD ( Truncated Singular Value Decomposition ). Methoden wie Tikhonov Regularization und Total Variation Regularizationsind einen Besuch wert. Die Tikhonov-Regularisierung (und ihre allgemeine Form) hat eine sehr elegante gestapelte Form, die in Matlab einfach zu implementieren ist. Lesen Sie das Buch: Lineare und nichtlineare inverse Probleme mit praktischen Anwendungen von Samuli Siltanen und Jennifer Mueller.

Eigentlich ist die Frage nicht klar. Aber die Antworten haben bestätigt, wonach Sie gefragt haben. Sie können ein System linearer algebraischer Gleichungen erstellen, wie einige Leute raten, das ist richtig, aber die Matrix, die auf einem bekannten Signal basiert, ist so genannt, schlecht konditioniert. Das heißt, wenn Sie versuchen, es zu invertieren, töten die Kürzungsfehler die Lösung und Sie erhalten Zufallszahlen als Ergebnis. Der übliche Ansatz ist ein eingeschränktes Extremum. Sie minimieren die Norm der Lösung || x || mit Einschränkung || Ax - y || <Delta. Sie suchen also nach einem solchen x mit der kleinsten Norm, bei dem der Unterschied zwischen Ax und y nicht groß sein darf. Es ist sehr einfach, der Hauptdiagonale der Matrix, die bei Anwendung der kleinsten Quadrate erhalten wird, einen sogenannten Regularisierungsparameter hinzuzufügen. Es heißt Tikhonov-Regularisierung. Ich habe Codierungsbeispiele, die das tun,