Entfaltung (Dekonvolution) mit C++

Hallo Leute,

ich habe ein Messsignal, welches mit der Intrument Response Function (IRF) gefaltet ist.

Als Kleine Erklärung: Durch eine Faser wird ein kurzer Lichtimpuls auf eine Fläche "geschossen", welche das Licht aufnimmt (absorbiert) und kurz darauf wieder abgiebt und zurückreflektiert. Das zurückreflecktierende Licht wird erneut durch eine Faser zurück zur Elektronik geleitet. Mit dem Senden des Lichtimpulses wird eine Stoppuhr gestartet. Die Stoppuhr wird wieder gestoppt, sobald das Licht durch die zweite Faser bei der Elektronik ankommt.

Über die Zeit erhält man also eine Kurve, die vereinfacht bei Null beginnt, dann sehr stark aufsteigt, sobald das Licht von der Fläche reflektiert wird und dann langsam abflacht, bis die Fläche nichts mehr reflektiert.

Ich will nun messen, wie viel Zeit vergeht zwischen Absorption des Lichtes auf der Fläche und der Rückreflektion (es handelt sich dabei um wenige Nanosekunden, sehr genaue Stoppuhr ). Nun vergeht natürlich Zeit, in welcher das Licht durch die Fasern geleitet und in der Elektronik verarbeitet wird und welche das Ergebnis verfälschen. Dieser Weg bzw. diese Beeinflussung durch Fasern und Elektronik nennt man IRF und kann man mit einem Spiegel messen. Das Ergebnis ist eine Funktion, eben die besagt Instrument Response Function.

So, ich will Euch mit dem ganzen Käse nicht nerven, aber es geht darum, dass ihr nachvollziehen könnt, wo mein Problem ist.
Ich habe also das gemessene Gesamtergebnis und ich habe die gemessene IRF.
Das Gesamtergebnis ist mit der IRF gefaltet.

Diese beiden Kurven/Funktionen müsste ich nun "Entfalten". Ich bin mir nicht sicher, ob es dafür bei boost etwas gibt und bin mir auch mathematisch trotz recherchen nicht ganz sicher und bitte um Eure Hilfe, ob jemand soetwas schon einmal gemacht hat. Ich brauche das für die Uni und hänge in wenig fest.

Falls ich zu verwirrend erklärt habe, tut mir das leid, fragt einfach nach, dann versuche ich klarer zu werden. Im Prinzip geht es einfach um zwei gefaltete Funktionen die gemeinsam eine weitere Funktion bilden, wobei eine bekannt ist und ich die andere errechnen muss.

knivil

Faltung wird im Fourierraum zur Multiplikation. Aber ich weiss nicht, ob man das machen darf. Dann gibt es noch Wiener deconvolution, habe ich aber noch nicht ausprobiert.

krümelkacker

Nö, Boost hat sowas sicherlich nicht. Jedenfalls nicht, dass ich wüsste. Das ist schon etwas sehr Spezielles aus dem Signalverarbeitungsbereich.

Ich würde das folgendermaßen machen
- Impulsantwort des inversen Filters bestimmen
- Dein Messsignal mit dieser Impulsantwort falten, ggf per "Schnelle Faltung" (siehe dspguide.com für die schnelle FFT-basierte Faltung).

Die Impulsantwort des inversen Filters bekommst du, indem du deine IRF mit Nullen links und rechts auffüllst (zero adding) und entsprechend so zyklisch verschiebst, dass der Nullpunkt der Zeitachse dem Element mit Index 0 entspricht (also das erste), dann darauf eine FFT machst, dann vom Ergebnis komponentenweise die Kehrwerte bestimmst, dann darauf die iFFT anwendest, dann es wieder so zyklisch zurückschiebst, dass der Nullpunkt in der Mitte liegt und dann bei dem, was da rauskommt noch eine nette Fensterfunktion draufmultiplizierst.

Das mit dem Zero Padding ist einigermaßen wichtig, weil du damit den sogenannten "wrap-around"-Fehler vermeidest. Am Ende ist wahrscheinlich nur der mittlere Ausschnitt einigermaßen von Interesse, daher das Fenstern.

Hinterher am besten "die Probe" machen und die IRF mit der neu berechneten Impulsantwort für den inversen Filter falten und gucken, ob da auch etwas rauskommt, was genügend nah an einem "Unit-Impuls" liegt -- sozusagen a*b=1 wobei b das Inverse von a sein soll.

Wenn du dir die Beträge der FFT-Koeffizienten anguckst, kannst du damit auch die Konditionszahl berechnen, die zum äquivalenten Gleichungssystem gehört, welches dein Faltungsproblem auch lösen würde. Die Konditionszahl kann je nach IRF auch leider unendlich sein und zwar dann, wenn mindestens eine Frequenz durch die Faltung per IRF komplett gedämpft wird. Das kannst du natürlich nicht mehr wiederherstellen. Einmal die Amplitude mit 0 multipliziert und du kommst nicht mehr drauf, was die Amplitude vorher mal war, quasi.

Matlab oder Octave bietet sich für sowas eigentlich prima an. Da kann man dann erstmal ein bisschen rumspielen, sich Dinge plotten lassen, bevor man das dann in C++ gießt (wenn überhaupt).