Wiener-Dekonvolution

Von links: Originalbild, unscharfes Bild, mittels Wiener-Entfaltung geschärftes Bild

In der Mathematik stellt die Wiener-Dekonvolution eine Anwendung des Wiener-Filters für Rauschprobleme in der Entfaltung dar. Sie versucht, bei der Entfaltung den Einfluss von Rauschen im Frequenzraum zu minimieren und wird daher meist bei schlechten Signal-Rausch-Verhältnissen angewendet.

Die Wiener-Entfaltung ist in Entfaltungsanwendungen im Fotobereich weit verbreitet, da das Frequenzspektrum von Bildern im sichtbaren Bereich vergleichsweise einfach zu bestimmen ist.

Die Wiener-Entfaltung ist nach Norbert Wiener benannt.

Definition

Sei

y ( t ) = h ( t ) x ( t ) + n ( t ) {\displaystyle y(t)=h(t)*x(t)+n(t)} ,

wobei {\displaystyle *} die Faltung bezeichnet und

  • x ( t ) {\displaystyle x(t)} das (unbekannte) Eingangssignal zur Zeit t {\displaystyle t} .
  • h ( t ) {\displaystyle h(t)} die bekannte Impulsantwort eines linear zeitinvarianten Systems
  • n ( t ) {\displaystyle n(t)} ein unbekanntes Rauschen, das unabhängig von x ( t ) {\displaystyle x(t)} ist
  • y ( t ) {\displaystyle y(t)} das beobachtete Signal.

Das Ziel ist, g ( t ) {\displaystyle g(t)} zu bestimmen, sodass sich x ( t ) {\displaystyle x(t)} wie folgt ergibt:

x ^ ( t ) = g ( t ) y ( t ) {\displaystyle {\hat {x}}(t)=g(t)*y(t)}

wobei x ^ ( t ) {\displaystyle {\hat {x}}(t)} eine Abschätzung von x ( t ) {\displaystyle x(t)} mit minimiertem quadratischen Fehler ist.

Der Wiener-Filter liefert ein solches   g ( t ) {\displaystyle \ g(t)} . Er kann am einfachsten im Frequenzraum beschrieben werden:

G ( f ) = H ( f ) S ( f ) | H ( f ) | 2 S ( f ) + N ( f ) {\displaystyle G(f)={\frac {H^{*}(f)S(f)}{|H(f)|^{2}S(f)+N(f)}}}

wobei

  • G ( f ) {\displaystyle G(f)} und H ( f ) {\displaystyle H(f)} die Fourier-Transformation von   g {\displaystyle \ g} bzw. h {\displaystyle h} bei der Frequenz f {\displaystyle f} sind.
  • S ( f ) {\displaystyle S(f)} die mittlere Spektrale Leistungsdichte des Eingangssignals x ( t ) {\displaystyle x(t)} ist
  • N ( f ) {\displaystyle N(f)} die mittlere Spektrale Leistungsdichte des Rauschens n ( t ) {\displaystyle n(t)} ist
  • H ( f ) {\displaystyle H^{*}(f)} die komplex Konjugierte von H ( f ) {\displaystyle H(f)} bezeichnet.

Die Filter-Operation kann, wie oben, im Zeitbereich, oder im Frequenzraum durchgeführt werden:

X ^ ( f ) = G ( f ) Y ( f ) {\displaystyle {\hat {X}}(f)=G(f)Y(f)}

(wobei X ^ ( f ) {\displaystyle {\hat {X}}(f)} die Fourier-Transformation von x ^ ( t ) {\displaystyle {\hat {x}}(t)} ) ist. Eine inverse Fourier-Transformation von X ^ ( f ) {\displaystyle {\hat {X}}(f)} liefert x ^ ( t ) {\displaystyle {\hat {x}}(t)} .

Es ist zu beachten, dass bei Bildern die Argumente t {\displaystyle t} und f {\displaystyle f} zweidimensional werden; Das Ergebnis bleibt aber das gleiche.

Interpretation

Die Anwendung des Wiener-Filters zeigt sich, wenn die obige Gleichung umgeschrieben wird:

G ( f ) = 1 H ( f ) [ | H ( f ) | 2 | H ( f ) | 2 + N ( f ) S ( f ) ] = 1 H ( f ) [ | H ( f ) | 2 | H ( f ) | 2 + 1 S N R ( f ) ] {\displaystyle {\begin{aligned}G(f)&={\frac {1}{H(f)}}\left[{\frac {|H(f)|^{2}}{|H(f)|^{2}+{\frac {N(f)}{S(f)}}}}\right]\\&={\frac {1}{H(f)}}\left[{\frac {|H(f)|^{2}}{|H(f)|^{2}+{\frac {1}{\mathrm {SNR} (f)}}}}\right]\end{aligned}}}

Hierbei ist   1 / H ( f ) {\displaystyle \ 1/H(f)} das Inverse des Ausgangssystems und   S N R ( f ) = S ( f ) / N ( f ) {\displaystyle \ \mathrm {SNR} (f)=S(f)/N(f)} ist das Signal-Rausch-Verhältnis. Ohne Rauschen (d. h. unendliches Signal-zu-Rausch-Verhältnis) ist der Term innerhalb der eckigen Klammern gleich 1, was bedeutet, dass der Wiener-Filter einfach das Inverse des Systems ist, wie man es erwarten kann. Wenn das Rauschen bei bestimmten Frequenzen steigt, das Signal-zu-Rausch-Verhältnis also fällt, nimmt der Term innerhalb der eckigen Klammern ebenfalls ab. Das heißt, der Wiener-Filter dämpft die Frequenzen in Abhängigkeit von ihrem Signal-zu-Rausch-Verhältnis.

