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 );
void VF_convolvewEdit( fVector Y, fVector Flt, fVector X, fVector Rsp, ui size; fComplex thresh );
C++ VecObj#include <OptiVec.h>
void vector<T>::convolve( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp );
void vector<T>::convolvewEdit( vector<T> Flt, const vector<T>& X, const vector<T>& Rsp, complex<T> thresh );
Pascal/Delphiuses VFstd;
procedure VF_convolve( Y, Flt, X, Rsp:fVector; size:UIntSize );
procedure VF_convolvewEdit( Y, Flt, X, Rsp:fVector; size:UIntSize; thresh:fComplex );
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 );
int cudaVF_convolvewEdit( fVector d_Y, fVector d_Flt, fVector d_X, fVector d_Rsp, ui size, fComplex thresh );
void VFcu_convolvewEdit( fVector h_Y, fVector h_Flt, fVector h_X, fVector h_Rsp, ui size, fComplex thresh );
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 );
function cudaVF_convolvewEdit( d_Y, d_Flt, d_X, d_Rsp:fVector; size:UIntSize; thresh:fComplex ): IntBool;
procedure VFcu_convolvewEdit( h_Y, h_Flt, h_X, h_Rsp:fVector; size:UIntSize; thresh:fComplex );
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 muss so in Rsp übergeben werden, dass 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, dass Rsp dieselbe Anzahl von Elementen aufweisen muss wie X.

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

X, Y, Rsp und Flt müssen alle dieselbe Größe size besitzen; diese muss 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.

Eine Response-Funktion, bei der manche Frequenzen so stark gedämpft werden, dass für sie jegliche Information verlorengeht (also z.B. die hohen Frequenzen bei einem Tiefpass-Filter), erkennt man an sehr kleinen Werten von Flt für die betreffenden Frequenzen. Sehr klein" bedeutet dabei, dass sie relativ zu dem Spitzenwert im Bereich der Genauigkeit des jeweiligen Datentypes liegen, also in der Größenordnung von VF_absmax(Flt, size)*epsilon. Zur Minimierung von Rundungsfehlern bei der Faltung ersetzt VF_convolve solche sehr kleinen Werte in Flt durch 0. Der standardmäßig verwendete Schwellenwert hierfür kann durch VF_setRspEdit ausgelesen werden. Um für sämtliche Aufrufe von VF_convolve und VF_deconvolve einen anderen Schwellenwert einzustellen, kann die Funktion VF_setRspEdit aufgerufen werden. Diese Methode ist aber nicht fiber-sicher und daher nicht geeignet, um den Schwellenwert bei verschiedenen Aufrufen von VF_convolve bzw. VF_deconvolve unterschiedlich einzustellen. Hier muss die Variante VF_convolvewEdit verwendet werden, die den gewünschten Schwellenwert als Argument thresh übernimmt (und die standardmäßig eingestellte Schwelle ignoriert). Da Flt aus komplexen Zahlen besteht und es gelegentlich wünschenswert ist, Real- und Imaginärteile unterschiedlich zu behandeln, ist thresh ebenfalls komplex.

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, dass 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, muss X gegebenenfalls an beiden Enden so extrapoliert werden, dass 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 muss 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, dass Flt als fVector und nicht als cfVector deklariert ist, obwohl die in Flt gespeicherte Information aus komplexen Zahlen besteht. Der Grund ist, dass 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