Преобразования Фурье для обработки сигналов с помощью Python

Данный материал представляет собой незначительно сокращенный перевод статьи Кэмерона МакЛаудаFourier Transforms With scipy.fft: Python Signal Processing.

***

Преобразование Фурье– повсеместно используется для анализа сигналов – от обработки звука до сжатия изображений, от инженерных расчетов до Data Science. Популярная Python-библиотека SciPy предоставляет готовую реализацию преобразования Фурье в модулеscipy.fft.

Сам модуль поначалу выглядит устрашающе. Виной тому множество однотипно названных функций и то, чтодокументацияоперирует технической терминологией без каких-либо пояснений. Но есть и хорошие новости: чтобы начать использовать модуль, достаточно усвоить лишь несколько основных концепций.

Не волнуйтесь, если не чувствуете себя уверенно в математических изысканиях – мы рассмотрим алгоритм на конкретном практическом примере. Чтобы представить преобразование Фурье визуально, обязательно посмотрите следующийролик 3Blue1Brown(есть русские субтитры).

Из этого руководства вы узнаете:

• как и когда использовать преобразование Фурье;
• как для вашей задачи выбрать правильную функцию изscipy.fft;
• как посмотреть и изменить частотный спектр сигнала;
• примеры доступных вscipy.fftпреобразований.

Модуль scipy.fft

Установка SciPy и Matplotlib

Прежде чем начать, необходимо установить SciPy, NumPy (библиотека для работы с массивами) и Matplotlib (библиотека для визуализации данных). Вы можете сделать это одним из двух способов:

1. С помощью Anaconda: загрузите и установитеAnaconda Individual Edition. В этот набор инструментов уже включены перечисленные библиотеки.
2. С помощьюpipвы можете установить (или обновить) библиотеки посредством следующей команды:

         python -m pip install -U numpy scipy matplotlib     

Вы можете убедиться, что установка прошла успешно, запустив следующий код:

         import numpy, scipy, matplotlib print(numpy.__version__) print(scipy.__version__) print(matplotlib.__version__)     

Этот код импортирует NumPy, SciPy, Matplotlib и выведет версии модулей, если они установлены в системе.

Разница между scipy.fft и другими модулями

Очертим различия между модулемscipy.fftи другими модулями со схожими названиями. При просмотре документации SciPy вы обнаружите два родственных модуля:

• scipy.fft
• scipy.fftpack

Модульscipy.fftновее и предпочтительнее, чемscipy.fftpack:

• scipy.fftимеет улучшенный интерфейс;
• scipy.fftпозволяет использовать несколько воркеров, что в некоторых ситуациях может повысить скорость.

Реализация быстрого преобразования Фурье (англ. Fast Fourier transform, FFT) в SciPy содержит больше функций и более вероятно будет исправлена в случае обнаружения ошибки, чем реализация NumPy (numpy.fft), которая поддерживает реализацию FFT лишь для обеспеченияобратной совместимости.

Преобразование Фурье

Анализ Фурье – область математического анализа, отвечающая на вопрос, как можно представить математическую функцию в виде комбинации простых тригонометрических функций. Преобразование Фурье – это инструмент, который позволяет в исследуемом сигнале увидеть вклад каждой из этих гармонических составляющих, характеризуемых определенной частотой. В этом смысле говорят, что преобразование Фурье позволяет разложить функцию по частотам.

Итак, мы будем говорить о трех компонентах:

1. Сигнал– некоторая информация, которая меняется со временем. Например, аудиосигнал, видеосигнал, изменение разности электрических потенциалов – всё это примеры сигналов.
2. Частота(англ. frequency) – это скорость, с которой что-то повторяется. Например, часы тикают с частотой один герц (1 Гц) или, иначе говоря, совершают одно колебание в секунду.
3. Мощность(англ. power) – в данном случае просто мощность сигнала для каждой частоты.

Следующее изображение иллюстрирует примеры гармонических сигналов различной частоты и мощности.

Представьте, что вы использовали преобразование Фурье для записи того, как кто-то играет на фортепиано аккорд из трёх нот.

Результирующий частотный спектр покажет три пика – по одному для каждой ноты. Если человек играл одну ноту мягче, мощность для частоты этой ноты будет меньше, чем для двух других.

