Harmoniques

[passiion 0]

Voilà je voudrais savoir, plutôt simplement dans un premier temps, qu'est ce que c'est une décomposition en harmoniques en maths et comment on implique celà a la météo ? ( J'avais vu François Lott en parlé dans une de ses prestations ) .

Merci

[paix 0]

http://savoirsenmultimedia.ens.fr/expose.php?id=730

J'en parlais ici en fait :

http://www.forums.meteobelgium.be/index.ph...st&p=450198

Sur une fonction périodique, on peut modifier sa période pour avoir des harmoniques. Pour une onde zonale, avec x la longitude, tu peux écrire ainsi :

Z500 = a₀ + [ a₁*cos(1*ω*x) + b₁*sin(1*ω*x) ] + [ a₂*cos(2*ω*x) + b₂*sin(2*ω*x) ] + [ a₃*...

C'est une décomposition en série de Fourrier tout ce qu'il y a de plus classique. Les a_k et les b_k du nombre d'onde k. Le nombre d'onde est 1 pour le premier gros paquet, 2 pour le deuxième, 3 pour le troisième,... et a0 représente l'état de base. Ici ω est la pulsation :

ω = 2 * π * f avec f, la fréquence ( 1/T ). Si tu utilise par exemple une grille de longitude en degré :

ω = 2 * π * (1/360)

ω = 0.017453 et des brouettes

Les harmoniques, ce sont les multiples de ta pulsation de base. Visuellement (j'ai tiré cela d'Internet ^^ )

Ta fréquence fondamentale est la courbe violette, avec un maximum et un minimum sur ta période. Après tu as le signal bleu, qui admet deux maximums et deux minimums sur ta période, c'est la première harmonique. Puis la deuxième harmonique en vert. Le signal total est en orange.

Mathématiquement, il y a toujours moyen de revenir à l'exponentielle complexe évidemment ( exp(i*x) = cos(x) + i*sin(x) pour mémoire ), cela est plus compact Mais la décomposition en série de Fourier rend le truc plus évident. Bien sûr, tu peux toujours sortir l'amplitude alors :

Amplitude = Racine(a_k²+b_k²)

Phase = arctg(a_k/b_k) -> C'est là que tu as envie d'exploser ton ordi/ta caltoche en général vu que arctg est défini sur [ -pi/2 ; pi/2 ] XD

Si on se fait une petite manipulation sur les données du Z500 du 18 Décembre (récupéré du NCEP/NCAR). Je te passe les codes de R avec, si tu veux faire mumuse

> local({pkg <- select.list(sort(.packages(all.available = TRUE)),graphics=TRUE)
+ if(nchar(pkg)) library(pkg, character.only=TRUE)})
> path<-"C:/Users/aslan/Documents/ERS/hgt.2012.nc"
> hgt1 = open.ncdf(con=path)
> tape1=t(get.var.ncdf(hgt1,start=matrix(c(1,13,6,352)),count=matrix(c(144,1,1,1))))
> tape1<-t(tape1)

Tu ouvres ton fichier et tu balance les données du Z500. Tu prend les 144 points de longitude, la latitude 60°N (la 13ème en partant du Nord dans la grille du NCEP NCAR), le niveau 500 (6ème niveau, le NCEP/NCAR a les niveaux 1000 / 925 / 850 / 700 / 600 / 500 / 400 / 300 /... Le 500 c'est bien le 6ème ) et le jour du 18 Décembre, le 352ème jour. Tu plot pour voir si cela a une gueule qui ressemble à quelque chose ou pas :

>plot(x,tape1,main="Z500 60°N 18 Déc",xlab="Longitude",ylab="Z500, m")

Par rapport aux cartes, tu vois que tu es bon :

Tu as les bas Z500 qui s'enfoncent sur l'Atlantique, puis le pic de la poussé subtropicale vers Kara, puis la banane de bas Z500 du VP troposphérique sibérien, puis vers 200° sur le Pacifique, tu as une petite hausse qui est le PNA- (pas très vaillant il est vrai), puis encore une bonne tartine de Z500 avant d'avoir les bulles anticycloniques sur l'Est américain qui fait fantasmer le forum sur un très hypothétique changement de synop'

Tu peux sortir la régression linéaire alors de la série de Fourier :

> Banane<-lm(tape1~sin(omega*x)+cos(omega*x)+sin(omega*2*x)+cos(omega*2*x)+sin(omega*3*x)+cos(omega*3*x)+sin(omega*4*x)+cos(omega*4*x)+sin(omega*5*x)+cos(omega*5*x)+sin(omega*6*x)+cos(omega*6*x)+sin(omega*7*x)+cos(omega*7*x)+sin(omega*8*x)+cos(omega*8*x))
> summary(Banane)

Le premier bloc, c'est la première harmonique, celle en "omega*x". Là on va chercher 8 harmoniques, ou 8 nombre d'ondes zonal

Le résumé a cette tronche :

Call:
lm(formula = tape1 ~ sin(omega * x) + cos(omega * x) + sin(omega * 
    2 * x) + cos(omega * 2 * x) + sin(omega * 3 * x) + cos(omega * 
    3 * x) + sin(omega * 4 * x) + cos(omega * 4 * x) + sin(omega * 
    5 * x) + cos(omega * 5 * x) + sin(omega * 6 * x) + cos(omega * 
    6 * x) + sin(omega * 7 * x) + cos(omega * 7 * x) + sin(omega * 
    8 * x) + cos(omega * 8 * x))

Residuals:
    Min      1Q  Median      3Q     Max 
-32.734  -7.722   0.366   6.635  31.657 

Coefficients:
                   Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        5284.181      1.042 5072.202  < 2e-16 ***
sin(omega * x)       46.757      1.473   31.736  < 2e-16 ***
cos(omega * x)      195.005      1.473  132.358  < 2e-16 ***
sin(omega * 2 * x)   41.704      1.473   28.307  < 2e-16 ***
cos(omega * 2 * x)  -97.303      1.473  -66.044  < 2e-16 ***
sin(omega * 3 * x)    1.336      1.473    0.907 0.366148    
cos(omega * 3 * x)  -54.034      1.473  -36.675  < 2e-16 ***
sin(omega * 4 * x)   31.808      1.473   21.589  < 2e-16 ***
cos(omega * 4 * x)    9.001      1.473    6.109 1.14e-08 ***
sin(omega * 5 * x)  -14.448      1.473   -9.807  < 2e-16 ***
cos(omega * 5 * x)   -3.476      1.473   -2.360 0.019822 *  
sin(omega * 6 * x)   -8.491      1.473   -5.763 5.91e-08 ***
cos(omega * 6 * x)   -2.613      1.473   -1.774 0.078511 .  
sin(omega * 7 * x)   -5.080      1.473   -3.448 0.000767 ***
cos(omega * 7 * x)   18.896      1.473   12.825  < 2e-16 ***
sin(omega * 8 * x)   -1.830      1.473   -1.242 0.216383    
cos(omega * 8 * x)   -9.732      1.473   -6.605 9.82e-10 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

Residual standard error: 12.5 on 127 degrees of freedom
Multiple R-squared: 0.9951,     Adjusted R-squared: 0.9945 
F-statistic:  1619 on 16 and 127 DF,  p-value: < 2.2e-16

Cela charge surtout sur l'onde 1 et 2

Graphiquement :

En bleu, la courbe des données, en rouge le modèle avec 8 harmoniques. Cela se tient pas trop mal Dans la troposphère en géénral, les 5 premiers nombre d'ondes sont suffisant pour modéliser. En stratosphère, à 3 cela suffit (plus grande stabilité donc plus de filtrage, seules les ondes 1 et 2, et il peut rester un petit résidu en 3).

Par exemple, l'onde 1 a pour amplitude :

Amplitude(onde 1) = racine ( 195.055² + 46.757²) = 200.53m

Et pour phase :

Phase(onde 1) = 90 - arctg(195.055/46.757) = 13.5°

Tu le vérifie sur le graphique, le maximum de ta fonction pour la première harmonique est à 13.5°, et de crête à crête on a 400 mètres environ.

Et ainsi de suite, pour l'onde 2 tu as une amplitude d'un peu plus de 100m, environ 50m pour l'onde 3,...

Ce qui est intéressant, comme tu me le disais, les ondes de grandes amplitudes sont quasiment stationnaire. Si on regarde 5 jours plutôt, le 13 Décembre :

Coefficients:
                    Estimate Std. Error  t value Pr(>|t|)    
(Intercept)        5267.9514     0.9921 5309.972  < 2e-16 ***
sin(omega * x)       46.0504     1.4030   32.822  < 2e-16 ***
cos(omega * x)      108.5641     1.4030   77.379  < 2e-16 ***
sin(omega * 2 * x)   84.3468     1.4030   60.118  < 2e-16 ***
cos(omega * 2 * x)   43.7171     1.4030   31.159  < 2e-16 ***
sin(omega * 3 * x)  -58.0551     1.4030  -41.379  < 2e-16 ***
cos(omega * 3 * x)   -0.1152     1.4030   -0.082  0.93470    
sin(omega * 4 * x)   10.8386     1.4030    7.725 2.93e-12 ***
cos(omega * 4 * x)  -94.2523     1.4030  -67.178  < 2e-16 ***
sin(omega * 5 * x)  -45.5645     1.4030  -32.476  < 2e-16 ***
cos(omega * 5 * x)   -3.7924     1.4030   -2.703  0.00781 ** 
sin(omega * 6 * x)  -37.7764     1.4030  -26.925  < 2e-16 ***
cos(omega * 6 * x)    2.3356     1.4030    1.665  0.09844 .  
sin(omega * 7 * x)  -17.2723     1.4030  -12.311  < 2e-16 ***
cos(omega * 7 * x)   15.6027     1.4030   11.121  < 2e-16 ***
sin(omega * 8 * x)  -28.5112     1.4030  -20.321  < 2e-16 ***
cos(omega * 8 * x)    5.6458     1.4030    4.024 9.77e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Les ondes de grandes amplitudes n'ont pratiquement pas bougé malgré que le Z500 a une tronche apparemment différent. C'est donc un moyen utile pour faire de la prévision. Par contre, les données des ondulations ne font pas partie des sorties des modèles ^^ cela eprmet pourtant de sortir la variabilité basse fréquence et ne pas s'attacher à un Z500 qui n'est que le produit final de tout ceci et n'a pas grande utilité à moyen terme. Tu remarque un phasage de l'onde 2 par exemple avec le PNA- (qui tend à s'affaiblir cependant) et qui nous place un creux planétaire juste à l'Ouest de chez nous, sur l'Atlantique. Et là, c'est stable à l'échelle de la semaine.

Modifié 19 décembre 2012 par paix

[passiion 0]

Ouai je vois un peu le truc, on décompose une onde en signaux plus précis en gros

Merci

[paix 0]

Ouai je vois un peu le truc, on décompose une onde en signaux plus précis en gros

Merci

Oui c'est une manière de décomposer un gros machin ingérable qu'est le Z500 en différents signaux. Si tu arrive à choper le signal basse fréquence, c'est gagné