Я использовал функцию fft
в numpy, которая привела к сложному массиву. Как получить точные значения частоты?
Как извлечь частоту, связанную с значениями fft в python
Ответ 1
np.fft.fftfreq
указывает частоты, связанные с коэффициентами:
import numpy as np
x = np.array([1,2,1,0,1,2,1,0])
w = np.fft.fft(x)
freqs = np.fft.fftfreq(len(x))
for coef,freq in zip(w,freqs):
if coef:
print('{c:>6} * exp(2 pi i t * {f})'.format(c=coef,f=freq))
# (8+0j) * exp(2 pi i t * 0.0)
# -4j * exp(2 pi i t * 0.25)
# 4j * exp(2 pi i t * -0.25)
ОП спрашивает, как найти частоту в Герце.
Я считаю, что формула frequency (Hz) = abs(fft_freq * frame_rate)
.
Вот какой код, который демонстрирует это.
Сначала мы создаем волновой файл с частотой 440 Гц:
import math
import wave
import struct
if __name__ == '__main__':
# http://stackoverflow.com/questions/3637350/how-to-write-stereo-wav-files-in-python
# http://www.sonicspot.com/guide/wavefiles.html
freq = 440.0
data_size = 40000
fname = "test.wav"
frate = 11025.0
amp = 64000.0
nchannels = 1
sampwidth = 2
framerate = int(frate)
nframes = data_size
comptype = "NONE"
compname = "not compressed"
data = [math.sin(2 * math.pi * freq * (x / frate))
for x in range(data_size)]
wav_file = wave.open(fname, 'w')
wav_file.setparams(
(nchannels, sampwidth, framerate, nframes, comptype, compname))
for v in data:
wav_file.writeframes(struct.pack('h', int(v * amp / 2)))
wav_file.close()
Это создает файл test.wav
.
Теперь мы читаем в данных, FFT, находим коэффициент с максимальной мощностью,
и найти соответствующую частоту fft, а затем преобразовать в Hertz:
import wave
import struct
import numpy as np
if __name__ == '__main__':
data_size = 40000
fname = "test.wav"
frate = 11025.0
wav_file = wave.open(fname, 'r')
data = wav_file.readframes(data_size)
wav_file.close()
data = struct.unpack('{n}h'.format(n=data_size), data)
data = np.array(data)
w = np.fft.fft(data)
freqs = np.fft.fftfreq(len(w))
print(freqs.min(), freqs.max())
# (-0.5, 0.499975)
# Find the peak in the coefficients
idx = np.argmax(np.abs(w))
freq = freqs[idx]
freq_in_hertz = abs(freq * frate)
print(freq_in_hertz)
# 439.8975
Ответ 2
Частоты, связанные с значениями DFT (в python)
По fft, Быстрому преобразованию Фурье, мы понимаем элемент большого семейства алгоритмов, которые позволяют быстро вычислять ДПФ, Дискретное преобразование Фурье, равносильного сигнала.
A DFT преобразует список N комплексных чисел в список из N комплексных чисел, при том понимании, что оба списка являются периодическими с периодом N.
Здесь мы имеем дело с реализацией numpy
fft.
Во многих случаях вы думаете о
- сигнал x, определенный во временной области длины N, отбираемый на постоянный интервал dt,
- его DFT X (здесь конкретно
X = np.fft.fft(x)
), элементы которого отбираются на частотной оси с частотой дискретизации dw.
Некоторое определение
-
период (aka duration) сигнала
x
, выбранный вdt
сN
образцами, равенT = dt*N
-
основные частоты (в Гц и в рад/с)
x
, ваш DFTdf = 1/T dw = 2*pi/T # =df*2*pi
-
верхняя частота - это Nyquist частота
ny = dw*N/2
(и это не
dw*N
)
Частоты, связанные с конкретным элементом в DFT
Частоты, соответствующие элементам в X = np.fft.fft(x)
для заданного индекса 0<=n<N
, можно вычислить следующим образом:
def rad_on_s(n, N, dw):
return dw*n if n<N/2 else dw*(n-N)
или в одной развертке
w = np.array([dw*nif n<N/2 else dw*(n-N) for n in range(N)])
если вы предпочитаете рассматривать частоты в Гц, s/w/f/
f = np.array([df*n if n<N/2 else df*(n-N) for n in range(N)])
Используя эти частоты
Если вы хотите изменить исходный сигнал x
→ y
, применяя оператор в частотной области только в форме функции только частоты, то вычисляем w
и
Y = X*f(w)
y = ifft(Y)
Представляем np.fft.fftfreq
Конечно, numpy
имеет функцию удобства np.fft.fftfreq
, которая возвращает безразмерные частоты, а не мерные, но это так же просто, как
f = np.fft.fftfreq(N)*N*df
w = np.fft.fftfreq(N)*N*dw
Ответ 3
Частота - это только индекс массива. При индексе n частота равна 2 πn/длина массива (радианы на единицу). Рассмотрим:
>>> numpy.fft.fft([1,2,1,0,1,2,1,0])
array([ 8.+0.j, 0.+0.j, 0.-4.j, 0.+0.j, 0.+0.j, 0.+0.j, 0.+4.j,
0.+0.j])
результат имеет ненулевые значения при индексах 0, 2 и 6. Имеется 8 элементов. Это означает
2πit/8 × 0 2πit/8 × 2 2πit/8 × 6
8 e - 4i e + 4i e
y ~ ———————————————————————————————————————————————
8