Легке обчислення та візуалізація FFT у Python за допомогою UliEngineering

UliEngineering — це бібліотека аналізу даних змішаного типу в Python. Одна з утиліт, які вона надає, — це простий у використанні пакет для обчислення FFT. На відміну від інших пакетів, ця бібліотека орієнтована на практичне застосування і дозволяє виконати FFT лише одним рядком коду! Не потрібні знання математики. Спочатку встановіть UliEngineering.

Початок роботи

Спочатку ми згенеруємо деякі тестові дані. Див. цей попередній пост для детальнішої інформації про те, як згенерувати синусоїдальні тестові дані:

generate_test_data.py

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 Гц плюс синусоїди 400 Гц (з половинною амплітудою). Цей сигнал дискретизовано з частотою 1000 Гц.

Тепер ми можемо обчислити та візуалізувати FFT за допомогою matplotlib:

compute_fft.py

# Обчислюємо 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, amplitudes та phase. Вона виконує FFT розміру вхідних даних (тобто, оскільки data — масив довжини 1000, FFT буде розміру 1000).

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

Як ви можете побачити на графіку вище, максимальна частота, яку можна виявити, завжди становить половину частоти дискретизації, тобто для нашої частоти дискретизації $f_s = 1000,\text{Hz}$ вона дорівнює $500,\text{Hz}$. Див. цей більш математичний опис FFT, якщо ви хочете дізнатися більше деталей.

Внутрішньо compute_fft() виконує таке обчислення:

$$2 \cdot \frac{\text{abs}\left(\text{FFT}(\text{data} \cdot \text{Window})\right)}{\text{len(data)}}$$

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

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

Використовуючи API UliEngineering, вибір діапазону частот FFT — дуже просто: просто використовуйте fft[lowfreq:highfreq]. Ви можете використовувати fft[lowfreq:], щоб вибрати все починаючи з lowfreq, або fft[:highfreq], щоб вибрати все до highfreq.

fft_frequency_range.py

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 не має комірки для цієї конкретної частоти, ви отримаєте розумні результати.

fft_point.py

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 на них

simple_serial_fft.py

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

Побудова графіка даних як вище дає

Графік FFT синусоїд 100 Гц та 400 Гц з використанням послідовного FFT reduce

який виглядає майже точно як наш графік compute_fft() раніше — саме так, як ми і очікували.

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

Угода про найменування тут має значення:

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

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

Якщо ви маєте великий набір даних, ви можете використовувати simple_parallel_fft_reduce() так само, як simple_serial_fft_reduce():

simple_parallel_fft.py

from UliEngineering.SignalProcessing.FFT import *
fft = simple_parallel_fft_reduce(data, samplerate=1e3, fftsize=1024)
# Використовуйте fft.frequencies, fft.amplitudes тощо

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

simple_parallel_fft_executor.py

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$