Kreuzkorrelation und Signallängen

[Christoph a. (chrichri)]

Hi
ich habe zwei Signale (A und B). Signal B soll in A gefundne werden. 
Dabei gilt A ist länger (also mehr Messwerte) als B.

Dazu verwende ich in C++ die FFT-Library FFTW und berechne die 
Korrelation.

- A und B auf Länge länge(A)+länge(B)-1 bringen ("rechtsseitiges 
zeropadding")
- FFT von A und B bilden
- Beide FFTs elementweise multiplizieren, eines aber komplex konjugiert 
Ergebnis = FFT(A) * conj(FFT(B))
- Inverse FFT des Ergebnis IFFT(Ergebnis)
- Maxima in Inverse suchen. Maximaposition ist die gesuchte Position.

Das funktioniert soweit wie erwartet.

Mein Problem ist aber dass die Länge von A variablen ist (aber stehts 
größer/gleich der Länge von B) und ich so schnell wie möglich sein will.
D.h. ich will, wenn möglich, FFT(B) nur einmalig berechnen. Das Pronlem 
ist ja dass ich mit länge(A)+länge(B)-1 rechne.

Kann ich das irgendwie machen?
Gruß
Christoph

[Andreas S. (andreas) (Admin)]

"Schnelle Faltung" mit Overlap-Save oder Overlap-Add 
(http://inst.eecs.berkeley.edu/~ee123/fa11/docs/FastConv.pdf).

Was sind denn typische Längen von A und B?

[Christoph a. (chrichri)]

Danke.

A typisch 300-1200 messwerte
B typisch 100-300  messwerte

nur so eine größenordnung ;-)

[Florian F. (ultrazauberer)]

Wie sieht denn der Quellcode dazu aus? Finde ich ziemlich interessant, 
da ich etwas ähnliches vorhabe (Laufzeitmessung).
Vielleicht finde ich noch etwas dazu.

[Christoph a. (chrichri)]

Also das PDF von Florian konnte ich noch nicht bearbeiten, geschweige 
denn in Code umsetzen. Habe gerade auch andere Probleme an denen ich 
hänge, deswegen muss die Performance-Geschichte erstmal warten.

Meinen bisherigen Code darfst du aber gerne haben ;-)
Es ist aber nur ein Testprogramm! Ich verwende die Library FFTW für die 
Korrelation und die Library OpenCV weil ich später auch mit Bildern 
arbeite. (Deswegen is der Code auch umfangreicher als nötig)

Ich nehme ein Bild X und betrachte es Zeilenweise als Signal X_Zeile.
Dann nehme ich den festen Ausschnitt X_Template aus X_Zeile (z.B. Pixel 
100 bis 600) und korreliere X_Template mit X_Zeile (also z.B. 800 
einzelne Korrelationen bei einem Bild der Höhe 800 px). Für jede 
Korrelation wird das Maximim gesucht und als Ergebnis gespeichert. Ideal 
sollte überall 100 herauskommen (mit den Beipielwerten). (Tut es aber 
nicht)

Anschließend bilde suche ich den Median aller Ergebnisse. Er entspricht 
i.d.R. den gesuchten 100. (Wenn das Template zu ist dann nicht mehr)

Verbesserungswürdig ist die Art der Korrelation da nicht mittelwertfrei 
oder normiert. (Zum testen reicht erstmal die "normale" 
Kreuzkorrelation)

Was ich später nicht haben will ist dass ich jedes mal alle Daten 
(fft-daten etc) allokiere, die Datansätze fülle etc. Das kostet Zeit.

#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <opencv2/opencv.hpp>
#include <fftw3.h>
#include <fstream>
#include "fftw_samples.h"
#include <time.h>

void fft_1d();
int compare(const void * a, const void * b);

using namespace std;

int main (int argc, char * const argv[]) {
  
  fft_1d();
  
  return 0;
}

