Scipy Spectrogram против нескольких Numpy FFT

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

def calc_old(raw_data):
    FFT_old = list()

    for i in range(0, len(raw_data), bf.WINDOW_STRIDE_LEN):
        if (i + bf.WINDOW_LEN) >= len(raw_data):
            # Skip the windows that would extend beyond the end of the data
            continue

        data_tmp = raw_data[i:i+bf.WINDOW_LEN]
        data_tmp -= np.mean(data_tmp)
        data_tmp = np.multiply(data_tmp, np.hanning(len(data_tmp)))
        fft_data_tmp = np.fft.fft(data_tmp, n=ZERO_PAD_LEN)
        fft_data_tmp = abs(fft_data_tmp[:int(len(fft_data_tmp)/2)])**2
        FFT_old.append(fft_data_tmp)

И новый код:

def calc_new(raw_data):
    data = np.array(raw_data)  # Required as the data is being handed in as a list
    f, t, FFT_new = spectrogram(data,
                                fs=60.0,
                                window="hann",
                                nperseg=bf.WINDOW_LEN,
                                noverlap=bf.WINDOW_OVERLAP,
                                nfft=bf.ZERO_PAD_LEN,
                                scaling='spectrum')

Таким образом, старый код окна временного ряда, удаляет среднее значение, применяет оконную функцию Ханна, принимает БПФ (с заполнением нулями, как ZERO_PAD_LEN>WINDOW_LEN), а затем берет абсолютное значение действительной половины и возводит его в квадрат для получения спектр мощности (Единицы V ** 2). Затем он сдвигает окно на WINDOW_STRIDE_LEN и повторяет процесс до тех пор, пока окно не выйдет за пределы конца данных. Это перекрытие WINDOW_OVERLAP.

Насколько я могу судить, спектрограмма должна делать то же самое с приведенными мной аргументами. Однако результирующие размеры БПФ различаются на 1 для каждой оси (например, старый код - MxN, новый код - (M + 1) x (N + 1)), а значение в каждом интервале частот существенно отличается - на несколько порядков. величины, в некоторых случаях.

Что мне здесь не хватает?


person esilk    schedule 17.08.2018    source источник
comment
Является ли один результат масштабированной версией другого? Может быть разница в нормализации. Но это не важно, важна форма вывода.   -  person Cris Luengo    schedule 17.08.2018


Ответы (2)


Масштабирование

Реализация в calc_old использует вывод из np.fft.fft напрямую без масштабирования.

С другой стороны, реализация calc_new использует scipy.signal.spectrogram который в конечном итоге использует _5 _ но также масштабирует результаты на основе полученных аргументов scaling и return_onesided. Более конкретно:

  • Для значения по умолчанию return_onesided=True (поскольку вы не указали явное значение в calc_new), значение каждой ячейки удваивается для подсчета общей энергии, включая симметричную ячейку.
  • Для предоставленного scaling='spectrum' значения дополнительно масштабируются с коэффициентом 1.0/win.sum()**2. Для выбранного окна Ханна это соответствует 4/N**2, где N=bf.WINDOW_LEN - длина окна.

Таким образом, вы можете ожидать, что новая реализация cald_new даст вам результат, который масштабируется на 8/bf.WINDOW_LEN**2 общий коэффициент по сравнению с calc_old. В качестве альтернативы, если вы хотите, чтобы ваша вторая реализация давала то же масштабирование, что и calc_old, вам следует умножить результат scipy.signal.spectrogram на 0.125 * bf.WINDOW_LEN**2.

Количество интервалов частоты

Учитывая четное число точек nperseg, ваша первоначальная реализация calc_old сохраняет только nperseg//2 интервалы частоты.

С другой стороны, полный неизбыточный половинный спектр должен дать вам nperseg//2 + 1 частотный бин (есть nperseg-2 бина с соответствующей симметрией плюс 2 асимметричных бина при 0 Гц и частоте Найквиста, поэтому сохранение неизбыточной части оставит вас с (nperseg-2)//2 + 2 == nperseg//2 + 1). Вот что возвращает scipy.signal.spectrogram.

Другими словами, в вашей первоначальной реализации calc_old отсутствует частотный интервал Найквиста.

Количество временных шагов

Реализация в calc_old пропускает последний временной шаг, если для вычисления последнего временного шага осталось менее bf.WINDOW_LEN выборок. Эти образцы не пропускаются только тогда, когда len(raw_data)-bf.WINDOW_STRIDE_LEN является точным кратным bf.WINDOW_LEN. Я предполагаю, что это не относится к вашей конкретной входной последовательности.

Напротив, scipy.signal.spectrogram дополняет данные дополнительными выборками, если это необходимо, так что все входные выборки используются во время вычисления спектрограммы, и это может привести к одному дополнительному временному шагу по сравнению с вашей реализацией calc_old.

person SleuthEye    schedule 18.08.2018
comment
Это подтверждает все, что мы с коллегой подозревали и на что смотрели, особенно в отношении заполнения временного шага, которое я нигде не вижу как задокументированное? Теперь мы получаем совпадающие результаты (и намного быстрее), спасибо! - person esilk; 20.08.2018

Наверное, кто-то мог знать причину, почему я получаю небольшое расхождение при сравнении результатов спектрограммы с ручным БПФ?

# Parametrs of signal and its preprocessing
sample_rate = 4
window_size = 512 * sample_rate
detrend = 'linear'
tukey_alpha = 0.25
[![enter image description here][1]][1]
# Spectrogram
f, t, S = scipy.signal.spectrogram(signal, sample_rate, nperseg=window_size, noverlap=sample_rate, scaling='spectrum', mode='magnitude', detrend=detrend)

# FFT on the leftmost window of signal
windowed_signal = signal[:window_size]
windowed_signal = scipy.signal.detrend(windowed_signal, type=detrend)
windowed_signal *= scipy.signal.windows.tukey(window_size, tukey_alpha)
A = np.fft.rfft(windowed_signal)
freqs = np.fft.rfftfreq(window_size) * sample_rate
positive_freqs_n = int(np.ceil(window_size / 2.))
freqs_slice = slice(0, positive_freqs_n)
magnitudes = np.abs(A)[freqs_slice] / window_size

# Plotting
plt.plot(f, S.T[0], label='First window from scipy.signal.spectrogram')
plt.plot(freqs[freqs_slice], magnitudes, alpha=0.5, label='np.fft on first window')
plt.yscale('log')
plt.xscale('log')
plt.legend()
plt.xlabel('Frequency, Hz')
plt.ylabel('Magnitude')

введите описание изображения здесь

Откуда могло возникнуть это небольшое вертикальное несоответствие? Спасибо!

person Valentin Melnikov    schedule 21.05.2020