[Alg/C++] Lokale extremen vinden in noisy data

[Zoijar]

Ik zoek een makkelijke maar vooral robuste manier om de lokale maxima/minima te vinden van een gesampelde functie. In dit geval is het een sinus-achtige golf.

Input: array van floats met integer index (std::vector<float>)
Output: array van integer indices voor alle lokale maxima en minima

Wat ik nu doe:

C++:

void find_peaks(const std::vector<double>& ys, std::vector<int>& mins, std::vector<int>& maxs) {
    int last_dir = 0;

    for (int i=1; i<ys.size(); ++i) {
        int dir;
        
        if (ys[i] > ys[i-1]) {dir = +1;}
        else if (ys[i] < ys[i-1]) {dir = -1;}
        else {dir = last_dir;}

        if (dir != 0 && dir != last_dir) {
            if (dir > 0) {
                mins.push_back(i);
            }
            if (dir < 0) {
                maxs.push_back(i);
            }
        }
        last_dir = dir;
    }
    
}

Het probleem daarmee is dat er soms iets van ruis in de data zit, en er dan hele kleine lokale extremen worden gevonden. Bijvoorbeeld hier:

code:

-0.969114
-0.976812
-0.987645
-0.993775
-0.995343
-1
-0.999953 // dit stukje gaat mis
-1
-0.996959
-0.989308
-0.983179
-0.969398
-0.950962

Maar soms is dat over iets meer samples verspreid: het is niet altijd alleen de volgende die fout gaat. Een ander probleem is als het maximum gelijk is, dan zou hij het midden moeten proberen te vinden:

code:

0.992729
0.997196
0.997196
1 // dubbel max, gaat mis
1
1
0.998812
0.995675
0.989688

Nu kan ik wel de uiteindelijke lijst weer doorlopen om dit soort dingen op te sporen en dan samen te voegen, maar het wordt allemaal een beetje ingewikkeld en ad-hoc. Dat moet toch makkelijker kunnen? Hebben jullie ideeen?

[pkuppens]

Kun je nog meer vertellen over je probleem?
Is het altijd een sinus-functie en wil je dit gebruiken om de periode te bepalen?
Of is het een sinus met harmonischen?
Orde grootte ruis of signaal ruis verhouding?
Je sinus lijkt een heel grote sampling rate te hebben, dan verwacht ik ook dat je bv. maar per 10 punten een lokaal extreem wilt hebben?

Een beetje afhankelijk van de antwoorden op die vragen zou je bv. eerst het signaal door een filter kunnen gooien, bv (gewogen) gemiddelde van +/- 5 punten. Het filter moet de curve volgen, en niet de ruis. En dan jouw algoritme erop loslaten.

Signaal-theoretisch gezien zou je de Fourier transformatie kunnen nemen, en de hogere orden (die de ruis bevatten) weggooien. Evt. terugtransformeren, maar de kans is er dat je antwoorden op het onderliggend probleem ook uit de fase informatie kunt halen.

[Zoijar]

Uiteindelijk wil ik het gebruiken om het fase verschil tussen twee slingers te meten. Voorbeeld data is dit:



Eerste idee was om het verschil in periode tussen alle pieken en dalen te meten, en dat nemen als gemiddeld fase verschil tussen de twee signalen. Dit is 60s signaal gesampled op 60Hz. Fase verschil zal zo ergens in de buurt van de 80ms liggen.