void fft_1d(){
  
  clock_t c_start, c_init, c_fft, c_log, c_tmp;
  
  int matchloc_should_be = 100+1; //+1 wegen matlab
  int width_of_template  = 600;
  cv::Mat image = cv::imread("abc.bmp",0);
  cv::Mat templ = cv::Mat(image,cv::Rect(matchloc_should_be,0,width_of_template,image.size().height)).clone();
  
  cv::Mat image_sm;
  cv::Mat templ_sm;
  pyrDown(image, image_sm);
  pyrDown(templ, templ_sm);
  imshow("image", image_sm);
  imshow("templ", templ_sm);
  cvWaitKey(0);
  cv::destroyAllWindows();
  
  printf("Matchloc should be: %d \n", matchloc_should_be);
  printf("image - width: %d height: %d \n", image.size().width, image.size().height);
  printf("templ - width: %d height: %d \n", templ.size().width, templ.size().height);
  
  
  c_tmp   = clock();
  c_start = clock();
  
  // length of data (both, image and template data must have same length)
  int size_w = image.size().width + templ.size().width - 1;
  
  // number of lines must be the same
  int size_h = image.size().height;
  
  //fftw plans
  fftw_plan    plan_forward, plan_backward;
  
  // Pointer to current data
  // fftw_complex    * curr_image_data, *cur_templ_data;
  
  // new array of pointers for image data
  fftw_complex    ** image_data = new fftw_complex*[size_h];
  fftw_complex    ** image_fft  = new fftw_complex*[size_h];
  
  // new array of pointers for template data
  fftw_complex    ** templ_data = new fftw_complex*[size_h];
  fftw_complex    ** templ_fft  = new fftw_complex*[size_h];
  
  // new array of pointers for result data
  fftw_complex    ** corr_data  = new fftw_complex*[size_h];
  fftw_complex    ** ifft_data  = new fftw_complex*[size_h];
  
  // allocate data for image and template
  for (int i = 0; i < size_h; i++) {
    image_data[i] = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
    image_fft[i]  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
    
    templ_data[i] = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
    templ_fft[i]  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
    
    corr_data[i]  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
    ifft_data[i]  = ( fftw_complex* ) fftw_malloc( sizeof( fftw_complex ) * size_w );
  }
  
  //zero padding image, template and result data
  for (int i = 0; i < size_h; i++) {
    for (int e = 0; e < size_w ; e++) {
      image_data[i][e][0] = 0.0;    // real      part
      image_data[i][e][1] = 0.0;    // imaginary  part
      
      templ_data[i][e][0] = 0.0;    // real      part
      templ_data[i][e][1] = 0.0;    // imaginary  part
      
      corr_data[i][e][0]  = 0.0;    // real      part
      ifft_data[i][e][1]  = 0.0;    // imaginary  part
    }
  }
  
  // populate image data
  for (int i = 0; i < size_h; i++) {
    for (int e = 0; e < image.size().width ; e++) {
      image_data[i][e][0] = (double)(image.at<uchar>(i,e));    // real      part
      //image_data[i][e][1] = 0.0;                // imaginary  part
    }
  }
  
  // populate template data
  for (int i = 0; i < size_h; i++) {
    for (int e = 0; e < templ.size().width ; e++) {      
      templ_data[i][e][0] =  (double)(templ.at<uchar>(i,e));;    // real      part
      //templ_data[i][e][1] = 0.0;                // imaginary  part
    }
  }
  
  c_init = clock() - c_tmp;
  c_tmp  = clock();
  
  // calculate fft of image and template data
  for (int i = 0; i < size_h; i++) {
    plan_forward  = fftw_plan_dft_1d(size_w, image_data[i], image_fft[i], FFTW_FORWARD, FFTW_ESTIMATE );
    fftw_execute( plan_forward );
    
    plan_forward  = fftw_plan_dft_1d(size_w, templ_data[i], templ_fft[i], FFTW_FORWARD, FFTW_ESTIMATE );
    fftw_execute( plan_forward );
  }
  
  //calculate correlation result in frequency domain
  double a,b,c,d;
  for (int i = 0; i < size_h; i++) {
    for (int e = 0; e < size_w; e++) {
      a = image_fft[i][e][0];  // real
      b = image_fft[i][e][1];  // imaginary
      c = templ_fft[i][e][0];  // real
      d = templ_fft[i][e][1];  // imaginary
      
      corr_data[i][e][0] = a*c + b*d;
      corr_data[i][e][1] = b*c - a*d;      
    }
  }
  
  //inverse transform result
  for (int i = 0; i < size_h; i++) {
    plan_backward = fftw_plan_dft_1d(size_w, corr_data[i], ifft_data[i], FFTW_BACKWARD, FFTW_ESTIMATE );
    fftw_execute( plan_backward );
  }
  
  //find maxima in result
  double * max = new double[size_h];
  int * maxpos = new int[size_h];
  for (int i = 0; i < size_h; i++) {
    max[i]    = 0;
    maxpos[i] = 0;
  }  
  for (int i = 0; i < size_h; i++) {
    for (int e = 0; e < size_w; e++) {
      if(ifft_data[i][e][0] > max[i]) {
        max[i]    = ifft_data[i][e][0];
        maxpos[i] = e;
      }
    }
  }
  
  c_fft = clock() - c_tmp;
  c_tmp = clock();
  
  // save to disk for matlab
  ofstream myfile;
  myfile.open ("log.txt");
  for (int e = 0; e < size_w; e++){
    for (int i = 0; i < size_h; i++) {
      myfile << ifft_data[i][e][0] << " ";
    }
    myfile << endl;
  }
  myfile.close();
  
  myfile.open ("log_corr_cpp.txt");
  for (int i = 0; i < size_h; i++){
    myfile << maxpos[i] << " " << max[i] << endl;
  }
  myfile.close();

  c_log = clock() - c_tmp;
  c_tmp = clock();
  
  // find median of maxima positions
  qsort((void*)maxpos, size_h, sizeof(int), compare);
  int median = 0;
  if(size_w%2){ // uneven
    median = maxpos[size_h/2 + 1];
  }
  else {
    median = (maxpos[size_h/2] + maxpos[size_h/2+1])/2;
  }
  printf("median found at: %d",median);
  
  // log sorted correlation result
  myfile.open ("log_corr_sorted_cpp.txt");
  for (int i = 0; i < size_h; i++){
    myfile << maxpos[i] << " " << max[i] << endl;
  }
  myfile.close();
  
  printf("\ntime: init %f - fft %f - log %f \n", 
       (float)c_init/(float)CLOCKS_PER_SEC, 
       (float)c_fft/(float)CLOCKS_PER_SEC,
       (float)c_log/(float)CLOCKS_PER_SEC);
}

