Простой расчёт и визуализация FFT в Python с использованием UliEngineering

UliEngineering — это библиотека смешанной аналитики данных в Python — одна из утилит, которую она предоставляет, — это простой в использовании пакет для вычисления FFT. В отличие от других пакетов, эта библиотека ориентирована на практическое применение и позволяет вам выполнить FFT всего в одной строке кода! Знание математики не требуется. Сначала установите UliEngineering.

Начало работы

Сначала мы сгенерируем некоторые тестовые данные. См. этот предыдущий пост для более подробной информации о том, как генерировать синусоидальные тестовые данные:

from UliEngineering.SignalProcessing.Simulation import *

# Сгенерировать тестовые данные: тон 100 Гц + 400 Гц
data = sine_wave(frequency=100.0, samplerate=1000, amplitude=1.) \
       + sine_wave(frequency=400.0, samplerate=1000, amplitude=0.5)

Тестовые данные состоят из синусоиды 100 Hz плюс синусоиды 400 Hz (с половинной амплитудой). Этот сигнал сэмплирован с частотой дискретизации 1000 Hz.

Теперь мы можем вычислить и визуализировать FFT, используя matplotlib:

# Вычислить FFT. Используйте ту же частоту дискретизации
# ПРИМЕЧАНИЕ: Окно по умолчанию "blackman"!
from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Построить график
from matplotlib import pyplot as plt
plt.style.use("ggplot")

plt.gcf().set_size_inches(10, 5) # Используйте (20, 10) для получения большего графика
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Частота")
plt.ylabel("Амплитуда")

compute_fft(data, samplerate=1e3) возвращает объект FFT, содержащий поля, такие как frequencies, amplitude и phase. Он выполняет FFT размером входных данных (т.е. поскольку data — массив длиной 1000, FFT будет размера 1000).

fft.frequencies — массив частот (в Hz), соответствующих значениям в fft.amplitudes. Вы также можете использовать fft.angles, чтобы получить относительные углы в градусах, но это не рассматривается в этом посте.

Как видно на графике выше, максимальная частота, которая может быть обнаружена, всегда половина частоты дискретизации, т.е. для нашей частоты дискретизации $f_s = 1000,\text{Hz}$ это $500,\text{Hz}$. См. это более математическое объяснение FFT, если вы хотите узнать больше деталей.

Внутренне compute_fft() выполняет это вычисление:

• $2$ — поправочный коэффициент, учитывающий, что мы отбрасываем вторую половину сырого результата FFT (поскольку мы делаем real FFT)
• $\frac{1}{\text{len(data)}}$ нормализует результаты FFT, чтобы они были независимы от длины данных (т.е. если вы передадите более длинный сэмпл той же синусоиды, вы всё равно получите тот же результат)
• $\text{abs}\left(\cdots\right)$ преобразует комплексный фазово-зависимый результат FFT в спектр, который легче читать и визуализировать.
• Window (по умолчанию blackman) — это окно, которое применяется к данным для смягчения некоторых математических эффектов в начале и конце набора данных. См. статью в Википедии о оконных функциях для более подробной информации. UliEngineering в настоящее время предлагает этот список оконных функций:
  • blackman
  • bartlett
  • hamming
  • hanning
  • kaiser (параметр фиксирован на 2.0)
  • none

Выбор диапазонов частот

Используя API UliEngineering, выбрать диапазон частот FFT тривиально просто: просто используйте fft[lowfreq:highfreq]. Вы можете использовать fft[lowfreq:], чтобы выбрать всё начиная с lowfreq, или использовать fft[:highfreq], чтобы выбрать всё до highfreq.

from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Выбрать диапазон частот: от 50 до 200 Гц
fft = fft[50.0:200.0]

# Построить график
from matplotlib import pyplot as plt
plt.style.use("ggplot")

plt.gcf().set_size_inches(10, 5) # Используйте (20, 10) для получения большего графика
plt.plot(fft.frequencies, fft.amplitudes)
plt.xlabel("Частота")
plt.ylabel("Амплитуда")
plt.savefig("/ram/fft-frequency-range.svg")

Частотное окно

Извлечение амплитуды и угла на определённой частоте

Используя [frequency], т.е. оператор getitem с одним значением, вы получаете объект FFTPoint(), содержащий частоту, амплитуду и относительный угол для данной частоты. Библиотека автоматически выбирает ближайший FFT-бакет, поэтому, даже если ваш FFT не имеет бакета для этой конкретной частоты, вы получите разумные результаты.

from UliEngineering.SignalProcessing.FFT import compute_fft
fft = compute_fft(data, samplerate=1e3)

# Показать значение на определённой частоте
print(fft[30]) # FFTPoint(frequency=30.0, value=2.91e-08, angle=0.0)
print(fft[100]) # FFTPoint(frequency=100.0, value=0.419, angle=0.0)

Короткие FFT для длинных данных

Используя compute_fft(), если у нас есть крайне длинный массив данных, это означает, что мы вычислим крайне длинный FFT. Во многих случаях это нежелательно, и вы хотите вычислить FFT фиксированного размера (чаще всего FFT степени двойки, например, 1024, 2048, 4096 и т.д.).

Давайте сгенерируем некоторые длинные тестовые данные и предположим, что мы хотим вычислить FFT размера 1024 на них

from UliEngineering.SignalProcessing.FFT import *
fft = simple_serial_fft_reduce(data, samplerate=1e3, fftsize=1024)
# fft.frequencies, fft.amplitudes и т.д.

Построение графика данных, как указано выше, даёт

что выглядит почти точно как наш график compute_fft() ранее — как мы и ожидали.

simple_serial_fft_reduce() заботится обо всей магии и нормализации для нас, включая разделение данных на перекрывающиеся чанки, добавление FFT и правильную нормализацию результата.

Соглашение об именовании здесь значимо:

• simple_...._reduce означает, что это вариант с разумными значениями по умолчанию (simple) для использования функции reduction (по умолчанию: sum) на нескольких FFT.
• serial означает отдельные FFT

Параллелизация FFT

Если у вас огромный набор данных, вы можете использовать simple_parallel_fft_reduce() идентично simple_serial_fft_reduce():

from UliEngineering.SignalProcessing.FFT import *
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024)
# Использовать fft.frequencies, fft.amplitudes и т.д.

Однако в большинстве случаев вы хотите инициализировать исполнитель вручную, чтобы вы могли повторно использовать его позже:

from UliEngineering.SignalProcessing.FFT import *
from concurrent.futures import ThreadPoolExecutor
executor = ThreadPoolExecutor() # Без аргумента => использовать num_cpus потоков
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024, executor=executor)

Мы можем использовать ThreadPoolExecutor(), поскольку scipy.fftpack (который UliEngineering использует для выполнения тяжёлых вычислений) разблокирует Python GIL.

Обратите внимание, что из-за необходимости выполнять много задач по ведению хозяйства, simple_parallel_fft_reduce() намного медленнее, чем simple_serial_fft_reduce(), если у вас набор данных настолько мал, что параллелизация неэффективна. Моя первоначальная рекомендация — рассмотреть использование параллельного варианта, если общее время выполнения последовательного варианта больше $0.5s$