Как с помощью fft ускорить звук

Практическое введение в анализ частотного диапазона

В этом примере показано, как выполнить и интерпретировать основной анализ сигнала частотного диапазона. Пример обсуждает преимущества использования частотного диапазона по сравнению с представлениями временного интервала сигнала и иллюстрирует фундаментальные понятия с помощью симулированных и действительных данных. Пример отвечает на основные вопросы, такие как: каковы значение величины и фаза БПФ? Действительно ли мой сигнал является периодическим? Как я измеряю степень? Есть ли один, или больше чем один сигнал в этой полосе?

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

Сигнал может быть преобразован между временным и частотным диапазоном с парой математических операторов, названных преобразованием. Примером является преобразование Фурье, которое разлагает функцию на сумму (потенциально бесконечного) количества частотных составляющих синусоиды. ‘Спектр’ частотных составляющих является представлением частотного диапазона сигнала. Обратное преобразование Фурье преобразует функцию частотного диапазона назад в функцию времени. fft и ifft функции в MATLAB позволяют вам вычислять Дискретное преобразование Фурье (ДПФ) сигнала и инверсию этого преобразования соответственно.

Величина и информация о фазе БПФ

Представление частотного диапазона сигнала несет информацию о величине сигнала и фазе на каждой частоте. Поэтому выход расчета БПФ является комплексным. Комплексное число , имеет действительную часть $x_r$ и мнимую часть $x_i$ , такую что $x = x_r + jx_i$ . Величина $x$ вычисляется как $$\sqrt<(x_r^2+x_i^2)>$» /> , и фаза <img decoding=$ вычисляется как $$\arctan<(x_i/x_r)>$» /> . Можно использовать функции MATLAB abs и angle соответственно получить величину и фазу любого комплексного числа.</p> <p>Используйте аудио пример, чтобы разработать некоторое понимание на том, какую информацию несут величина и фаза сигнала. Для этого загрузите звуковой файл, содержащий 15 секунд музыки акустической гитары. Частота дискретизации звукового сигнала составляет 44,1 кГц.</p> <pre>Fs = 44100; y = audioread(<span style=$ 'guitartune.wav');

Используйте fft наблюдать содержимое частоты сигнала.

NFFT = length(y); Y = fft(y,NFFT); F = ((0:1/NFFT:1-1/NFFT)*Fs).';

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

Постройте величину и компоненты фазы спектра частоты сигнала. Величина удобно построена в логарифмическом масштабе (дБ). Фаза развернута с помощью unwrap функционируйте так, чтобы мы видели непрерывную функцию частоты.

magnitudeY = abs(Y); % Magnitude of the FFT phaseY = unwrap(angle(Y)); % Phase of the FFT helperFrequencyAnalysisPlot1(F,magnitudeY,phaseY,NFFT)

Можно применять обратное преобразование Фурье к вектору частотного диапазона, Y, чтобы восстановить сигнал времени. ‘Симметричный’ флаг говорит ifft то, что вы имеете дело с сигналом времени с действительным знаком, таким образом, он обнулит маленькие мнимые компоненты, которые появляются на обратном преобразовании из-за числовых погрешностей в расчетах. Заметьте, что исходный сигнал времени, y, и восстановленный сигнал, y1, является практически тем же самым (норма их различия находится порядка 1e-14). Очень небольшая разница между этими двумя происходит также из-за числовых упомянутых выше погрешностей. Вопроизведите и послушайте непреобразованный сигнал y1.

y1 = ifft(Y,NFFT,'symmetric'); norm(y-y1)

БПФ. Как получить точность в 1Гц быстрее чем за секунду?

Добрый день.
Заранее прошу прощения, если буду плавать в сабже — в универе не преподавали преобразование Фурье.
В общем, я пытаюсь написать гитарный тюнер на verilog. Имеется аудио-кодек, который оцифровывает звук с частотами дискретизации 8-96кГц. Есть блок БПФ, в котором можно настроить сколько точек он будет принимать. Мне нужно получить точность в 1Гц (или лучше) в диапазоне частот 60-700Гц. При этом желательно иметь минимальную задержку. Насколько я понимаю, точность зависит от частоты дискретизации и числа точек в БПФ: error = fs/num_points. Получается, чтобы иметь погрешность Можно ли ускорить этот процесс? И да, мне известно про принцип неопределенности.

Отслеживать

20.2k 6 6 золотых знаков 37 37 серебряных знаков 81 81 бронзовый знак

задан 25 фев 2018 в 11:10

Андрей Солодовников Андрей Солодовников

718 5 5 серебряных знаков 23 23 бронзовых знака

Минутка «взрослого» программирования: как растягивать звук и как менять тональность.

Ни для кого не секрет, что если воспроизводить звук ускоренно, то изменится и тон его звучания и, разумеется, уменьшится длительность.

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

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

Чтобы было понятнее о чём пойдёт речь далее, небольшая видеозарисовка:

Ясно дело, без хитростей тут не обойдётся. На помощь нам приходит штука, известная как «Быстрое преобразование Фурье» (БПФ, или по-буржуйски FFT: Fast Fourier Transform).
О том как реализовать БПФ, я тут расписывать не буду, благо весь интернет кишит статьями с исходниками. Важно лишь то, что БПФ работает над массивом комплексных чисел и позволяет получить спектр сигнала.

Вот тут я накидал поясняющую картинку:

Я думаю по картинке всё и так понятно, но на всякий случай напишу пояснения:

1. Исходный сигнал разбивается на фрагменты, причём фрагменты перекрываются большей своей частью. Важно подойти к выбору длины каждого фрагмента: если фрагмент будет слишком коротким, то преобразование не затронет низкие частоты звука, не попавшие во фрагмент, если слишком длинным, то появится неприятный трубный звук, из-за того что ухо будет успевать различать изменения внутри фрагментов. Идеальная длительность фрагмента — около 50 мс. Но применение БПФ сильно ограничивает нас выбором длительности фрагментов кратным степени двойки. По моему опыту, для частоты воспроизведения 44100 семплов в секунду наилучший результат получается с фрагментами длиной 2048 сэмплов.

Как говорилось, фрагменты перекрывают друг друга, это значит что каждый последующий цикл преобразования захватывает большую часть исходных данных, уже участвовавших в предыдущих итерациях. Степень перекрытия фрагментов сильно влияет на качество звучания. По моему опыту, хороший результат получается при 16-кратном перекрытии (такое используется в демонстрационном видео), это означает что начала фрагмента для очередного цикла преобразования отстоит на 1/16 длины фрагмента от начала фрагмента прошлого цикла. При длине фрагмента 2048 сэмплов этот шаг составляет 128 сэмплов.

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

2. Итак, выбран один из фрагментов исходного звука, с него и начинается цикл преобразований.

3. Поскольку наибольшие ошибки БПФ накапливаются по краям окна, сигнал нужно отфильтровать, уменьшив его амплитуду на краях путём умножения на функцию оконного фильтра.

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

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

f(i) = sin((Π / 2) * sin(Π * (i / N))²)

где Π — это пи (3,14159…); N — длина окна, i — номер сэмпла в окне (от 0 до N — 1)

4. БПФ работает над комплексными числами, а исходный сигнал у нас вещественный. Для подготовки массива для преобразования исходный сэмпл просто копируется в вещественную часть комплексного числа, а мнимая часть остаётся нулевой.

Теперь можно выполнять преобразование.

5. Поскольку все мнимые части у нас были нулевыми, полученный спектр получился симметричным. Только мнимая часть в правой половине имеет противоположный знак. Поэтому правую часть спектра можно смело отбросить.

В общем и целом, вещественная часть i-го элемента массива (гармоники) означает амплитуду косинуса частоты 2 * Π * i / N, а мнимая — синуса.

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

6. Для каждой гармоники спектра комплексное число необходимо преобразовать в полярные координаты. Теперь получилась l — амплитуда соответствующей частотной составляющей, и α — смещение фазы, т.е. насколько верхняя точка синусоиды отстоит от начала окна.

7. Представим такую ситуацию: у нас есть длинный протяжный звук, к примеру синусоида. Если частота синусоиды совпадает с частотой одной из гармоник, то на спектре будет присутствовать только эта гармоника, а вот фаза α будет меняться на каждой итерации, в зависимости от того на какую часть синусоиды пришлось начало окна. Не трудно догадаться что разница между фазами в двух соседних итерациях будет постоянной.

Читать:

Какие ключевые слова используются в javascript для того чтобы объявить переменную

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

Представим теперь, что амплитуда гудящей синусоиды нарастает со временем. Теперь растягивание её во времени будет означать интерполяцию значений амплитуд (l), но разница между фазами в двух соседних итерациях опять таки будет изменяться с прежней скоростью.

Собственно эта идея и лежит в основе алгоритма:
1) для каждой гармоники на каждой итерации смотрим скорость изменения фазы — т.е. разницу с фазой из предыдущей итерации.
2) полученное значение скорости прибавляется к значению фазы такой же гармоники из предыдущей итерации выходного спектра, вне зависимости от того ускоряем мы, или замедляем воспроизведение.
3) значение амплитуды гармоники в выходном спектре получается путём интерполяции.

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

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

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

10. Затем выходной сигнал нужно домножить на оконный фильтр. Здесь я использую ту же формулу, что и в п.3, с последующим делением на половину степени перекрытия.

11. Полученный сигнал складывается с учётом степени перекрытия, с результатами предыдущих и последующих итераций

12. Выходной сигнал готов, можно его выводить на аудиоустройство, или сохранять в файл — кому что ближе.

Изменение тона

Для изменения тона в N раз выполняется растяжение (т.е. замедление) звука в эти же N раз, а затем ускорение воспроизведения в N раз, что приводит к изменению тона при сохранении исходного темпа воспроизведения.

Модуль scipy.fft

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

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

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

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

Схематическое представление аккорда и соответствующего ему частотного спектра

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

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

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

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

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

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

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

Аудиосигнал во временной области

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

Тот же аудиосигнал в частотной области

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

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

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

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

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

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

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

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

Вид смикшированного сигнала

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

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

Результат FFT-преобразования

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

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

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

Форма спектра сигнала до фильтрации

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

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

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

Форма спектра сигнала после фильтрации

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

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

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

Форма сигнала после фильтрации

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

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

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

Примеры четной и нечетной функций – соответственно квадратичная и кубическая функции

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

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

Представление конечного дискретного сигнала в случае полного, косинусного и синусоидального преобразований Фурье

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

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

Заключение

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

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

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

Больше полезной информации вы можете получить на нашем телеграм-канале «Библиотека питониста». Рекомендуем также обратить внимание на учебный курс по Python от «Библиотеки программиста».