int compare (const void * a, const void * b)
{
  return ( *(int*)a - *(int*)b );
}

[Florian F. (ultrazauberer)]

Danke für den Code! Scheint ja mehr als reichlich zu sein.

Da wir hier gerade beim Thema sind: Ich möchte ein einfaches Echolot / 
Sonar auf MCU Basis realisieren.

Dazu soll die Korrelation zum Einsatz kommen. Ich bin gerade dabei mich 
in das Thema einzulesen. Jedenfalls hab ich mittlerweile begriffen, dass 
man mittels Kreuzkorrelation die Laufzeit zwischen 2 Signalen bestimmen 
kann.

Gedanklich hab ich folgendes vor:
- Impuls erzeugen und aussenden (erstes Signal ist bekannt)
- Echo empfangen und Signal puffern
- jetzt das erste Signal im zweiten suchen
- Differenz ist dann die doppelte Laufzeit

Das ähnelt ja sehr dem Verfahren von Christoph.

Jetzt habe ich gelesen, dass durch die Korrelation auch bildgebende 
Algorithmen möglich sind, um beispielsweise eine Darstellung wie auf 
einem richtigen Sonar zu haben. Weiß einer wie das funktioniert? Das 
kann ich mir gerade nämlich nicht vorstellen.

[Christoph a. (chrichri)]

