Scipy / Numpy FFT Частотный анализ

Question

Я ищу, как превратить ось частоты в БПФ (взятую через scipy.fftpack.fftfreq) в частоту в герцах, а не в двоичные или дробные бины.

Я пытался кодировать нижедля проверки БПФ:

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

Частота дискретизации должна составлять 4000 выборок / 120 секунд = 33,34 выборки / сек.

Сигнал имеет сигнал 2,0 Гц, сигнал 8,0 Гц,и некоторый случайный шум.

Я беру БПФ, беру частоты и строю график.Цифры довольно бессмысленные.Если я умножу частоты на 33,34 (частоту дискретизации), то получу пики примерно на 8 Гц и 15 Гц, что кажется неправильным (также частоты должны быть с коэффициентом 4, а не 2!).

Есть мысли о том, что я здесь делаю не так?

Answer 1

Я думаю, вам не нужно выполнять fftshift (), и вы можете передать период выборки в 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()

Из графика видно, что есть два пика при 2 Гц и 8 Гц.

Answer 2

scipy.fftpack.fftfreq (n, d) дает вам частоты напрямую. Если вы установите d=1/33.34, это скажет вам частоту в Гц для каждой точки БПФ.

Answer 3

Ширина частоты каждого бина равна (sampling_freq / num_bins).

Более фундаментальная проблема заключается в том, что вашей частоты дискретизации недостаточно для ваших сигналов интереса. Ваша частота дискретизации составляет 8,3 Гц; вам нужно по крайней мере 16 Гц для захвата входного тона 8 Гц. ¹

---

_{1. Для всех экспертов DSP; Я знаю, что это действительно BW, а не максимальная частота. Но я предполагаю, что ОП не хочет делать сбор данных с низкой выборкой.}

Answer 4

Ваше уравнение испорчено.

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

Это дает вам 4000 выборок вашей функции с частотой дискретизации 33,33 Гц, что составляет 120 секунд данных.

Теперь возьми свой БПФ. Корзина 0 будет содержать результат DC. Корзина 1 будет 33,33, корзина 2 - 66,66 и т. Д.

Редактировать: я забыл упомянуть, что, поскольку ваша частота дискретизации составляет 33,33 Гц, максимальная частота, которая может быть представлена, будет fs / 2 или 16,665 Гц.