Зачем может понадобиться преобразование Фурье?

Преобразование Фурье полезно во многих приложениях. Например, Shazam и другие службы распознавания музыки используют преобразование Фурье для идентификации песен. Алгоритм сжатияJPEGпредставляет собой вариант преобразования Фурье, применяемый для удаления высокочастотных компонент изображений. В распознавании речи преобразование Фурье и связанные с ним преобразования служат для восстановления произнесенных слов.

Задача преобразования Фурье возникает всякий раз, когда нужно как-либо работать с сигналом, представляемым в пространстве частот.

Временная область против частотной области

Далее мы будем иметь дело свременно?йи частотной областями] – двумя подходами к представлению сигнала: как информации, которая изменяется во времени и информации, отображенной в виде набора частот и соответствующих им амплитуд.

Ниже представлено характерное изображение аудиосигнала – классического примера сигнала во временной области. Горизонтальная ось соответствует времени, вертикальная ось – амплитуде.

Тот же звуковой сигнал можно представить разложенным по составляющим его частотам. Горизонтальная ось на рисунке ниже представляет частоту, вертикальная ось – мощность.

Классификация преобразований Фурье

Преобразование Фурье подразделяют на категории по нескольким признакам. В первую очередь – по типу функций, с которыми работает преобразование: непрерывные или дискретные. В этом руководстве мы рассматриваемдискретное преобразование Фурье(DFT).

Термины DFT и FFT нередко используются как взаимозаменяемые. Однако это не совсем одно и то же: быстрое преобразование Фурье (FFT) – лишь один из алгоритмов вычисления дискретного преобразования Фурье.

Еще одна линия раздела в терминологии, с которым вы столкнетесь при использованииscipy.fft,– разные типы ввода. Например, функцияfft()принимает комплексные числа, аrfft()работает только с действительными числами. В дальнейшем мы обсудим это подробнее.

Практический пример: удаление нежелательного шума из аудиофайла

Чтобы лучше понять преобразование Фурье и то, как его можно применить, решим задачу фильтрации звука. Намеренно создадим звуковой сигнал с высокочастотным шумом, а затем удалим шум с помощью преобразования Фурье.

Создание сигнала

Одиночное гармоническое (синусоидальное) колебание представляют одну частоту и в музыкальном отношении является чистым тоном. Воспользуемся свойством таких волн для генерации звука:

         import numpy as np from matplotlib import pyplot as plt  SAMPLE_RATE = 44100  # Гц DURATION = 5  # Секунды  def generate_sine_wave(freq, sample_rate, duration):     x = np.linspace(0, duration, sample_rate*duration, endpoint=False)     frequencies = x * freq     # 2pi для преобразования в радианы     y = np.sin((2 * np.pi) * frequencies)     return x, y  # Генерируем волну с частотой 2 Гц, которая длится 5 секунд x, y = generate_sine_wave(2, SAMPLE_RATE, DURATION) plt.plot(x, y) plt.show()     

После импорта NumPy и Matplotlib мы определили две константы:

1. SAMPLE_RATE(частота дискретизации) определяет, сколько точек используется для представления синусоидальной волны на интервале 1 с. Если бы сигнал имел частоту дискретизации 10 Гц и представлял пятисекундную синусоидальную волну, то он содержал бы 50 точек данных.
2. DURATION– длина сгенерированной выборки.

Затем мы определяем функцию для генерации синусоидальной волны – позже мы воспользуемся ей несколько раз. Функция принимает частотуfreqи возвращает значенияxиy, которые далее будут использоваться для построения изображения сигнала.

Координатыxсинусоидальной волны равномерно распределены между 0 иDURATION. Установкаendpoint = Falseв функцииnp.linspace()важна для правильной работы преобразования Фурье – предполагается, что сигнал является периодическим.

Осьxпредставляет время в секундах – обратите внимание, что синусоидальная волна действительно совершает два колебания в секунду. Эта синусоида имеетслишком низкую частоту, чтобы ее можно было слышать, поэтому в следующем разделе мы сгенерируем несколько высокочастотных синусоид и рассмотрим, как их можно смешивать.

Смешивание аудиосигналов

Микширование аудиосигналов состоит всего из двух этапов:

• cложение сигналов;
• нормализация результата.

         _, nice_tone = generate_sine_wave(400, SAMPLE_RATE, DURATION) _, noise_tone = generate_sine_wave(4000, SAMPLE_RATE, DURATION)  noise_tone = noise_tone * 0.3 mixed_tone = nice_tone + noise_tone     

Символ подчеркивания (_) мы используем, чтобы отбросить значенияx, возвращаемые функциейgenerate_sine_wave()– нам не нужно складывать значения времени.

Следующий шаг – нормализация, масштабирование сигнала под целевой формат. В нашем случае это 16-битное целое число в диапазоне от-32768до32767:

         normalized_tone = np.int16((mixed_tone / mixed_tone.max()) * 32767)  plt.plot(normalized_tone[:1000]) plt.show()     

Делениеmixed_toneна максимальное значение масштабирует его в интервале от-1до1. Умножение на32767масштабирует сигнал между-32767и32767, что примерно соответствует диапазонуnp.int16. Код отображает только первые1000точек, чтобы мы могли четче проследить структуру сигнала. Видимая нами синусоидальная волна – это сгенерированный тон 400 Гц, искаженный тоном 4000 Гц.

Чтобы прослушать звук, необходимо сохранить его в формате, который может прочитать аудиоплеер. Воспользуемся методом SciPywavfile.writeи сохраним результат в файле форматаWAV. Выбранное нами 16-битное целочисленное представление является стандартным типом данных для wav-файлов.

         from scipy.io.wavfile import write  write("mysinewave.wav", SAMPLE_RATE, normalized_tone)     

Этот код запишет данные в файлmysinewave.wavв директории, где мы запускаем этот скрипт Python. Файл можно прослушать с помощью любого медиаплеера.

Используем быстрое преобразование Фурье для удаления шума

Быстрое преобразование Фурье (FFT) – алгоритм, который позволяет вычислить частотный спектр сигнала:

         from scipy.fft import fft, fftfreq  # число точек в normalized_tone N = SAMPLE_RATE * DURATION  yf = fft(normalized_tone) xf = fftfreq(N, 1 / SAMPLE_RATE)  plt.plot(xf, np.abs(yf)) plt.show()     

На построенном спектре видны два пика на положительных частотах и два их зеркальных отражения в отрицательной области. Пики положительных частот находятся на позициях 400 и 4000 Гц.

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

Самый важный раздел в этом небольшом скрипте – вычисление преобразования Фурье:

         yf = fft(normalized_tone) xf = fftfreq(N, 1/SAMPLE_RATE)     

Код вызывает две функции:

1. fft()вычисляет само преобразование.
2. fftfreq()находит частоты в центре каждого «бина» на выходеfft(). Без этого не было бы возможности построить осьxнашего спектра.

Под бином здесь понимается интервал значений, сгруппированных аналогично гистограмме. В рамках этого руководства достаточно рассматривать их как отдельные значения.

Интересной частью кода является обработка, выполняемая сyfперед построением – вызовnp.abs()для yf вызван лишь тем, что значенияyf– комплексные числа.

Комплексное число– это число, состоящее из двух частей: действительной и мнимой. Такие необычные числа полезны во многих приложениях, но если вы столкнулись с ними впервые, то сейчас достаточно знать лишь то, что они существуют.

Математики обычно записывают комплексные числа в формеa + bi, гдеa– действительная часть,b– мнимая часть,i–мнимая единица.

Поскольку комплексные числа состоят из двух компонент, построение графика их зависимости от частоты на двумерной оси требует, преобразовать два значения в одно. На помощь приходитnp.abs(). Эта функция вычисляет ?(a?+b?).

Обрабатываем сигнал еще быстрее с помощью rfft()

Частотный спектр, выдаваемыйfft(), зеркально отражался относительно оси y. Эта симметрия вызвана вводом в преобразование действительных чисел. Эту симметрию можно использовать, чтобы ускорить преобразование Фурье, вычислив лишь половину с помощью функцииrfft().

         from scipy.fft import rfft, rfftfreq  # обратите внимание на r в начале имён функций yf = rfft(normalized_tone) xf = rfftfreq(N, 1/SAMPLE_RATE)  plt.plot(xf, np.abs(yf)) plt.show()     