Een andere optie is om op allebei een sinus functie te fitten (bv (1-c*x)*sin(a*x+b) en dan het verschil tussen de 'a' parameters te berekenen... maar dat ging niet helemaal lekker omdat het niet exact sinus golven zijn.

(maar ik wil eigenlijk geen signaal-specifieke parameters hoeven te gebruiken; het moet gewoon altijd werken op soortgelijke signalen, ookal is de sampling anders, of de shift, of de amplitude etc)

[ Voor 13% gewijzigd door Zoijar op 05-07-2008 14:49 ]

[Soultaker]

Kun je niet handig gebruikmaken van het gegeven dat het een slingerbeweging om de oorsprong is? Aangezien je alternerend toppen en dalen wil vinden, en de toppen boven de x-as liggen en de dalen eronder, kun je zoiets doen:

code:

begin := 0
dir := -1
for i := 0 to aantal_samples
begin
  if (dir == -1 and data[i] > 0) or  (dir == +1 and data[i] < 0)
  then
  begin
     eind := i
     // vind top in [begin,eind)
     begin := eind
     dir := -dir
  end
end

Op deze manier splits je de data op in halve golven (is daar een term voor?); in elk deel wil je dan precies één top/dal aanwijzen (ook als er misschien toevallig meerdere zijn). Je kunt dat doen door b.v. het middelste element met de grootste absolute waarde te nemen o.i.d.

[pkuppens]

Blijft je faseverschil ook constant over het signaal? Bij een klein verschil in frequentie van de signalen verloopt wat je denkt dat de fase is.

Ik ben niet helemaal fris meer met m'n Laplace transformaties, maar dit signaal laat zich prima beschrijven als het reele deel van e^(s*t+s0), met daarin frequentie en fase, amplitude en demping opgesloten in complexe getallen s en s0.
Of de Z transformatie als discrete variant. Als je al een idee hebt over bemonsteringsfrequentie, en signaal frequentie, komt daar met z = e^(2*pi*i*t*f/b) o.i.d een complex getal uit. Het imaginaire deel vertel wat over de fase.
De FFT geeft ook een prima benadering, als je de eerste piek in het spectrum gevonden hebt.
Het imaginaire gedeelte geeft dan informatie over de fase. Al dan niet met een factor 2*pi.

[Zoijar]

Ok, bedankt. Ik heb over de sugegsties nagedacht. Na het fourier verhaal ben ik toch even gaan kijken naar een FFT. Die ziet er zo uit voor beide functies:



Beide FFT pieken liggen op index 33, ofwel F=0.469Hz. De waardes uit de FFT data hierbij zijn:

p1 = 0.2734 + 0.1792i
p2 = 0.3000 + 0.1231i

Het verschil in hoeken da= (angle(p1)-angle(p2)) is 0.1907 graden. Als ik nu doe (da/360)*(1/F)*60, waar 60 de sampling rate in hz is, da het hoekverschil, F de piek frequentie, dan kom ik uit op 0.0678, ofwel 67.8ms fase verschil. Dat klopt aardig. Als ik deze methode gebruik op syntetische sinus golven krijg ik ook de juiste fase shift eruit, met een kleine fout (is dat FFT zero-padding oid?)

Klopt deze methode? Mag/kan dit?

Follow-up: Nu heb ik wel leuk een getal, maar hoe groot is mijn fout? Kan ik dat op de een of andere manier bepalen? Ik zou graag een 95% confidence interval oid krijgen. Niet 67ms zonder enige maat of dat [66,68]ms is of [47,87]ms... Ideeen?

[Verwijderd]

Waarom gebruik je geen kalman filter? In de regeltechniek worden die dingen heel vaak gebruikt om (witte) ruis weg te filteren.

[Zoijar]

Verwijderd schreef op zondag 06 juli 2008 @ 12:58:
Waarom gebruik je geen kalman filter? In de regeltechniek worden die dingen heel vaak gebruikt om (witte) ruis weg te filteren.

Volgens mij werkt de FFT wel goed; dan hoef ik helemaal geen pieken en dalen te zoeken. Er komt dit uit voor bovenstaande signalen gesampled op 60Hz over ~60sec. Volgens mij klopt dat. (latency is fase verschil)

code:

Read 3498 frames.
   1: Abs = 0.324258 f = 0.46875Hz T = 2.13333 Arg = 1.18135
   2: Abs = 0.326939 f = 0.46875Hz T = 2.13333 Arg = 0.990604
   Phase shift = 0.190746
   Latency = 64.7641ms
--------- blocks of 256 frames
Frames 0 - 255 :
   1: Abs = 0.496717 f = 0.46875Hz T = 2.13333 Arg = 1.0707
   2: Abs = 0.495867 f = 0.46875Hz T = 2.13333 Arg = 0.879942
   Phase shift = 0.190754
   Latency = 64.7669ms
Frames 256 - 511 :
   1: Abs = 0.456923 f = 0.46875Hz T = 2.13333 Arg = 1.10897
   2: Abs = 0.458417 f = 0.46875Hz T = 2.13333 Arg = 0.918807
   Phase shift = 0.190163
   Latency = 64.5663ms
Frames 512 - 767 :
   1: Abs = 0.421701 f = 0.46875Hz T = 2.13333 Arg = 1.15056
   2: Abs = 0.423841 f = 0.46875Hz T = 2.13333 Arg = 0.959573
   Phase shift = 0.190991
   Latency = 64.8474ms
Frames 768 - 1023 :
   1: Abs = 0.390809 f = 0.46875Hz T = 2.13333 Arg = 1.17787
   2: Abs = 0.393564 f = 0.46875Hz T = 2.13333 Arg = 0.986997
   Phase shift = 0.19087
   Latency = 64.8062ms
Frames 1024 - 1279 :
   1: Abs = 0.363178 f = 0.46875Hz T = 2.13333 Arg = 1.19954
   2: Abs = 0.366359 f = 0.46875Hz T = 2.13333 Arg = 1.00813
   Phase shift = 0.191403
   Latency = 64.987ms
Frames 1280 - 1535 :
   1: Abs = 0.338972 f = 0.46875Hz T = 2.13333 Arg = 1.21418
   2: Abs = 0.341903 f = 0.46875Hz T = 2.13333 Arg = 1.02235
   Phase shift = 0.191833
   Latency = 65.1333ms
Frames 1536 - 1791 :
   1: Abs = 0.316494 f = 0.46875Hz T = 2.13333 Arg = 1.22241
   2: Abs = 0.319711 f = 0.46875Hz T = 2.13333 Arg = 1.03127
   Phase shift = 0.191144
   Latency = 64.8991ms
Frames 1792 - 2047 :
   1: Abs = 0.296411 f = 0.46875Hz T = 2.13333 Arg = 1.22977
   2: Abs = 0.299473 f = 0.46875Hz T = 2.13333 Arg = 1.03797
   Phase shift = 0.191792
   Latency = 65.1194ms
Frames 2048 - 2303 :
   1: Abs = 0.277622 f = 0.46875Hz T = 2.13333 Arg = 1.23382
   2: Abs = 0.280949 f = 0.46875Hz T = 2.13333 Arg = 1.04077
   Phase shift = 0.193055
   Latency = 65.5479ms
Frames 2304 - 2559 :
   1: Abs = 0.260369 f = 0.46875Hz T = 2.13333 Arg = 1.23657
   2: Abs = 0.264004 f = 0.46875Hz T = 2.13333 Arg = 1.03994
   Phase shift = 0.196629
   Latency = 66.7616ms
Frames 2560 - 2815 :
   1: Abs = 0.245394 f = 0.46875Hz T = 2.13333 Arg = 1.22891
   2: Abs = 0.248671 f = 0.46875Hz T = 2.13333 Arg = 1.03365
   Phase shift = 0.195259
   Latency = 66.2965ms
Frames 2816 - 3071 :
   1: Abs = 0.231116 f = 0.46875Hz T = 2.13333 Arg = 1.21722
   2: Abs = 0.23435 f = 0.46875Hz T = 2.13333 Arg = 1.02281
   Phase shift = 0.194412
   Latency = 66.0089ms
Frames 3072 - 3327 :
   1: Abs = 0.217829 f = 0.46875Hz T = 2.13333 Arg = 1.20315
   2: Abs = 0.22123 f = 0.46875Hz T = 2.13333 Arg = 1.00989
   Phase shift = 0.193259
   Latency = 65.6172ms

[ Voor 25% gewijzigd door Zoijar op 06-07-2008 14:22 ]

[Ruben314]

> maar hoe groot is mijn fout?

Je kunt het verschil bepalen tussen je (met fft) gefilterde signalen en je originele signaal. Bijvoorbeeld met de kleinste kwadraten methode oid.
Je kunt dan een getal geven voor hoe goed je meetsignaal past op een bepaalde slingerbeweging.

Vervolgens verschuif je je gefitte (fft) signaal iets in fase. Dat vergelijk je dan weer met je oorspronkelijke meetsignaal. Daar komt ook weer een getal uit. Dat getal kun je vergelijken met het resultaat van je eerste vergelijking.

Doe dat een aantal keer en je kunt uitzetten hoe goed je gemeten signaal past op een 60Hz signaal met verschillende fasen. Dat zou een redelijke schatter kunnen opleveren voor je nauwkeurigheid.

Sorry als het een beetje een wazig verhaal geworden is. Ben ook geen statisticus. Hopelijk heb je d'r wat aan.

[pkuppens]

Zoijar schreef op zaterdag 05 juli 2008 @ 20:03:
Ok, bedankt. Ik heb over de sugegsties nagedacht. Na het fourier verhaal ben ik toch even gaan kijken naar een FFT. Die ziet er zo uit voor beide functies:

[afbeelding]

Beide FFT pieken liggen op index 33, ofwel F=0.469Hz. De waardes uit de FFT data hierbij zijn:

p1 = 0.2734 + 0.1792i
p2 = 0.3000 + 0.1231i

Het verschil in hoeken da= (angle(p1)-angle(p2)) is 0.1907 graden. Als ik nu doe (da/360)*(1/F)*60, waar 60 de sampling rate in hz is, da het hoekverschil, F de piek frequentie, dan kom ik uit op 0.0678, ofwel 67.8ms fase verschil. Dat klopt aardig. Als ik deze methode gebruik op syntetische sinus golven krijg ik ook de juiste fase shift eruit, met een kleine fout (is dat FFT zero-padding oid?)

Klopt deze methode? Mag/kan dit?

Follow-up: Nu heb ik wel leuk een getal, maar hoe groot is mijn fout? Kan ik dat op de een of andere manier bepalen? Ik zou graag een 95% confidence interval oid krijgen. Niet 67ms zonder enige maat of dat [66,68]ms is of [47,87]ms... Ideeen?

Het probleem wat jij denkt dat zero-padding zou kunnen zijn is misschien iets anders, namelijk dat je complete bemonstering niet een geheel veelvoud is van de periode van het signaal. De FFT is dan ook in theorie niet exact. (Het is wel hetzelfde als bv 500 punten wel de complete periode is, maar je een FFT op 512 punten doet met zero padding, terwijl een DFT op 500 punten had gewerkt.)

De frequentie ligt dan niet precies op een top, maar zou tussen twee toppen in vallen.
Leuk is wel dat de FFT het spectrum in de buurt neerzet.
De piek zit bij 33, maar kijk ook eens naar de buurman. Wat ik ook voor een eindopdracht Signaalverwerking TU/e Wiskunde, 2e jaar, heb gedaan is om om de frequentie beter te bepalen is gewogen middelen met de buurman. Je piek zou dan bv. 33.3 kunnen zijn. Dat gaf al een heel aardig resultaat, maar was zeker geen eindoplossing omdat het gevoelig was voor harmonischen en voor faseverschuiving.

Het verschil tussen de sinus bij 0.469 Hz en 0.500 Hz die bij 33.3 zou kunnen horen is een laagfrequent signaal. Het heeft wel een spectrum met wat lagere waarden, maar je ziet het wel aan die extra piekjes.

Ik heb destijds niet gekeken hoe de fases van die twee pieken zich verhielden, maar je huidige strategie lijkt me een heel goede, hooguit zou je die middeling ook op fase eens kunnen bekijken.

Als je een schatting hebt van de ruis, bekijk dan eens het spectrum van puur ruis (evt. gesimuleerd).
Fase-informatie zegt niks in ruis, je kunt wel kijken hoe nauwkeurig je je frequentie kunt krijgen.

Qua follow-up denk ik dat je echt een serieus probleem aanboort. De invloed van ruis, harmonischen, of een faseverschuiving is hiervoor al te groot.

[MSalters]

Minima en maxima vinden schiet inderdaad niet op. Veel stabiler is het vinden van nuldoorgangen; daarvoor kun je meer punten gebruiken. Juist omdat je daar snel van -0.5 naar +0.5 gaat maakt een fout van 0.01 weinig uit - je ruis is waarschijnlijk voornamelijk in de y richting?

Een FFT is nauwkeuriger, maar vooral als je de demping eerst compenseert. Anders is er een alternatief: demping is exponentieel; je sinusfunctie is in het complexe domein ook een exponentiële functie. Je moet dus een verzameling complexe exponenten schatten.

[roy-t]

MSalters schreef op maandag 07 juli 2008 @ 20:39:
Minima en maxima vinden schiet inderdaad niet op.

Door eerst een afgeleide functie te maken (f ' ) zou je minima en maxima makkelijk kunnen vinden door kruisingen met de 0 as. Hoewel dit het ruis probleem niet oplost natuurlijk.

[Zoijar]

Bedankt voor de suggesties, ik zal er nog eens wat aan prutsen. In ieder geval werkt de basis methode nu heel goed. Je ziet ook dat over het hele signaal de amplitude van de sterkste harmonische slechts 0.64 is door de demping, terwijl dat dus 1 zou moeten zijn. Maar als ik het in kleinere stukejs hak, dan zijn die twee amplitudes (abs(c) waar c de complexe fourier piek is) wel steeds nagenoeg gelijk. In beiden gevallen is het fase verschil vrijwel gelijk, dus het FFT maximmum geeft een hele redelijke schatting. Ik heb het vandaag ook op echte data geprobeerd door er 50ms fase-shift bij te simuleren en dat kwam precies terug in de metingen. Misschien kan het iets beter, of iets nauwkeuriger, maar het werkt voor de praktijk nu vrij goed. In ieder geval was een FFT een heel goed idee: veel beter dan zelf pieken of zero-crossings zoeken.

[MSalters]

roy-t schreef op maandag 07 juli 2008 @ 20:52:
[...]
Door eerst een afgeleide functie te maken (f ' ) zou je minima en maxima makkelijk kunnen vinden door kruisingen met de 0 as. Hoewel dit het ruis probleem niet oplost natuurlijk.

Eh, nee. Je kunt geen afgeleide functie maken van een meetreeks. Je kunt wel de verschillen berekenen. Die verschillen hebben alleen een dubbele ruisterm: x1+e1-x0-e0 = dx/dt + e1 + -e0.
Het voordeel van mijn methode is dat je de fout e(t)+e(t+1) deelt door dx/dt.

[pkuppens]

MSalters schreef op maandag 07 juli 2008 @ 20:39:
Minima en maxima vinden schiet inderdaad niet op. Veel stabiler is het vinden van nuldoorgangen; daarvoor kun je meer punten gebruiken. Juist omdat je daar snel van -0.5 naar +0.5 gaat maakt een fout van 0.01 weinig uit - je ruis is waarschijnlijk voornamelijk in de y richting?

Een FFT is nauwkeuriger, maar vooral als je de demping eerst compenseert. Anders is er een alternatief: demping is exponentieel; je sinusfunctie is in het complexe domein ook een exponentiële functie. Je moet dus een verzameling complexe exponenten schatten.

Dat is ook leuk! (Wel wat ingewikkelder?)
In het complexe vlak het optimum van de Laplace getransformeerde zoeken (Z transformatie).
Heb je zowel frequentie als demping te pakken, waar het optimum zit, en de waarde in het optimum geeft amplitude en fase.

[Verwijderd]

Op basis van het idee dat MSalters al heeft aangedragen: als je signaal een sinus-achtige golf is met een frequentie en -variatie die aanzienlijk geringer is dan je sample frequentie, dan zou je de volgende strategie kunnen overwegen:

1) bepaal uit de metingen eerst een aantal nuldoorgangen, bijvoorbeeld 50, en zet de tijdsverschillen tussen twee opvolgende doorgangen uit op een tijdschaal;
2) bepaal op basis hiervan een indicatie van de min of meer constant geachte tussenaankomsttijd;
3) filter dan de nuldoorgangs-spikes eruit, die eenvoudig te vinden zijn omdat de gemiddelde tussenaankomsttijd tussen twee nuldoorgangen afwijkt, bijvoorbeeld omdat de tijd tussen twee doorgangen aanzienlijk kleiner is dan verwacht. Hierna kun je de gemiddelde frequentie op het gekozen interval beter schatten, en daarmee een exactere schatting van de 'echte' nuldoorgangen;
4) maxima en minima liggen waarschijnlijk ongeveer in het midden van twee 'echte' nuldoorgangen. Ga daar zoeken in je originele data. Skip wel de gevonden ruis: indien de variaties in signaalsterkte binnen een lokale bandbreedte liggen, clip (ten gevolge van ruis) sterk afwijkende extremen dan aan die bandbreedte.

Merk op dat deze simpele methode voor wat betreft het bepalen van de frequentie behoorlijk robuust is bij niet al te grote variaties in je signaalsterkte, maar (iig initieel) wel een samplefrequentie vereist die flink hoger is dan de signaalfrequentie.

Indien er frequentie variatie in je signaal is, dan zou je kunnen overwegen om dit verhaal opnieuw uit te voeren na een aantal nuldoorgangen - wel met enige overlap natuurlijk.

[Zoijar]

Ik vind die Z-transform wel interessant; ga ik naar kijken. Nog nooit gebruikt, dus dan leer ik weer eens wat nieuws