Korrelationskoeffizient zweier Arrays

[Enrico K. (ekoeck)]

Hallöchen!

Ich möchte zwei Arrays (256 Werte, 8 Bit) auf Ähnlichkeit überprüfen. 
Mit genügend Rechenpower würde ich einfach stumpf eine 
Korrelationsfunktion drüber jagen, habe jedoch nur einen 8 Bit uC mit 
begrenzten Ressourcen zur Verfügung und recht wenig Zeit für das ganze 
(jippie!). Daher bin ich auf der Suche nach einem effektiven 
Algorithmus, meine DSV-Büchlein haben dazu aber allesamt keine großartig 
hilfreichen Beiträge zu leisten.
Meine Ideen sind bis dato so weit fortgeschritten, dass wahrscheinlich 
eine Art Clustering sinnvoll wäre, bin aber gerade recht unmotiviert das 
ganze groß mathematisch durchzukauen und würde daher gerne schon 
erworbenes Wissen recyclen. Daher meine Frage: hat jemand eine 
Literaturempfehlung für mich oder kennt entsprechende Papers o.ä.? 
Freund G hat mir nicht viel verwertbares liefern können, da die meisten 
Algorithmen (z.B. Bildverarbeitung) schlicht zu komplex sind und in 
meinem Fall einfach überdimensioniert wären.

Wäre super, falls jemand einen hilfreichen Tip hat, schönes sonniges 
Wochenende wünsche ich des Weiteren ;)

[Falk B. (falk)]

@ Enrico Koeck (ekoeck)

>Ich möchte zwei Arrays (256 Werte, 8 Bit) auf Ähnlichkeit überprüfen.
>Mit genügend Rechenpower würde ich einfach stumpf eine
>Korrelationsfunktion drüber jagen, habe jedoch nur einen 8 Bit uC mit
>begrenzten Ressourcen zur Verfügung und recht wenig Zeit für das ganze
>(jippie!). Daher bin ich auf der Suche nach einem effektiven
>Algorithmus, meine DSV-Büchlein haben dazu aber allesamt keine großartig
>hilfreichen Beiträge zu leisten.

Für Korellation gibt es wie für DFT eine schnelle Variante (FFT).
Wie die jetzt aber genau heißt und funktioniert, weiß ich im Moment 
nicht.

[Andreas S. (andreas) (Admin)]

Enrico Koeck schrieb:
> Ich möchte zwei Arrays (256 Werte, 8 Bit) auf Ähnlichkeit überprüfen

Was für Signale sind das, und wie definierst du Ähnlichkeit (z.B., ist 
Verschiebung, Skalierung erlaubt)?

[Enrico K. (ekoeck)]

Falk Brunner schrieb:

> Für Korellation gibt es wie für DFT eine schnelle Variante (FFT).
> Wie die jetzt aber genau heißt und funktioniert, weiß ich im Moment
> nicht.

Das habe ich auch schon gefunden, erscheint mir aber zu aufwändig 
besonders bei der geringen Datenmenge.

Andreas Schwarz schrieb:
> Enrico Koeck schrieb:
>> Ich möchte zwei Arrays (256 Werte, 8 Bit) auf Ähnlichkeit überprüfen
>
> Was für Signale sind das, und wie definierst du Ähnlichkeit (z.B., ist
> Verschiebung, Skalierung erlaubt)?

Es handelt sich um ein von einem RF-Detektor aufgenommenes Spektrum. 
Durch verschiedene Spektren lassen sich verschiedene Sender 
identifizieren. Skalierung in der Höhe ist selbstverständlich gegeben 
(durch Umwelteinflüsse, Empfindlichkeit, Sendestärke, ...). In 
Y-Richtung (Frequenzachse) ist der Bereich aber fix, also bei jedem 
Sender eindeutig zuzuordnen. Ich habe auch schon überlegt mir markante 
Punkte herauszupicken und die zu speichern, das scheitert aber daran, 
dass zeitweilig einige Träger abgeschaltet sein können, es sich aber 
trotzdem um den selben Sender handelt. Ähnlichkeit definiere ich daher 
so, wie sie ein Mensch sehen würde so quasi: "Ja die Huckel und Täler 
sind ähnlich, da fehlt zwar etwas, es sieht aber schon recht passend 
aus." Ist halt recht schwierig mathematisch auszudrücken ;)

[Enrico K. (ekoeck)]

Falls es wen interessiert, ich habe mir einen Pearson-Algorithmus auf 
Integer umgebaut und etwas optimiert. Aktuell benötigt der Algorithmus 
bei 256 8-Bit Werten 6.5ms (PIC18F252, 20Mhz, MPLAB X IDE v2.00, XC8 
v1.33 free), was für mein Problem schon einigermaßen i.O. ist. 
Vielleicht hat ja jemand von den erfahreneren Programmierern noch einen 
Tip für mich, wo da noch etwas Zeit herauszuholen ist. Mit Reduzierung 
auf 4-Bit Datenbreite gewinnt man etwa 20%, die Genauigkeit sinkt aber 
auch erheblich.

TIA!

#define baseN 4
#define N 256

#define PEARSON8BIT

/*              (sumXY - N * meanX * meanY)^2
 * r^2 = ------------------------------------------
 *       (sumXX - N * meanX^2)(sumYY - N * meanY^2)
 */
uint8_t pearson(uint8_t* dataX, uint8_t* dataY)
{
#ifdef PEARSON8BIT
    // pre-calculated values for reference
    const int16_t meanX = 1226;
    const int24_t denomX = 64283;
    int16_t sumY = 0;
    int24_t sumXY = 0, sumYY = 0;
    int16_t meanY;
#else
    // meanX = sumX / N * sqrt(N)
    const int8_t meanX = 1;
    // denomX = sumXX - N * meanX ^ 2
    const int16_t denomX = 1;
    int8_t sumY = 0;
    int16_t sumXY = 0, sumYY = 0;
    int8_t meanY;
#endif

    uint16_t i;

    // add up sums
    for (i = 0; i < N /* SPECPOINTS*/; i++)
    {
        // sum of all y-values
        sumY += dataY[i];
        // sum of squares of y-values
        sumYY += (uint16_t) dataY[i] * (uint16_t) dataY[i];
        // sum of product of x- and y-values
        sumXY += (uint16_t) dataX[i] * (uint16_t) dataY[i];
    }

    // calculate meanY-value as sqrt(N) * meanY
    meanY = sumY >> baseN;
    // numerator equals sumXY - N * meanX * meanY
    sumXY -= (uint24_t) meanX * (uint24_t) meanY;
    // the y-denominator equals sumYY - meanY ^ 2
    sumYY -= (uint24_t) meanY * (uint24_t) meanY; 

    //TODO: check for overflow

    // r ^ 2 should be mapped to 0 .. 255 so divide denominators by 16 equals
    // overall multiplication by 256 but increases accuracy
    sumYY >>= 4;
    // increase numerator to increase accuracy
    sumXY <<= 0;

    // check for divide by zero
    if (sumYY == 0)
    {
        return 255;
    }
    else
    {
        //          num      num
        // r ^ 2 = ------ x ------
        //         denomX   denomY

        sumYY = sumXY / sumYY;
        sumXY = sumXY / denomX;
        sumXY *= sumYY;
        // scale back down and return
        return (uint8_t) (sumXY >> 0);
    }
}