Как извлечь частоту, связанную со значениями FFT в Python - PullRequest
Новая
       66

Как извлечь частоту, связанную со значениями FFT в Python

35 голосов
/ 12 сентября 2010

Я использовал функцию fft в numpy, что привело к сложному массиву.Как получить точные значения частоты?

Ответы [ 3 ]

57 голосов
/ 12 сентября 2010

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__':
    # /2963226/kak-pisat-stereo-wav-faily-v-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. Теперь мы читаем в данных, БПФ это, найти коэффициент с максимальной мощностью, и найти соответствующую частоту БПФ, а затем преобразовать в герц: 

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
30 голосов
/ 28 ноября 2014

Частоты, связанные со значениями DFT (в питоне)

Под fft , быстрым преобразованием Фурье, мы понимаем член большого семейства алгоритмов, которые позволяют быстрое вычисление ДПФ, дискретного преобразования Фурье, равноосэмплированного сигнала.

A DFT преобразует список N комплексных чисел в список N комплексных чисел с пониманием, что оба списка являются периодическими с периодом N .

Здесь мы имеем дело с numpy реализацией fft .

Во многих случаях вы думаете о

  • сигнал x , определенный во временной области длины N , сэмплированный на постоянный интервал дт ,
  • его DFT X (здесь конкретно X = np.fft.fft(x)), элементы которого дискретизируются по оси частот с частотой дискретизации dw .

Некоторое определение

  • период (он же длительность) сигнала x, сэмплированный при dt с N сэмплов равен 

    T = dt*N

  • основные частоты (в Гц и в рад / с) X, ваше ДПФ равно 

    df = 1/T
dw = 2*pi/T # =df*2*pi

  • верхняя частота Найквиста частота 

    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
5 голосов
/ 12 сентября 2010

Частота - это просто индекс массива.При индексе 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
...