Korrelationskoeffizient zweier Arrays

@ 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.

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)?

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 ;)

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!

1

#define baseN 4

2

#define N 256

3

4

#define PEARSON8BIT

5

6

/*              (sumXY - N * meanX * meanY)^2

7

 * r^2 = ------------------------------------------

8

 *       (sumXX - N * meanX^2)(sumYY - N * meanY^2)

9

 */

10

uint8_t pearson(uint8_t* dataX, uint8_t* dataY)

11

{

12

#ifdef PEARSON8BIT

13

    // pre-calculated values for reference

14

    const int16_t meanX = 1226;

15

    const int24_t denomX = 64283;

16

    int16_t sumY = 0;

17

    int24_t sumXY = 0, sumYY = 0;

18

    int16_t meanY;

19

#else

20

    // meanX = sumX / N * sqrt(N)

21

    const int8_t meanX = 1;

22

    // denomX = sumXX - N * meanX ^ 2

23

    const int16_t denomX = 1;

24

    int8_t sumY = 0;

25

    int16_t sumXY = 0, sumYY = 0;

26

    int8_t meanY;

27

#endif

28

29

    uint16_t i;

30

31

    // add up sums

32

    for (i = 0; i < N /* SPECPOINTS*/; i++)

33

    {

34

        // sum of all y-values

35

        sumY += dataY[i];

36

        // sum of squares of y-values

37

        sumYY += (uint16_t) dataY[i] * (uint16_t) dataY[i];

38

        // sum of product of x- and y-values

39

        sumXY += (uint16_t) dataX[i] * (uint16_t) dataY[i];

40

    }

41

42

    // calculate meanY-value as sqrt(N) * meanY

43

    meanY = sumY >> baseN;

44

    // numerator equals sumXY - N * meanX * meanY

45

    sumXY -= (uint24_t) meanX * (uint24_t) meanY;

46

    // the y-denominator equals sumYY - meanY ^ 2

47

    sumYY -= (uint24_t) meanY * (uint24_t) meanY; 

48

49

    //TODO: check for overflow

50

51

    // r ^ 2 should be mapped to 0 .. 255 so divide denominators by 16 equals

52

    // overall multiplication by 256 but increases accuracy

53

    sumYY >>= 4;

54

    // increase numerator to increase accuracy

55

    sumXY <<= 0;

56

57

    // check for divide by zero

58

    if (sumYY == 0)

59

    {

60

        return 255;

61

    }

62

    else

63

    {

64

        //          num      num

65

        // r ^ 2 = ------ x ------

66

        //         denomX   denomY

67

68

        sumYY = sumXY / sumYY;

69

        sumXY = sumXY / denomX;

70

        sumXY *= sumYY;

71

        // scale back down and return

72

        return (uint8_t) (sumXY >> 0);

73

    }

74

}