Discussion :

[Julien N]

Bonjour,

N'ayant pour une fois pas trouvé ce que je cherche dans les anciens posts, je me permet d'en créer un nouveau.

Voici mon problème.

J'ai en ma possession un script VBS qui récupère, entre autre, un jeu de points depuis un fichier texte correspondant à des positions en fonction du temps. Le signal discret constitué de l'ensemble de ces points est ensuite post-traité.

Je souhaiterais améliorer ce programme en lui permettant de faire une analyse spectrale du signal. Le but étant de dégager les différentes périodes de mon signal (qui est forcément la somme de fonctions périodiques).

Ce que je demande n'a pas l'air très courant. En même temps, aujourd'hui avec des logiciels comme MatLab, Scilab, ou même Excel (il me semble), c'est beaucoup plus simple.

La programmation n'est pas mon fort, c'est juste un outil que j'utilise ponctuellement. Donc, si personne ne peut m'éclairer, je migrerai vers Excel.


Merci par avance!

ju

[Jerome Briot]

Sous MATLAB, la transformé de Fourier est effectuée grâce à la bibliothèque FFTW (www.fftw.org). Essaie de voir si tu ne peux pas "lier" cette bibliothèque à ton code ?

[Julien N]

Je ne suis malheuresement pas assez bon en informatique pour faire ça...

J'ai par contre écrit un bout de code fonctionnel réalisant la transformée de fourier discrete d'une liste:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
Function FDFT (X,N)
	Dim FX()
	TWOPI = 6.2831853071
	w = TWOPI/N
 
	for i = 0 to N-1
		Redim preserve FX(i)
		for j = 0 to N-1
			FX(i) = FX(i) + X(j)*cos(w*j*i)
		next
	next
	FDFT = FX
end function

C'est très basique et ne prend en compte que des nombres réels (de toute façon c'est ce qui m'intéresse)

Je peux aussi faire la transformée inverse:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
Function BDFT (FX,N)
	Dim X()
	TWOPI = 6.2831853071
	w = TWOPI/N
 
	for i = 0 to N-1
		Redim preserve X(i)
		for j = 0 to N-1
			X(i) = X(i) + (FX(j)*cos(w*j*i)/N
		next
	next
	BDFT = X
end function

Si je traite une sinusoïde, le résultat n'est pas trop mal. Par contre c'est lent...

Etant donnée que cette fonction sera intégrée un programme de post traitement plus gros et que j'aurais à l'appeler une bonne centaine de fois, j'aimerais coder une FFT en me basant sur la formulation de Cooley-Tukey que l'on peut trouver ici : http://en.wikipedia.org/wiki/Cooley%..._FFT_algorithm

Voici ce que j'ai fait :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
Function FFT (X,N,way)
	'with the Cooley-Tukey method, N must be a power of 2
	'if way = 1 the forward fourier transform is computed
	'if way = 0 the backward fourier transform is computed
 
	TWOPI = 6.2831853071
	Dim FX()
	Dim E()
	Dim O()
	'Compute the even and odd parts of the fourier transform
	for k = 0 to N-1
		redim preserve E(k)
		redim preserve O(k)
		for m = 0 to (N/2)-1
			c = cos(TWOPI*2*m*k/N)
			Ekm = Ekm + X(2*m)*c
			Okm = Okm + X(2*m+1)*c
		next
		E(k) = Ekm
		O(k) = Okm
	next
	'assembly
	for k = 0 to N-1
		redim preserve FX(k)
		if (way  = 1) then
			if (k<(N/2)) then
				FX(k) = E(k) + cos(TWOPI*k/N)*O(k)
			else
				FX(k) = E(k-(N/2)) + cos(TWOPI*(k*N/2)/N)*O(k-(N/2))
			end if
		else
			if (k<(N/2)) then
				FX(k) = (E(k) + cos(TWOPI*k/N)*O(k))/N
			else
				FX(k) = (E(k-(N/2)) + cos(TWOPI*(k*N/2)/N)*O(k-(N/2)))/N
			end if
		end if
	next
	FFT = FX
end function

Bien que plus rapide (heureusement), ça ne marche pas du tout. Les résultats sont incohérents : à savoir spectre discontinu sur toute la plage de fréquence (au lieu d'un seul pic pour une sinusoïde), et signal inversé à côté de la plaque..

Je ne sais pas si c'est moi qui fait une erreure, ou si je comprend mal l'algorithme.

[Jerome Briot]

Je ne sais pas si c'est moi qui fait une erreure, ou si je comprend mal l'algorithme.

Si ce n'est pas une erreur de syntaxe, ou si tu n'obtiens pas de réponse ici, tu pourras toujours poser ta question dans le forum Algorithmes > Mathématiques en mettant un lieu vers la présente discussion.