[NumPy] Fonction fft() issue de numpy/numarray

Version imprimable

Bonjour,

En complément du topic traitant de la création du .wav monofréquence, je passe actuellement à l'étape clé de la FFT. Grace à l'aide de oiffrig , la création du fichier .wav se fait de la sorte:
Code:

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 import math, numarray as NA, wave, numarray as NA2 freq = 200. freqech = 44100 duree = 10. omega = 2 * math.pi * freq son = NA.arange(0., duree, 1./freqech, NA.Float64) son = NA.sin(omega * son) son *= 32768. son = son.astype(NA.Int16) wavfile = wave.open('AutreAmplitude200.wav','w') wavfile.setparams((1, 2, freqech , len(son), 'NONE', 'not compressed')) wavfile.writeframes(son.tostring()) wavfile.close() Valeurs = son[0:222]
Dans la varable Valeur, je sauvegarde les valeurs prises par la fonction sinus (avec f = 200Hz) sur 1 période environ.
En appliquant la fonction fft issue du module numarray, je devrais donc obtenir quelque chose me disant qu'effectivement j'ai un signal de l'ordre de 200Hz.
Code:

1 2 from numarray import fft print fft.fft(Valeurs,450).real
Or le résultat obtenu est le suivant:
Code:

1 2 3 4 5 6 7 8 9 3.51000000e+02 3.02301824e+06 2.29319478e+05 -1.96406314e+06 -1.23762282e+04 -4.48999295e+05 -1.02429368e+04 -2.03945581e+05 -9.59045653e+03 -1.15423303e+05 -9.26517011e+03 -7.30000882e+04 -9.02794648e+03 -4.93465665e+04 -8.82407680e+03 -3.47937512e+04 -8.63322297e+03 -2.52527407e+04 -8.44808818e+03 -1.86411046e+04 -8.24596904e+03 -1.39443011e+04 -8.04285242e+03 -1.04709028e+04 -7.82704600e+03 -7.87677142e+03 -7.59881077e+03 -5.88953325e+03 -7.36427447e+03 -4.36747628e+03 -7.11766236e+03 -3.17780670e+03 ETC...
La première valeur correspond à la composante continue du signal, autrement sa valeur moyenne, OK
Mais Quid des autres valeurs? Amplitudes de quelles harmoniques, fréquences(cos, sin)? Pourquoi des valeurs négatives?

Si vous avez des suggestions....

19/06/2007, 15h21
Svart26

Bon, j'ai debusqué ça sur le net:

Citation:

Results of FFT are in the form of complex vector of the same length as
data. Starting from position 0 it gives you an constant (DC factor), position
1 an amplitude and phase (remember - complex number gives you both amplitude
and phase) of wave of the length table/2 and so on. If you take real value of
this, you discard part of the information. AFAIR - taking real part gives you
sine component, taking imaginary part gives you cosine component, taking
absolute value gives you total amplitude and taking angle gives you phase

Autrement, pour récuperer l'amplitude des sinusoidale, il suffit de prendre le module du comple correspondant, mais quand à connaître la fréquence qui lui correspond...:roll:
19/06/2007, 15h52
Svart26

A titre de simplication:

la fft de [ 1. 0. 1. 0.] avec 8 samples donne:
[ 2.+0.j 1.-1.j 0.+0.j 1.+1.j 2.+0.j 1.-1.j 0.+0.j 1.+1.j]

la fft de [ 1. 0. 1. 0.] avec 4 samples donne:
[2 0 2 0]

Si quelqu'un comprend le résultat... :mrgreen:

:boulet:
19/06/2007, 18h08
parp1

Et bien oui c'est simple:

Avec 8 samples tu as donc 8 echantillons
Avex 4 Samples tu en as donc 4 echantillons. (oui je sais que jusqu'ici je sais traduire.)

Mais lorsque tu regardes ton resultat avec huit sample et que tu prends un echantillon sur deux.

Ca te donne :

Code:

[2.+0.j 1.-1.j 0.+0.j 1.+1.j 2.+0.j 1.-1.j 0.+0.j 1.+1.j]

Ce qui donne:

Code:

2+0j,0+0j,2+0j,0+0j ou plus simplement 2,0,2,0
19/06/2007, 18h10
Svart26

Citation:

Envoyé par parp1

Et bien oui c'est simple:

Avec 8 samples tu as donc 8 echantillons
Avex 4 Samples tu en as donc 4 echantillons. (oui je sais que jusqu'ici je sais traduire.)

Mais lorsque tu regardes ton resultat avec huit sample et que tu prends un echantillon sur deux.

Ca te donne :

Code:

[2.+0.j 1.-1.j 0.+0.j 1.+1.j 2.+0.j 1.-1.j 0.+0.j 1.+1.j]

Ce qui donne:

Code:

2+0j,0+0j,2+0j,0+0j ou plus simplement 2,0,2,0

Oui, tout à fait. Mais plus précisément, je cherche à savoir à quelles fréquences correspondent ces amplitudes :aie:
21/06/2007, 10h50
fabrice-102

Pour récupérer les fréquence, sous numpy, il suffit de faire :
freq = fftfreq(nbr,d=eps)

avec nbr le nombre d'échantillons et eps l'inverse du nombre d'échantillons récupérés par seconde. Pour vous : 1.0/44100 (ne pas oublie 1.0 pour travailler avec une précision suffisante !!)

Vous pouvez ensuite, par exemple, calculer la fréquence dont l'énergie est la plus importante (soit pour une sinusoïde, la fréquence de la sinusoïde) :
maxfreq = max(freq)

et vous pouvez ainsi calculer l'énergie correspondante :
fft = fft(Valeurs,len(Valeurs))
energie = fft[argmax(freq)]
24/06/2007, 22h15
Svart26

Citation:

Envoyé par fabrice-102

Pour récupérer les fréquence, sous numpy, il suffit de faire :
freq = fftfreq(nbr,d=eps)

avec nbr le nombre d'échantillons et eps l'inverse du nombre d'échantillons récupérés par seconde. Pour vous : 1.0/44100 (ne pas oublie 1.0 pour travailler avec une précision suffisante !!)

Vous pouvez ensuite, par exemple, calculer la fréquence dont l'énergie est la plus importante (soit pour une sinusoïde, la fréquence de la sinusoïde) :
maxfreq = max(freq)

et vous pouvez ainsi calculer l'énergie correspondante :
fft = fft(Valeurs,len(Valeurs))
energie = fft[argmax(freq)]

Merci beaucoup, tout paraît évident après lecture de votre réponse! :D

Citation:

Envoyé par fabrice-102

Pour récupérer les fréquence, sous numpy, il suffit de faire :
freq = fftfreq(nbr,d=eps)

avec nbr le nombre d'échantillons et eps l'inverse du nombre d'échantillons récupérés par seconde. Pour vous : 1.0/44100 (ne pas oublie 1.0 pour travailler avec une précision suffisante !!)

Vous pouvez ensuite, par exemple, calculer la fréquence dont l'énergie est la plus importante (soit pour une sinusoïde, la fréquence de la sinusoïde) :
maxfreq = max(freq)

et vous pouvez ainsi calculer l'énergie correspondante :
fft = fft(Valeurs,len(Valeurs))
energie = fft[argmax(freq)]

Bien, donc petit bilan des essais que j'ai effectué aujourd'hui, histoire de donner un retour des choses:

Effectivement, fftfreq renvoie les fréquences, et je dois dire que cette fonction fait vraiment plaisir.

Mon besoin étant de déterminer la/les fréquences où l'amplitude est maximale, je réalise les étapes suivantes:
Code:

1 2 3 4 5 6 7 8 9 10 11 FFT = fft.fft(valeurs, 2048) FFTabs = [] for i in range(50): FFTabs.append(abs(FFT[i])) #calcul de la valeur absolu des amplitudes (*) (pour les 50 premiers echantillons ici) maxi = max(FFTabs) indice = FFTabs.index(maxi) freq = numpy.fft.fftfreq(2048,1.0/44100) print 'frequence avec amplitude max: ', freq[indice]
(*)d'apres ce lien, l'amplitude totale de chacune des frequences obtenues lors de la decomposition est obtenue en prenant la valeur absolue du complexe . http://mail.python.org/pipermail/tut...st/040263.html

Bref, d'après mes premiers essais, tout est K, j'espère que ca ne changera pas...
Merci ;)

26/06/2007, 20h17
Matthieu Brucher

A la place de numarray, commence par utiliser numpy complètement, ça résoudra peut-être quelques soucis !
Ensuite, si le retour est complexe, tu prends la valeur absolue de tes termes, et avec l'échelle qui t'a été donné, c'est bon.

Autre détail, la FFT fournie avec FFTW ne possède pas de terme correcteur, cf mon tuto : http://miles.developpez.com/tutoriels/algo/fft/
26/06/2007, 23h28
Svart26

Citation:

Envoyé par Miles

A la place de numarray, commence par utiliser numpy complètement, ça résoudra peut-être quelques soucis !
Ensuite, si le retour est complexe, tu prends la valeur absolue de tes termes, et avec l'échelle qui t'a été donné, c'est bon.

Autre détail, la FFT fournie avec FFTW ne possède pas de terme correcteur, cf mon tuto : http://miles.developpez.com/tutoriels/algo/fft/

Sympathique tutoriel. Sinon oui en effet, avec l'amplitude "valeur absolue" du complexe et l'echelle de frequence c'est parfait, le tout complété par un graphique sous Tkinter est le résultat est ma foi fort sympathique. (même si les valeurs des amplitudes rendues par FFt() me tourmentent un peu).
Bref, demain je jetterai également un oeil à FFTW, Thanks.