Florian F. schrieb:
> Jetzt habe ich gelesen, dass durch die Korrelation auch bildgebende
> Algorithmen möglich sind, um beispielsweise eine Darstellung wie auf
> einem richtigen Sonar zu haben. Weiß einer wie das funktioniert? Das
> kann ich mir gerade nämlich nicht vorstellen.

Das ist jetzt nur mein Verständnis der Sache ohne großartig Literatur zu 
lesen:
"Bildgebend" ist das Verfahren doch immer dann, wenn du eine Reihe von 
Daten, die du durch dein Verfahren bestimmt hast, interpretierts.

Z.B. Ein UBoot bewegt sich über den Meeresgrund, misst die 
Signallaufzeit zum Grund in festen Intervallen. Über die Laufzeit ist 
die Entfernung bestimmbar (z.B. Schallgeschwindigkeit unter Wasser). 
Somit kannst du ein 1-D Höhenprofil des Grundes (in Fahrtrichtung des 
UBootes) erstellen.

Was ich damit aber sagen will: IMHO ist das "bildgebende Verfahren" nur 
die Verwertung deiner Messdaten und hat erst einmal nichts mit dem 
Messverfahren zu tun.

Das ganze kann man natürlich noch sinnvoll verbessern. Z.B. wäre es bei 
dem UBoot angebracht den Grund in einer gewissen Breite abzutasten. 
Banale Lösung wäre: Einfach mehrer Sensoren parallel anordnen, dann hat 
man maherer 1-D Streifen und kann so ein 2-D Bild erstellen. Aber wieder 
bleibt das Messprinzip unberührt.
Hier wir aber der Nutzen der Korrelation deutlich: Die Sensoren können 
gleichzeitig Arbeiten, sie müssten aber nur jeder ein hinreichend 
unterschiedliches Signal aussenden, welches die anderen in iher 
Signalverarbeitung nicht stört.

[iPod-User (Gast)]

Christoph abc schrieb:
> Das kostet Zeit.

Das kostet das Scrollen durch deinen Code auch, insbesondere auf 
Mobilgeräten. Ist bei dir die Funktion zum Hochladen von Anhängen 
kaputt?

[Robert K. (Firma: Medizintechnik) (robident)]

Florian F. schrieb:
> - Echo empfangen und Signal puffern
>
> - jetzt das erste Signal im zweiten suchen

Wie willst du wissen, dass das, was Du empfängts, ein Echo ist und zu 
wissen, was du Puffern musst? Wenn Es das Echo ist, musst Du nicht 
suchen.

Wenn Du suchst, um es erst zu finden, dann musst Du alles aufnehmen und 
mehrfach suchen.

[Florian F. (ultrazauberer)]

Ich nehme alles auf und lass das Signal samplen. Wenn ich schon 
aufnehme, bevor ich einen Ping losschicke, habe ich bei der Aufnahme am 
Anfang ein ziemlich starkes Signal, was gleich als Muster herhalten kann 
und zur Synchronisation dient.

Danach mache ich eine diskrete, normierte Korrelation über das gesamte 
Signal und suche darin das Muster. Am Anfang werde ich somit ein Signal 
mit einem Korrelationskoeffizienten von 1 erhalten, was mein Muster ist. 
Alle Signale danach sind entweder Rauschen oder das Echo. Das Echo wird 
sicherlich auch irgendwas um die 0,5-0,9 als Korrelationskoeffizienten 
haben.

Wie sollte man es sonst am besten machen?

