Algorithme de Goertzel et Analyse du Signal

Algorithme de Goertzel et analyse du signal.

Dans un papier (MESA vs Goertzel), Dennis Meyers compare l'utilisation de l'algorithme de Goertzel avec la méthode de Ehlers pour l'analyse du signal : trouver le cycle dominant.  Dans l'article qui suit, je me base sur la publication de Meyers pour tester, comme il l'a fait, les 2 méthodes.  Cela m'a également été utile pour voir si les codes programmés sur PRT donnaient des résultats corrects.

J'ai donc créé, comme Meyers, un signal fictif :  1 sinusoïde de période 11 ajoutée à une sinusoïde de période 19 avec un décalage de phase, le tout avec une pincée de bruit.  Voici ce que cela donne, en fenêtre 1, le signal sans bruit (les 2 sinusoïdes mélangées) en fenêtre 2, le bruit (variant de -1 à +1), et en fenêtre 3, le signal avec le bruit.



Sur ce signal, j'applique l'algorithme de Goertzel et le programme de Ehlers pour essayer de retrouver les cycles dominants : soit les périodes 11 et 19.  Ci-dessous, le résultat avec l'algo de Goertzel.  Il y a bien une détection du cycle de période 11, mais elle n'est pas continuelle.  Pour le cycle de période 19, l'affichage est continu mais varie de 19 à 20.  Le programme retourne les deux fréquences qui ont la meilleure réponse.



Voici une vue du programme de Elhers "brut" (une seule fréquence est affichée), la détection des périodes 11 et 19 est moins évidente :



Et une autre vue avec le code bidouillé pour avoir 2 fréquences en sortie :



Le même, en élargissant la fenêtre d'analyse de 50 à 100 barres :



Conclusions : 
- Les codes ont l'air corrects, les cycles 11 et 19 apparaissent.
- La précision de la détection dépend pour les 2 programmes, de la longueur de le fenêtre de calcul.  Plus elle est longue, plus le calcul dure longtemps, et plus il est précis.  Meyers préconise 3 fois la longueur du cycle estimé, mais 5 ou 6 fois me paraît indispensable d'après les tests.
- Pour les 2 méthodes, il y a des périodes où le cycle 11 n'est pas détecté.  Par contre, le cycle 19 est toujours présent avec plus ou moins de précision.
- La structure des programmes est similaire, avec une subroutine, et donc c'est assez lourd à faire tourner.
- Reste à appliquer la détection sur les cours, préalablement "detrented", pour un prochain article.

Voici les différents codes pour Prorealtime, utilisés dans l'article :

/////////// signal sans bruit ////////////
b=barindex
sig=sin(360*b/11)+sin(360*b/19+90)
return sig,0

/////////////// bruit //////////////////
b=SQRT(barindex*(b+1))
b=abs((b-ROUND(b))*4)
c=b-1
return c,0

////////////////// sous prgm goertzel ////////////
a=barindex
sig=sin(360*a/11)+sin(360*a/19+90)
b=SQRT(barindex*(b+1))
b=abs((b-ROUND(b))*4)
c=b-1
sig=sig+c

pr=sig
alpha=2*cos(360*1/k)
q1=0
q2=0
for i=200 downto 0//fenêtre de 200 barres
    q3=pr[i]+alpha*q1-q2
    q2=q1
    q1=q3
next
amp=sqrt(square(q1)+square(q2)-alpha*q1*q2)
return amp

/////////////////// prgm goertzel ///////////////////
l=0
maxi=l
maxo=0
for p=5 to 25//détection de cycle de période de 5 à 25
    my = CALL "goertzel sous"[p]
    if my>maxi then
        per=p
        maxi=my
    endif
    if my>maxo and my<maxi then
        peri=p
        maxo=my
    endif
next
return max(per,peri),min(per,peri)

////////////// sous prgm ehlers ////////////////
a=barindex
sig=sin(360*a/11)+sin(360*a/19+90)
b=SQRT(barindex*(b+1))
b=abs((b-ROUND(b))*4)
c=b-1
sig=sig+c
cl=sig

cospart=0
sinpart=0
for n=0 to 99// fenetre de 100 barres
    cycper=(360*n)/peri
    cospart=cospart+cl[n]*cos(cycper)
    sinpart=sinpart+cl[n]*sin(cycper)
next
pwr=square(cospart)+square(sinpart)
return pwr

/////////////// prgm ehlers ///////////////
mix=0
maxpwr = CALL "sousprog dft"[8]
pero=8

for ki=9 to 30//détection de cycle de période de 9 à 30
    mys=call "sousprog dft"[ki]
    if mys>maxpwr then
        maxpwr=mys
        pero=ki
    endif
    if mys>mix and mys<maxpwr then
        mix=mys
        por=ki
    endif
next
return pero,por