Scipy/Numpy FFT Analiza częstotliwości

Question

Szukam, jak obrócić osi częstotliwości w fft (podjęte przez scipy.fftpack.fftfreq) na częstotliwość w hercach, a nie pojemniki lub frakcje częściowe.

że próbował kodu poniżej przetestować 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() 

Częstotliwość próbkowania powinna 4000 120 próbek/s = 33,34 próbek/sek.

Sygnał ma sygnał 2,0 Hz, sygnał 8,0 Hz i niektóre zakłócenia losowe.

Podejmuję FFT, łapię częstotliwości i kreślę. Liczby są dość bezsensowne. Jeśli pomnożę częstotliwości przez 33,34 (częstotliwość próbkowania), otrzymuję szczyty o wartości około 8 Hz i 15 Hz, co wydaje się błędne (również częstotliwości powinny być 4-krotnie oddalone, a nie 2!).

Jakieś myśli o tym, co robię źle tutaj?

Answer 1

myślę, że nie trzeba robić fftshift(), i można przejść okres próbkowania do 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() 

z wykresu widać istnieją dwa pik przy 2 Hz i 8 Hz.

Answer 2

Szerokość częstotliwości każdego pojemnika to (sampling_freq/num_bins).

Bardziej fundamentalnym problemem jest to, że częstotliwość próbkowania nie jest wystarczająca dla sygnałów będących przedmiotem zainteresowania. Twoja częstotliwość próbkowania wynosi 8,3 Hz; potrzebujesz co najmniej 16 Hz, aby uchwycić ton wejściowy 8Hz.

---

_{1. Do wszystkich ekspertów DSP; Mam świadomość, że to rzeczywiście BW ma znaczenie, a nie maksymalna częstotliwość. Ale zakładam, że OP nie chce wykonać nieskompresowanego akwizycji danych.}

Answer 3

Twoje równanie jest pomieszane.

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)] 

Daje to 4000 próbek swojej funkcji, próbkowany 33,33 Hz, reprezentujących 120 sekund danych.

Teraz weź swoją FFT. Kosz 0 będzie utrzymywać wynik DC. Bin 1 będzie 33,33, bin 2 będzie 66,66, itp.

Edytuj: Zapominam wspomnieć, że ponieważ częstotliwość próbkowania wynosi 33,33 Hz, maksymalna częstotliwość, która może być reprezentowana będzie fs/2, lub 16.665 Hz.

Answer 4

scipy.fftpack.fftfreq(n, d) daje częstotliwości bezpośrednio. Jeśli ustawisz d=1/33.34, to powie ci częstotliwość w Hz dla każdego punktu fft.