ADD: convolve

Author: PascalCorpsman

data_control/Vectormath/testcase1.pas

@@ -66,6 +66,15 @@
     Procedure LoadFromSaveToStream;
   End;
 
+  { TOtherTests }
+
+  TOtherTests = Class(TTestCase)
+  protected
+  published
+    Procedure Convolve1D;
+    Procedure ConvolveAVG;
+  End;
+
 Implementation
 
 { TMatrix4x4Tests }
@@ -421,6 +430,47 @@
   m.free;
 End;
 
+{ TOtherTests }
+
+Procedure TOtherTests.Convolve1D;
+Var
+  a, b, c: TVectorN;
+Begin
+  (*
+   * Falten mit "nicht symmetrischem" filter kern
+   *)
+  a := VN([1, 2, 3]);
+  b := VN([4, 5, 6]);
+  c := Convolve(a, b);
+  AssertEquals('Invalid len', 5, length(c));
+  AssertEquals('Invalid result [0]', 4, c[0]);
+  AssertEquals('Invalid result [1]', 13, c[1]);
+  AssertEquals('Invalid result [2]', 28, c[2]);
+  AssertEquals('Invalid result [3]', 27, c[3]);
+  AssertEquals('Invalid result [4]', 18, c[4]);
+End;
+
+Procedure TOtherTests.ConvolveAVG;
+Var
+  a, b, c: TVectorN;
+Begin
+  (*
+   * Faltung als Mittelwert filter
+   *)
+  a := vn([0, 0, 10, 10, 10, 0, 0]);
+  b := vn([0.5, 0.5]);
+  c := Convolve(a, b);
+  AssertEquals('Invalid len', 8, length(c));
+  AssertEquals('Invalid result [0]', 0, c[0]);
+  AssertEquals('Invalid result [1]', 0, c[1]);
+  AssertEquals('Invalid result [2]', 5, c[2]);
+  AssertEquals('Invalid result [3]', 10, c[3]);
+  AssertEquals('Invalid result [4]', 10, c[4]);
+  AssertEquals('Invalid result [5]', 5, c[5]);
+  AssertEquals('Invalid result [6]', 0, c[6]);
+  AssertEquals('Invalid result [7]', 0, c[7]);
+End;
+
 Initialization
   Randomize;
   DefaultFormatSettings.DecimalSeparator := '.';
@@ -429,5 +479,7 @@
   RegisterTest(TMatrix3x3Tests);
   RegisterTest(TMatrix4x4Tests);
   RegisterTest(TMatrixNxMTests);
+  RegisterTest(TOtherTests);
+
 End.
 

data_control/uvectormath.pas

@@ -1,7 +1,7 @@
 (******************************************************************************)
 (* uvectormat.pas                                                  12.11.2014 *)
 (*                                                                            *)
-(* Version     : 0.18                                                         *)
+(* Version     : 0.19                                                         *)
 (*                                                                            *)
 (* Author      : Uwe Schächterle (Corpsman)                                   *)
 (*                                                                            *)
@@ -60,6 +60,7 @@
 (*                     TMatrix3x3.getInverse, TMatrix2x2.getInverse           *)
 (*                    InvertMatrix2 for TMatrix4x4, TMatrix2x2                *)
 (*               0.18 CalculatePlumbFootPoint                                 *)
+(*               0.19 Convolve                                                *)
 (*                                                                            *)
 (******************************************************************************)
 Unit uvectormath;
@@ -580,6 +581,14 @@
 // Die Koeffizienten von a werden als Vector übergben
 Function CalculatePolynom(x: TBaseType; Const a: TVectorN): TBaseType;
 
+// Faltet 2 Vectoren mit einander und gibt das Ergebnis zurück, Achtung,
+// hier wird die "ineffiziente" Quadratische Faltung verwendet,
+// also die Vectoren lieber nicht zu groß machen!
+// Es gibt auch noch die Variante, bei der man zuerst die beiden Vektoren
+// durch eine FFT jagt, dann komponentenweise multipliziert und dann wieder
+// zurück transformiert, das ist bei Großen Vektoren tatsächlich schneller (siehe hier: https://www.youtube.com/watch?v=KuXjwB4LzSA )
+Function Convolve(Const a, b: tvectorN): TVectorN;
+
 Implementation
 
 {$IFDEF UseOperandOverloading}
@@ -3302,6 +3311,21 @@
   End;
 End;
 
+Function Convolve(Const a, b: tvectorN): TVectorN;
+Var
+  i, j: Integer;
+Begin
+  result := Nil;
+  setlength(result, length(a) + length(b) - 1);
+  For i := 0 To High(Result) Do Begin
+    Result[i] := 0;
+    For j := 0 To High(a) Do Begin
+      If (i - j >= 0) And (i - j <= High(b)) Then
+        Result[i] := Result[i] + a[j] * b[i - j];
+    End;
+  End;
+End;
+
 Function CalculateOrthoganlProjection(A, B, M: TVector2): TVector2;
 Var
   d_AB, d_M: TVector2;