[Robert K. (Firma: Medizintechnik) (robident)]

Florian F. schrieb:
> Wie sollte man es sonst am besten machen?

Hört sich plausibel an, so. Ist die Frage, wie engmaschig du suchst, 
also von sample 1 ... sample 1024  und dann sample 2 ... saple 1025 oder 
gröber.

im Frequenzbereich (FFT) müsste es einfacher sein, weil dann die Phase 
rausfällt, wegen / sin/cos.

[A. S. (rava)]

ich weiß jetzt nicht genau, was du machen möchtest, aber ein matched 
filter sollte doch genau das tun, was du suchst, oder?

der Vorteil ist, er arbeitet im Zeitbereich. Du musst nur eine Faltung 
programmieren: du faltest das gemessene Signal mit dem Signal, das du 
finden möchtest und machst eine Grenzwert- und/oder Spitzendetektierung.


http://en.wikipedia.org/wiki/Matched_filter


ich habe keine Ahnung von Kreuzkorrelation, aber wenn die 
Aufgabenstellung tatsächlich eine "Suche" ist, ist mir nicht klar, wofür 
der Frequenzbereich nötig ist.
das Phasenargument kann ich nicht ganz nachvollziehen, denn die Faltung 
arbeitet ja kontinuierlich. Was also in den Daten ist, wird gefunden.

Oder nicht?

[Harald M. (mare_crisium)]

Florian,

ich habe so etwas einmal ausprobiert; es funktionierte auch ganz gut 
(siehe letzter Beitrag in diesem Thread: 
Beitrag "Ultraschall Merkmalsvektoren").

Als Signal wird bei diesem Verfahren eine Impulsfolge gesendet, die 
pseudo-zufällig ist, eine sogenanne MLS. Gleichzeitig mit dem Senden 
wird das empfangene Signal digitalisiert und gespeichert; das direkte 
Übersprechen vom Sender auf den Empfänger hat bei diesem speziellen 
Verfahren keinen störenden Einfluss, solange es den Empfänger nicht 
übersteuert. In meinem Versuchsaufbau hatte ich die beiden US-Kapseln 
ungefähr parallel ausgerichtet; das reichte schon aus, um diese 
Bedingung zu erfüllen.

Die MLS haben eine besondere mathematische Eigenschaft, die sich 
trickreich ausnutzen lässt, um das gesendete Signal im empfangenen 
wiederzufinden. Wie und warum das funktioniert, habe ich in dem 
.pdf-Dokument beschrieben, das an dem zitierten Beitrag angehängt ist. 
Die Mathematik dahinter ist etwas ungewohnt, lässt sich aber mit einem 
einfachen ATmega durchführen, der über genügend viel RAM verfügt.

Ciao,

mare_crisium

[Florian F. (ultrazauberer)]

Danke für das viele Feedback. Ich werde mir auf jeden Fall das PDF mal 
zu Gemüte führen. Den Ansatz, den ich damals beschrieben habe, hat auch 
funktioniert. Ich habe es nicht weiter verkompliziert.

Ich muss nochmal in meinen Unterlagen schauen, wo ich die Formel für die 
Kreuzkorrelation ist. Damit ging es zumindest eine einfache 
Echolotanwendung zu realisieren. Im Grunde habe ich die Formeln 1:1 als 
Programm umgesetzt und erstmal so etwas rumgespielt und getestet. 
Nachdem alles den Erwartungen entsprochen hat, habe ich es dann mit der 
FFT realisiert und damit waren dann wirkliche Leistungssprünge zu 
erreichen. :-)

[Udo S. (urschmitt)]

Falls ihr euch für Korrelationsfunktionen interessiert kann ich euch das 
Werk von Hänsler ans Herz legen
"Statistische Signale, Grundlagen und Anwendungen"

Ich durfte seine Vorlesung als Pflichtvorlesung hören und war froh 
meinen Schein geschafft zu haben :-)