Die obige Gleichung setzt voraus, dass der spektrale Inhalt eines typischen Bildes sowie der des Rauschens bekannt ist. Meistens sind die beiden Größen nicht bekannt, können aber abgeschätzt werden. Zum Beispiel hat bei Fotos das Signal (das originale Bild) typischerweise starke Anteile von niedrigen und schwache Anteile hohen Frequenzen und die Rauschanteile verteilen sich gleichmäßig über alle Frequenzen.

Herleitung

Wie oben beschrieben soll eine Annäherung an das Originalbild erzeugt werden das den quadratischen Fehler minimiert. Dieser lässt sich durch

  ϵ ( f ) = E | X ( f ) X ^ ( f ) | 2 {\displaystyle \ \epsilon (f)=\mathbb {E} \left|X(f)-{\hat {X}}(f)\right|^{2}}

ausdrücken, wobei   E {\displaystyle \ \mathbb {E} } der Erwartungswertoperator ist.

Wird   X ^ ( f ) {\displaystyle \ {\hat {X}}(f)} ersetzt, lässt sich der Ausdruck umschreiben:

ϵ ( f ) = E | X ( f ) G ( f ) Y ( f ) | 2 = E | X ( f ) G ( f ) [ H ( f ) X ( f ) + V ( f ) ] | 2 = E | [ 1 G ( f ) H ( f ) ] X ( f ) G ( f ) V ( f ) | 2 {\displaystyle {\begin{aligned}\epsilon (f)&=\mathbb {E} \left|X(f)-G(f)Y(f)\right|^{2}\\&=\mathbb {E} \left|X(f)-G(f)\left[H(f)X(f)+V(f)\right]\right|^{2}\\&=\mathbb {E} {\big |}\left[1-G(f)H(f)\right]X(f)-G(f)V(f){\big |}^{2}\end{aligned}}}

Das Quadrat kann entwickelt werden und ergibt:

ϵ ( f ) = [ 1 G ( f ) H ( f ) ] [ 1 G ( f ) H ( f ) ] E | X ( f ) | 2 [ 1 G ( f ) H ( f ) ] G ( f ) E { X ( f ) V ( f ) } G ( f ) [ 1 G ( f ) H ( f ) ] E { V ( f ) X ( f ) } + G ( f ) G ( f ) E | V ( f ) | 2 {\displaystyle {\begin{aligned}\epsilon (f)&={\Big [}1-G(f)H(f){\Big ]}{\Big [}1-G(f)H(f){\Big ]}^{*}\,\mathbb {E} |X(f)|^{2}\\&{}-{\Big [}1-G(f)H(f){\Big ]}G^{*}(f)\,\mathbb {E} {\Big \{}X(f)V^{*}(f){\Big \}}\\&{}-G(f){\Big [}1-G(f)H(f){\Big ]}^{*}\,\mathbb {E} {\Big \{}V(f)X^{*}(f){\Big \}}\\&{}+G(f)G^{*}(f)\,\mathbb {E} |V(f)|^{2}\end{aligned}}}

Allerdings wird angenommen, dass das Rauschen unabhängig vom Signal ist, also:

  E { X ( f ) V ( f ) } = E { V ( f ) X ( f ) } = 0 {\displaystyle \ \mathbb {E} {\Big \{}X(f)V^{*}(f){\Big \}}=\mathbb {E} {\Big \{}V(f)X^{*}(f){\Big \}}=0}

Die spektrale Leistungsdichte wird definiert als:

  S ( f ) = E | X ( f ) | 2 {\displaystyle \ S(f)=\mathbb {E} |X(f)|^{2}}
  N ( f ) = E | V ( f ) | 2 {\displaystyle \ N(f)=\mathbb {E} |V(f)|^{2}}

Damit ergibt sich:

ϵ ( f ) = [ 1 G ( f ) H ( f ) ] [ 1 G ( f ) H ( f ) ] S ( f ) + G ( f ) G ( f ) N ( f ) {\displaystyle {\begin{aligned}\epsilon (f)&={\Big [}1-G(f)H(f){\Big ]}{\Big [}1-G(f)H(f){\Big ]}^{*}S(f)\\&{}+G(f)G^{*}(f)N(f)\end{aligned}}}

Um den minimalen Fehler zu finden wird eine   G ( f ) {\displaystyle \ G(f)} differenziert und gleich null gesetzt. Da das einen komplexen Wert liefert, ist   G ( f ) {\displaystyle \ G^{*}(f)} eine Konstante.

  d ϵ ( f ) d G ( f ) = G ( f ) N ( f ) H ( f ) [ 1 G ( f ) H ( f ) ] S ( f ) = 0 {\displaystyle \ {\frac {d\epsilon (f)}{dG(f)}}=G^{*}(f)N(f)-H(f){\Big [}1-G(f)H(f){\Big ]}^{*}S(f)=0}

Diese Gleichung kann umgeschrieben werden, um den Wiener-Filter zu erhalten.

Referenzen

  • Rafael Gonzalez, Richard Woods, and Steven Eddins: Digital Image Processing Using Matlab. Prentice Hall, 2003.
  • Vergleich verschiedener Entfaltungsmethoden (englisch).
  • Entfaltung mit einem Wiener-Filter (englisch)
  • Beispiel einer auf Bewegungsunschärfe angewandte Wiener-Entfaltung (Quellcode in MATLAB/GNU Octave)