Фильтрация сигнала

Самая замечательная вещь в преобразовании Фурье заключается в том, что оно обратимо. Любой сигнал, измененный в частотной области, можно преобразовать обратно во временную область. Воспользуемся этим, чтобы отфильтровать высокочастотный шум.

Возвращаемыеrfft()значения соответствуют мощности каждого частотного бина. Если мы установим мощность бина равной нулю, соответствующая частота перестанет присутствовать в результирующем сигнале во временной области:

         # Максимальная частота составляет половину частоты дискретизации points_per_freq = len(xf) / (SAMPLE_RATE / 2)  # Наша целевая частота - 4000 Гц target_idx = int(points_per_freq * 4000)     

Обнулим yf для индексов около целевой частоты:

         yf[target_idx-2:target_idx+2] = 0  plt.plot(xf, np.abs(yf)) plt.show()     

Остался только один пик. Применим обратное преобразование Фурье, чтобы вернуться во временную область.

Применение обратного преобразования Фурье

Применение обратного FFT аналогично применению FFT:

         from scipy.fft import irfft  new_sig = irfft(yf)  plt.plot(new_sig[:1000]) plt.show()     

Поскольку мы использовалиrfft(), для обратного преобразования нужно использоватьirfft(). Однако, если бы мы использовалиfft(), обратной функцией была быifft().

Как видите, теперь есть одна синусоида, колеблющаяся с частотой 400 Гц – мы успешно удалили шум на 4000 Гц.

Нормализуем сигнал и запишем результат в файл. Сделать это можно так же, как в прошлый раз:

         norm_new_sig = np.int16(new_sig * (32767 / new_sig.max())) write("clean.wav", SAMPLE_RATE, norm_new_sig)     

Проиграв файл, вы услышите, что раздражающий писк исчез.

Дискретные косинусные и синусоидальные преобразования

Туториал по модулюscipy.fftбыл бы неполным без рассмотрения дискретногокосинусного(DCT) и синусоидального (DST) преобразований. Эти два преобразования тесно связаны с преобразованием Фурье, но работают только с действительными числами. В библиотеке SciPy соответствующие преобразования реализованы в виде функцийdct()иdst(). Варианты этих функций с названиями, начинающимися сiиn, представляют соответственно обратные и n-мерные версии функций.

Упрощенно говоря, DCT и DST – как бы две половины преобразования Фурье, вычисляемые по отдельности быстрее, чем полное преобразование Фурье. Прежде чем вы научитесь выбирать между ними, нужно освежить в памяти понятиечетныхинечетныхфункций. Четные функции симметричны относительно осиy, а нечетные – относительно начала координат. Чтобы представить это наглядно, взгляните на следующие примеры.

При расчете полного преобразования Фурье (DFT) предполагается, что функция, по которой происходит вычисление, повторяется бесконечно. Однако преобразования DCT и DST позволяют учесть симметрию сигнала. Косинусное преобразование (DCT) предполагает, что функция продлевается за счет четной симметрии, а для DST – за счет нечетной симметрии.

На следующем изображении показано, как каждое преобразование представляет, как функция будет продолжаться в бесконечности.

На изображении выше полное преобразование повторяет функцию как есть. DCT отражает функцию по вертикали, а DST – по горизонтали. Обратите внимание, что симметрия DST приводит к существенным разрывам функции. Это вносит высокочастотные составляющие в результирующем частотном спектре. Если нет сведений о симметрии сигнала, лучше использовать DCT.

Естьмножество примеровиспользования DCT в различных задачах, требующих высокой скорости преобразования Фурье, в том числе в алгоритмах JPEG, MP3 и WebM.

Заключение

Преобразование Фурье – это мощная концепция, применяемая в самых разных областях – от чистой математики до аудиотехники и даже финансов. В этом уроке мы рассмотрели:

• как и когда используется преобразование Фурье
• как выбрать нужную функцию изscipy.fft
• в чем разница между временной и частотной областями
• как посмотреть и изменить частотный спектр сигнала
• как использоватьrfft(), чтобы преобразование выполнялось еще быстрее

Мы рассмотрели только базовую идею, но ее понимание поможет разобраться в других вопросах, связанных с преобразованием Фурье и представлением функций в виде частотных спектров.