VF_convolveVD_convolveVE_convolve
FunktionFaltung (Konvolution) mit einer Impuls-Antwortfunktion
Syntax C/C++#include <VFstd.h>
void VF_convolve( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size );
C++ VecObj#include <OptiVec.h>
void vector<T>::convolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp );
Pascal/Delphiuses VFstd;
procedure VF_convolve( Y, Flt, X, Rsp:fVector; size:UIntSize );
CUDA-Funktion C/C++#include <cudaVFstd.h>
int cudaVF_convolve( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size );
void VFcu_convolve( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size );
CUDA-Funktion Pascal/Delphiuses VFstd;
function cudaVF_convolve( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize ): IntBool;
procedure VFcu_convolve( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize );
BeschreibungDer Vektor X wird mit der Impulsantwortfunktion (engl: response function) Rsp gefaltet und das Ergebnis in Y gespeichert. Ein Filter Flt wird gleichzeitig berechnet. Sollen mehr als nur ein Vektor mit derselben Impulsantwortfunktion gefaltet werden, so benutze man VF_convolve lediglich für den ersten. Für die weiteren rufe man VF_filter mit dem durch VF_convolve berechneten Filter.
Die Response-Funktion muß so in Rsp übergeben werden, daß der Nullpunkt auf dem nullten Element liegt. Die Elemente für positive Zeiten (oder was immer die unabhängige Variable ist) folgen als Rsp1 bis Rspsize/2 und die Elemente für negative Zeiten, beginnend mit der am stärksten negativen Zeit, als Rspsize/2+1 bis Rspsize-1. Um diese Ordnung zu erreichen, bediene man sich gegebenenfalls der Funktionen VF_rotate oder VF_reflect.
Man beachte, daß Rsp dieselbe Anzahl von Elementen aufweisen muß wie X.

Das Ergebnis der Faltung erscheint skaliert mit der Summe aller Elemente von Rsp. Um eine Amplituden-neutrale Faltung durchzuführen, muß Rsp also auf 1.0 normiert sein.

X, Y, Rsp und Flt müssen alle dieselbe Größe size besitzen; diese muß eine ganzzahlige Potenz von 2 sein. X darf durch Y überschrieben werden und Rsp durch Flt, aber X und Flt sowie Y und Rsp müssen voneinander verschieden sein.

Die Behandlung von Rundungsfehlern in der Konstruktion von Flt aus Rsp kann über die Funktion VF_setRspEdit modifiziert werden.

Beispiel C/C++VF_ramp( Time, 1024, 0.0, 1.0 );
VF_Gauss( Rsp, Time, 513, 10.0, 0.0, 1.0 );
      /* Response-Funktion für positive Zeiten ab 0 */
VF_reflect( Rsp+1, 1023 );
      /* ... und für negative Zeiten */
VF_divC( Rsp, Rsp, 1024, VF_sum( Rsp, 1024 ) );
      /* Normierung */
VF_convolve( X, Rsp, X, Rsp, 1024 );
      /* Faltung; X wird durch das gesuchte Ergebnis über-
         schrieben und Rsp durch das Frequenzfilter */
VF_filter( Y, Y, Rsp, 1024 );
      /* nächste Faltung: Anstatt eines erneuten Aufrufs von
         VF_convolve wird Y mit dem eben aus Rsp gewonnenen
         Frequenzfilter gefiltert */
 
Beispiel Pascal/DelphiVF_ramp( Time, 1024, 0.0, 1.0 );
VF_Gauss( Rsp, Time, 513, 10.0, 0.0, 1.0 );
      (* Response-Funktion für positive Zeiten ab 0 *)
VF_reflect( VF_Pelement(Rsp,1), 1023 );
      (* ... und für negative Zeiten *)
VF_divC( Rsp, Rsp, 1024, VF_sum( Rsp, 1024 ) );
      (* Normierung *)
VF_convolve( X, Rsp, X, Rsp, 1024 );
      (* Faltung; X wird durch das gesuchte Ergebnis über-
         schrieben und Rsp durch das Frequenzfilter *)
VF_filter( Y, Y, Rsp, 1024 );
      (* nächste Faltung: Anstatt eines erneuten Aufrufs von
         VF_convolve wird Y mit dem eben aus Rsp gewonnenen
         Frequenzfilter gefiltert  *)

Mathematisch basiert diese Faltung auf der Annahme, daß X periodisch ist. Sie funktioniert auch dann noch gut, wenn X zwar nicht periodisch ist, aber an beiden Enden zu demselben Wert konvergiert (also X0 = Xsize-1). Ist dies nicht der Fall, so beeinflussen sich die Elemente an beiden Enden gegenseitig, und es kommt zu Rand-Effekten. Um diese zu vermeiden, muß X gegebenenfalls an beiden Enden so extrapoliert werden, daß sich der ursprüngliche Vektor X in der Mitte eines größeren Vektors wiederfindet. Es sind mindestens so viele Elemente an beiden Seiten anzufügen, wie der halben Breite der Impulsantwortfunktion entspricht (bei asymmetrischen Antwortfunktionen nehme man die breitere Flanke als Maß). Die Faltung ist dann für den größeren Vektor durchzuführen. Die gesuchte Faltung von X mit Rsp befindet sich in der Mitte des so erhaltenen Ergebnis-Vektors.

Sind die Werte an beiden Enden von X zwar nicht gleich, werden aber jeweils durch genügend "sanfte" Konvergenz erreicht, so kann der Unterschied zwischen den Endpunkten auch als linearer Trend betrachtet werden. Man muß diesen dann nur für die Faltung beseitigen und zum erhaltenen Ergebnis wieder hinzuaddieren.

Beispiel C/C++d = (X[size-1] - X[0]) / (size-1);
VF_ramp( Trend, size, 0.0, d );
VF_subV( Y, X, Trend, size );
VF_convolve( Y, Flt, Y, Rsp, size );
VF_addV( Y, Y, Trend, size );
Beispiel für die Vermeidung von Randeffekten in Pascal/Delphi d := (VF_element(X,size-1) - X^) / (size-1);
VF_ramp( Trend, size, 0.0, d );
VF_subV( Y, X, Trend, size );
VF_convolve( Y, Flt, Y, Rsp, size );
VF_addV( Y, Y, Trend, size );

Man wird bemerken, daß Flt als fVector und nicht als cfVector deklariert ist, obwohl die in Flt gespeicherte Information aus komplexen Zahlen besteht. Der Grund ist, daß diese Zahlen in dem bei VF_FFT beschriebenen gepackten Format vorliegen, das ausschließlich im Zusammenhang mit der Fourier-Transformation reeller Vektoren angewandt wird.

FehlerbehandlungWenn size nicht eine Potenz von 2 ist, meldet sich VF_FFT (worauf VF_convolve basiert) mit der Fehlermeldung "Size must be integer power of 2." und bricht das Programm ab.
Rückgabewertkeiner
QuerverweisVF_filter,   VF_deconvolve,   VF_FFT,   VF_autocorr,   VF_xcorr,   VF_spectrum

VectorLib Inhaltsverzeichnis  OptiVec Home