Почему значения частоты округляются в сигнале с использованием БПФ?

Question

Итак, я пытаюсь выяснить, как использовать DFT на практике для обнаружения преобладающих частот в сигнале.Я пытался обдумать, что такое преобразования Фурье и как работают алгоритмы DFT, но, видимо, у меня еще есть пути.Я написал некоторый код для генерации сигнала (так как целью было работать с музыкой, я сгенерировал мажорный аккорд С, отсюда и странные значения частоты), а затем попытался вернуться к частотам.Вот код, который у меня есть

sr = 44100 # sample rate
x = np.linspace(0, 1, sr) # one second of signal
tpi = 2 * np.pi
data = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)
freqs = np.fft.fftfreq(sr)
fft = np.fft.fft(data)
idx = np.argsort(np.abs(fft))
fft = fft[idx]
freqs = freqs[idx]
print(freqs[-6:] * sr)

Это дает мне [-262. 262. -330. 330. -392. 392.] , что отличается от частот, которые я кодировал (261.63, 329.63 и 392.0).Что я делаю не так и как мне это исправить?

Answer 1

Результирующие ячейки DFT разделены по частоте Fs / N, где N - длина FFT.Таким образом, продолжительность вашего окна DFT ограничивает разрешение в единицах частотного интервала между бинами результатов DFT.

Но для хорошо разнесенных частотных пиков в малошумных (высокие S / N) вместо увеличения длительностиданные, вы можете вместо этого оценить местоположения пиков частоты с более высоким разрешением путем интерполяции результата DFT между ячейками результатов DFT.Вы можете попробовать параболическую интерполяцию для грубой оценки местоположения пика частоты, но оконная интерполяция Синка (по существу, реконструкция Шеннона-Уиттекера) обеспечит гораздо лучшую точность и разрешение оценки частоты (при достаточно низком минимальном уровне шума вокруг интересующего пика (ей) частоты),например, нет синусоиды поблизости в вашем случае искусственной формы волны).

Answer 2

Действительно, если кадр длится T секунд, частоты ДПФ составляют k/T Гц, где k - целое число.Как следствие, передискретизация не улучшает точность расчетной частоты, если эти частоты идентифицированы как максимумы величины ДПФ.Напротив, рассмотрение более длинных кадров продолжительностью 100 с приведет к разнесению между частотами DFT 0,01 Гц, что может быть достаточно для получения ожидаемой частоты. Можно добиться гораздо лучшего, оценив частоту пика как его среднюю частоту с учетом плотности мощности.

Рисунок 1: дажепосле применения окна Тьюки ДПФ оконного сигнала не является суммой Дирака: все еще есть некоторая спектральная утечка в нижней части пиков.Эта мощность должна учитываться при оценке частот.

Другая проблема заключается в том, что длина кадра не кратна периоду сигнала, который в любом случае не может быть периодическим.Тем не менее, ДПФ вычисляется так, как если бы сигнал был периодическим, но прерывистым на границе кадра.Он вызывает скачкообразные частоты, описываемые как спектральная утечка .Оконное управление является эталонным методом для решения таких проблем и смягчения проблемы, связанной с искусственным разрывом.Действительно, значение окна непрерывно уменьшается до нуля вблизи краев кадра. Существует список оконных функций , и множество оконных функций доступно в scipy.signal .Окно применяется как:

tuckey_window=signal.tukey(len(data),0.5,True)
data=data*tuckey_window

. В этот момент частоты, наибольшие значения которых по-прежнему составляют 262, 330 и 392. Применение окна только делает пики более заметными: ДПФ оконного сигнала обладаеттри выделенных пика, каждый из которых имеет центральный лепесток и боковые лепестки, в зависимости от DFT окна. лепестки этих окон симметричны: поэтому центральная частота может быть вычислена как средняя частота пика относительно плотности мощности.

import numpy as np
from scipy import signal
import scipy

sr = 44100 # sample rate
x = np.linspace(0, 1, sr) # one second of signal
tpi = 2 * np.pi
data = np.sin(261.63 * tpi * x) + np.sin(329.63 * tpi * x) + np.sin(392.00 * tpi * x)

#a window...
tuckey_window=signal.tukey(len(data),0.5,True)
data=data*tuckey_window

data -= np.mean(data)
fft = np.fft.rfft(data, norm="ortho")

def abs2(x):
        return x.real**2 + x.imag**2

fftmag=abs2(fft)[:1000]
peaks, _= signal.find_peaks(fftmag, height=np.max(fftmag)*0.1)
print "potential frequencies ", peaks

#compute the mean frequency of the peak with respect to power density
powerpeak=np.zeros(len(peaks))
powerpeaktimefrequency=np.zeros(len(peaks))
for i in range(1000):
    dist=1000
    jnear=0
    for j in range(len(peaks)):
        if dist>np.abs(i-peaks[j]):
             dist=np.abs(i-peaks[j])
             jnear=j
    powerpeak[jnear]+=fftmag[i]
    powerpeaktimefrequency[jnear]+=fftmag[i]*i


powerpeaktimefrequency=np.divide(powerpeaktimefrequency,powerpeak)
print 'corrected frequencies', powerpeaktimefrequency

Результирующие расчетные частоты составляют 261,6359Гц, 329,637 Гц и 392,0088 Гц: он намного лучше, чем 262, 330 и 392 Гц, и удовлетворяет требуемой точности 0,01 Гц для такого чистого бесшумного входного сигнала.

Answer 3

Так как вы хотите получить разрешение 0,01 Гц, вам нужно будет сэмплировать данные не менее чем за 100 секунд.Вы сможете разрешать частоты примерно до 22,05 кГц.