Analyse de fréquence FFT Scipy/Numpy

Question

Je cherche comment transformer l'axe de fréquence dans un fft (pris via scipy.fftpack.fftfreq) en une fréquence en Hertz, plutôt qu'en bins ou bins fractionnés.

J'ai essayé le code ci-dessous pour tester la FFT :

t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
FFT = scipy.fftpack.fftshift(FFT)
freqs = scipy.fftpack.fftfreq(signal.size)

pylab.plot(freqs,FFT,'x')
pylab.show()

Le taux d'échantillonnage doit être de 4000 échantillons / 120 secondes = 33,34 échantillons/sec.

Le signal comporte un signal de 2,0 Hz, un signal de 8,0 Hz et un bruit aléatoire.

Je prends la FFT, je prends les fréquences et je les trace. Les chiffres sont assez absurdes. Si je multiplie les fréquences par 33,34 (la fréquence d'échantillonnage), j'obtiens des pics à environ 8 Hz et 15 Hz, ce qui semble erroné (de plus, les fréquences devraient être séparées par un facteur de 4, et non de 2 !)

Vous avez une idée de ce que je fais de mal ?

Answer 1

Je pense que vous n'avez pas besoin de faire fftshift(), et que vous pouvez passer la période d'échantillonnage à fftfreq() :

import scipy
import scipy.fftpack
import pylab
from scipy import pi
t = scipy.linspace(0,120,4000)
acc = lambda t: 10*scipy.sin(2*pi*2.0*t) + 5*scipy.sin(2*pi*8.0*t) + 2*scipy.random.random(len(t))

signal = acc(t)

FFT = abs(scipy.fft(signal))
freqs = scipy.fftpack.fftfreq(signal.size, t[1]-t[0])

pylab.subplot(211)
pylab.plot(t, signal)
pylab.subplot(212)
pylab.plot(freqs,20*scipy.log10(FFT),'x')
pylab.show()

A partir du graphique, vous pouvez voir qu'il y a deux pics à 2Hz et 8Hz.

Answer 2

scipy.fftpack.fftfreq(n, d) vous donne directement les fréquences. Si vous réglez d=1/33.34 Ceci vous indiquera la fréquence en Hz pour chaque point du fft.

Answer 3

La largeur de fréquence de chaque bin est (sampling_freq / num_bins).

Un problème plus fondamental est que votre fréquence d'échantillonnage n'est pas suffisante pour vos signaux d'intérêt. Votre fréquence d'échantillonnage est de 8,3 Hz ; vous avez besoin d'au moins 16 Hz pour capturer un signal d'entrée de 8 Hz. 1

~~

~~

_{1. Pour tous les experts en DSP, je suis conscient que c'est la largeur de bande qui est importante, et non la fréquence maximale. Mais je suppose que l'OP ne veut pas faire d'acquisition de données sous-échantillonnées.}

Answer 4

Votre équation est erronée.

fs = 33.33
df1 = 2*pi * (2.0/fs)
df2 = 2*pi * (5.0/fs)
x = [10*sin(n*df1) + 5*sin(n*df2) + 2*random.random() for n in range(4000)]

Vous obtenez ainsi 4000 échantillons de votre fonction, échantillonnés à 33,33 Hz, ce qui représente 120 secondes de données.

Maintenant, prenez votre FFT. La case 0 contiendra le résultat du DC. La case 1 sera 33.33, la case 2 sera 66.66, etc

Edit : J'ai oublié de mentionner que, puisque votre taux d'échantillonnage est de 33,33 Hz, la fréquence maximale qui peut être représentée sera fs/2, soit 16,665 Hz.