Ограниченность преобразования Фурье или почему стоит доверять своему слуху
Последние несколько десятков лет задача распознавания аккордов музыкальной композиции ставилась довольно часто. Казалось бы, этот не столь оригинальный сервис был и остается довольно распространенным среди приложений, работающих со звуком (Ableton, Guitar Pro), однако универсального, срабатывающего всегда алгоритма не существует до сих пор. В этой статье я постараюсь раскрыть одну из множества причин неидеальной работы подобных сервисов на примере алгоритмов, использующих в своей основе преобразование Фурье.
Большинство аудиоформатов хранит информацию в виде зависимости амплитуды от времени (например, .wav) или в виде коэффициентов частотного преобразования (.mp3, .aac, .ogg), однако современные сложные алгоритмы цифровой обработки сигналов работают с частотной составляющей звуков. Приходится иметь дело с двойным переходом, из амплитудной области в частотную, затем обратно. На данный момент осуществление такого переход без потерь в качестве является достаточно распространенной проблемой с множеством неоднозначных решений.
Для разложения сигнала на частоты используется распространенное во многих областях преобразование Фурье. Поскольку в большинстве алгоритмах цифровой обработки сигналов приходится иметь дело с некоторой выборкой данных об амплитуде сигнала в конкретный момент времени, удобно использовать одну из многочисленных модификаций преобразования Фурье — дискретное преобразование Фурье. В своей классической интерпретации преобразование Фурье представляет собой сумму вида:
где N — количество компонент разложения,
xn, n = 0…N — 1 — значение отсчета сигнала в момент времени с номером n,
Xk, k = 0. N — 1 — выходные значения преобразования, представляющие собой искомые значения частот,
k — индекс частот.
Полученные в результате прямого преобразования значения Xk представляют собой набор комплексных значений — пар Re(Xk) + i*Im(Xk), характеризующих амплитуду и начальную фазу гармонического сигнала. Применив к полученным значениям обратное преобразование Фурье
мы восстановим исходный сигнал с некоторой точностью. На практике же часто приходится работать с большими объемами данных. Для значительного ускорения вычисления используют быстрое преобразование Фурье, во многих современных библиотеках именуемое FFT (Fast Fourier Transform). FFT также работает с комплексными числами и отличается тем, что размер самого преобразования обязательно является степенью двойки, на практике это чаще всего мы видим числа 1024 или 2048. Варьируем размер преобразования и величину выборки, получаем приближенный к значениям реального сигнала набор отсчетов, воздействуем на него обратным преобразованием Фурье, и на выходе перед нами исходные частоты, а значит, получение звучавшей последовательности аккордов не должно представлять особой сложности.
Однако. Преобразование Фурье любого вида является сложным универсальным инструментом, сильно зависящим от поставляемых ему данных. Непрерывное преобразование Фурье является развитием идеи рядов Фурье, более обобщенной, результирующим значением частот которого являются комплексные числа. Мы смотрим на классическое прямое преобразование
и видим в его реализации интеграл на бесконечном промежутке, на дискретное преобразование — это бесконечная сумма ряда. И о справедливости равенства правой и левой частей мы можем говорить только в том случае, когда исходная функция является функцией бесконечной длины, только тогда, когда перед нами представлен бесконечный сигнал. Легко согласиться, что на практике такие сигналы встретишь не часто.
Создадим искусственный сигнал, соответствующий ноте Ля, взятой в двух разных октавах (частоты равны соответственно 440 Гц и 3440 Гц), пропустим его через прямое преобразование Фурье, а затем через обратное:
dt = 0.001; T = 1; t=0:dt:T; w = 440; q = 2^(1/12); //исходный сигнал x = cos(2*%pi*w*t*q^1) + cos(2*%pi*w*t*q^4); scf(0); plot2d(x, rect=[0,0,1200,2]); // сигнал после преобразования Фурье y = fft(x); scf(1); plot(y); // обратное преобразование Фурье z = fft(y,-1); z1 = z/(T/dt); scf(2); plot2d(x, style=[color("red")], rect=[0,0,1200,2]); plot2d(z1,style=[color("green")], rect=[0,0,1200,2]); Запускаем код Scilab и видим, что исходный сигнал (x) и сигнал, полученный обратным преобразованием (z1), мягко говоря, разные:
Рис.1 График функции x исходного сигнала (красный цвет) и восстановленного сигнала z1 (желтый цвет)
Как обходят подобные ограничения с бесконечностью? Поскольку перевод сигнала в частотное представление при помощи FFT производится только блоками, у каждого такого блока мы можем увеличить частоту дискретизации, тем самым увеличив размер FFT блока и, соответственно, количество коэффициентов преобразования. Это даст нам более точный частотный результат на единицу времени. Достаточно разумный подход.
Но как бы мы не приближали длину промежутка к бесконечности, неточности все равно будут иметь место. FFT не знает ничего об исходных гармониках функции. В нашем случае, имелась частота в 440 Гц, а FFT превратил ее в набор гармоник с искаженными частотами, хоть и близкими по значению к исходной. Поэтому полученные значения сигнала являются лишь близкими к первоначальному. Чтобы ликвидировать это искажение, применяются так называемые оконные функции. Перед разложением блок FFT перемножается на оконную функцию W(ω — t), где t — момент времени, соответствующий текущему отсчету. Основное влияние оконной функции заключается в уменьшении коэффициентов перед мнимыми частотами в разложении, тем самым уменьшая их вклад в конечную частотную картину.
Но как мы поймем, что мы отсекли действительно “незначительные” частотные пики? Конечно, можно пытаться экспериментировать с разбиениями и их длиной, длиной промежутков с нулевыми коэффициентами вне оконной функции, однако это лишь частично сгладит итоговую картину. Взять хотя бы случай, когда искомые ноты имеют достаточно низкие частоты, здесь мы можем получить мнимые значения частот, значительно превышающие исходный материал!
Здесь уже явно стоит попробовать эксперименты с оконными функциями разного вида. Например, классическая прямоугольная оконная функция
сохранит все исходные значения частот. Но если мы возьмем, скажем, треугольное окно
станет очевидно, что краевые выбивающиеся значения попросту отсекутся, с высокой долей вероятности останется наиболее чистая информация.
Однако. Снова сгенерируем простой ограниченный сигнал и вычислим его спектр:
scf(0); t=0:0.0001:3*%pi; x = sin(t); plot2d(t,x); y = simp(fft(x)); scf(1); plot2d(y, rect = [-100, -10000, 40000, 10000]); disp(y(1000)); Исходный сигнал х:
Рис.2 Исходный сигнал х
Рис.3 Спектр сигнала y
Что же мы видим? Спектр перестал быть ограниченным, по обоим его сторонам образовались лишние значения, более того, их бесконечное число! Случилось что? В тех двух точках, где произошло соприкосновение с краями оконной функции, функция перестала быть гладкой. На вход Фурье-преобразование ожидает увидеть сигнал, представляющий собой сумму простых синусоид, где точки разрыва можно получить только с помощью бесконечной суммы этих синусоид.
Вот мы уже напрямую столкнулись с самым настоящим принципом неопределенности для области обработки сигналов. Принцип здесь носит название Теоремы Бенедика: функция не может быть одновременно ограниченной в диапазоне времени и в диапазоне частоты. В более строгой формулировке теорема звучит так: функция в пространстве 
и ее отображение Фурье 
не могут одновременно иметь один и тот же носитель функции на покрытиях с ограниченной мерой Лебега. И если с сигналами бесконечной длины дело приходится иметь редко, и эта область для практики имеет меньшее значение, то что же со случаем, когда сигнал имеет конечную длину?
И здесь мы не обойдемся без самой известной теоремы в теории обработки звуковых сигналов — теоремы Котельникова-Найкваиста-Шеннона, которая в одной из свои интерпретаций звучит следующим образом: любой аналоговый сигнал может быть восстановлен с какой угодно точностью по своим дискретным отсчетам тогда, когда значения его отсчетов f > 2Fmax, где Fmax — максимальная частота, которой ограничен спектр исходного сигнала. Исходя из этой теоремы и теоремы Бенедика мы делаем вывод, что в данный момент реальных сигналов, которые можно с бесконечно большой точностью представить в дискретном и ограниченном виде, нет. Мы всегда либо будем иметь дело с сигналом бесконечно большой длины, либо будем получать в ограниченном сигнале неточности разной степени, которые в итоге могут значительно испортить картину исходной аудиозаписи.
Конечно, разработчики подобных сервисов, работающих со звуками, эту проблему в какой-то мере решают. Используют ступенчатую интерполяцию, разные виды линейной интерполяции, фильтры нижних частот для борьбы с отсчетами, превышающими половину частоты дискретизации. Их самыми распространенными недостатками являются данные, которые лишь грубо приближены к исходным значениям частот, а также наличие высокочастотного мусора. Однако и аудиозаписи, с которыми приходится иметь дело, тоже далеко не идеальны. Высокий уровень шума, скорость записи, качество звукозаписи. Здесь уже нужно уметь адаптироваться под каждый сигнал, вычислять скорость игры, длительность аккордов, гасить белый (в лучшем случае) шум, но это уже совсем другая история.
- преобразование фурье
- fft
- распознавание музыки
- алгоритмы
- принцип неопределённости
Документация
Преобразование Фурье является инструментом для выполнения частоты и анализа спектра мощности сигналов временной области.
Количества спектрального анализа
Спектральный анализ изучает спектр частоты, содержавшийся в дискретном, однородно выборочных данных. Преобразование Фурье является инструментом, который показывает частотные составляющие времени — или основанный на пробеле сигнал путем представления его в пространстве частоты. В следующей таблице перечислены общие количества, используемые, чтобы охарактеризовать и интерпретировать свойства сигнала. Чтобы узнать больше о преобразовании Фурье, смотрите преобразования Фурье.
Демонстрационная частота (выборки в единицу времени или пробел)
Шаг времени или пространства на выборку
Область значений времени или пространства для данных
Дискретное преобразование Фурье данных (ДПФ)
Частота Найквиста (средняя точка частотного диапазона)
Сигнал с шумом
Преобразование Фурье может вычислить частотные составляющие сигнала, который повреждается случайным шумом.
Создайте сигнал с частотами компонента на уровне 15 Гц и 40 Гц, и введите случайный Гауссов шум.
rng('default') fs = 100; % sample frequency (Hz) t = 0:1/fs:10-1/fs; % 10 second span time vector x = (1.3)*sin(2*pi*15*t) . % 15 Hz component + (1.7)*sin(2*pi*40*(t-2)) . % 40 Hz component + 2.5*randn(size(t)); % Gaussian noise;
Преобразование Фурье сигнала идентифицирует свои частотные составляющие. В MATLAB®, fft функция вычисляет преобразование Фурье с помощью алгоритма быстрого преобразования Фурье. Используйте fft вычислить дискретное преобразование Фурье сигнала.
y = fft(x);
Постройте спектр мощности в зависимости от частоты. В то время как шум маскирует частотные составляющие сигнала на основанном на времени пробеле, преобразование Фурье показывает их как скачки в степени.
n = length(x); % number of samples f = (0:n-1)*(fs/n); % frequency range power = abs(y).^2/n; % power of the DFT plot(f,power) xlabel('Frequency') ylabel('Power')
Во многих приложениях более удобно просмотреть спектр мощности, сосредоточенный на 0 частотах, потому что это лучше представляет периодичность сигнала. Используйте fftshift функция, чтобы выполнить циклический сдвиг на y , и постройте степень в центре 0.
y0 = fftshift(y); % shift y values f0 = (-n/2:n/2-1)*(fs/n); % 0-centered frequency range power0 = abs(y0).^2/n; % 0-centered power plot(f0,power0) xlabel('Frequency') ylabel('Power')
Звуковой сигнал
Можно использовать преобразование Фурье, чтобы анализировать спектр частоты аудиоданных.
Файл bluewhale.au содержит аудиоданные от Тихоокеанской вокализации голубого кита, зарегистрированной подводными микрофонами недалеко от берегов Калифорнии. Файл от библиотеки вокализаций животных, обеспеченных Программой исследований Биоакустики Корнелльского университета.
Поскольку вызовы голубого кита являются настолько низкими, они являются едва слышимыми людям. Масштаб времени в данных сжат на коэффициент 10, чтобы повысить тангаж и выполнить более явно слышимый вызов. Считайте и постройте аудиоданные. Можно использовать команду sound(x,fs) слушать аудио.
whaleFile = 'bluewhale.au'; [x,fs] = audioread(whaleFile); plot(x) xlabel('Sample Number') ylabel('Amplitude')
Первый звук является «трелью», сопровождаемой тремя «стонами». Этот пример анализирует один стон. Задайте новые данные, которые приблизительно состоят из первого стона, и откорректируйте данные времени с учетом factor-10 ускорения. Постройте усеченный сигнал в зависимости от времени.
moan = x(2.45e4:3.10e4); t = 10*(0:1/fs:(length(moan)-1)/fs); plot(t,moan) xlabel('Time (seconds)') ylabel('Amplitude') xlim([0 t(end)])
Преобразование Фурье данных идентифицирует частотные составляющие звукового сигнала. В некоторых приложениях, что процесс большие объемы данных с fft , распространено изменить размер входа так, чтобы количество отсчетов было степенью 2. Это может сделать расчет преобразования значительно быстрее, особенно для объемов выборки с большими простыми множителями. Задайте новую длину сигнала n это — степень 2, и используйте fft функция, чтобы вычислить дискретное преобразование Фурье сигнала. fft автоматически заполняет исходные данные нулями, чтобы увеличить объем выборки.
m = length(moan); % original sample length n = pow2(nextpow2(m)); % transform length y = fft(moan,n); % DFT of signal
Настройте частотный диапазон из-за фактора ускорения, и вычислите и постройте спектр мощности сигнала. График показывает, что стон состоит из основной частоты приблизительно 17 Гц и последовательности гармоник, где вторая гармоника подчеркнута.
f = (0:n-1)*(fs/n)/10; power = abs(y).^2/n; plot(f(1:floor(n/2)),power(1:floor(n/2))) xlabel('Frequency') ylabel('Power')
Смотрите также
Похожие темы
- Преобразования Фурье
- 2D преобразования Фурье
Открытый пример
У вас есть модифицированная версия этого примера. Вы хотите открыть этот пример со своими редактированиями?
Документация MATLAB
Поддержка
- MATLAB Answers
- Помощь в установке
- Отчеты об ошибках
- Требования к продукту
- Загрузка программного обеспечения
© 1994-2021 The MathWorks, Inc.
- Условия использования
- Патенты
- Торговые марки
- Список благодарностей
Для просмотра документации необходимо авторизоваться на сайте
Войти
Памятка переводчика
1. Если смысл перевода понятен, то лучше оставьте как есть и не придирайтесь к словам, синонимам и тому подобному. О вкусах не спорим.
2. Не дополняйте перевод комментариями “от себя”. В исправлении не должно появляться дополнительных смыслов и комментариев, отсутствующих в оригинале. Такие правки не получится интегрировать в алгоритме автоматического перевода.
3. Сохраняйте структуру оригинального текста — например, не разбивайте одно предложение на два.
4. Не имеет смысла однотипное исправление перевода какого-то термина во всех предложениях. Исправляйте только в одном месте. Когда Вашу правку одобрят, это исправление будет алгоритмически распространено и на другие части документации.
5. По иным вопросам, например если надо исправить заблокированное для перевода слово, обратитесь к редакторам через форму технической поддержки.
Минутка «взрослого» программирования: как растягивать звук и как менять тональность.
Ни для кого не секрет, что если воспроизводить звук ускоренно, то изменится и тон его звучания и, разумеется, уменьшится длительность.
Это одно из самых простых преобразований, которые можно сотворить со звуком. Но что делать, если нужно изменить темп воспроизведения звука, но тон оставить прежним? Или наоборот, изменить тон, не трогая длительность?
Поскольку можно изменять одновременно тон и длительность, регулируя темп воспроизведения, обе задачи можно свести к одной: если, например, получится произвольно изменить темп, сохранив тон, то обратно пропорциональное изменение скорости воспроизведения соответствующим образом изменит тон, вернув темп в начальное состояние.
Чтобы было понятнее о чём пойдёт речь далее, небольшая видеозарисовка:
Ясно дело, без хитростей тут не обойдётся. На помощь нам приходит штука, известная как «Быстрое преобразование Фурье» (БПФ, или по-буржуйски 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. Представим такую ситуацию: у нас есть длинный протяжный звук, к примеру синусоида. Если частота синусоиды совпадает с частотой одной из гармоник, то на спектре будет присутствовать только эта гармоника, а вот фаза α будет меняться на каждой итерации, в зависимости от того на какую часть синусоиды пришлось начало окна. Не трудно догадаться что разница между фазами в двух соседних итерациях будет постоянной.
Если мы хотим растянуть наш этот монотонный гул во времени, то картина не должна измениться, он по прежнему останется таким же гулом и скорость изменения фаз будет ровно та же самая.
Представим теперь, что амплитуда гудящей синусоиды нарастает со временем. Теперь растягивание её во времени будет означать интерполяцию значений амплитуд (l), но разница между фазами в двух соседних итерациях опять таки будет изменяться с прежней скоростью.
Собственно эта идея и лежит в основе алгоритма:
1) для каждой гармоники на каждой итерации смотрим скорость изменения фазы — т.е. разницу с фазой из предыдущей итерации.
2) полученное значение скорости прибавляется к значению фазы такой же гармоники из предыдущей итерации выходного спектра, вне зависимости от того ускоряем мы, или замедляем воспроизведение.
3) значение амплитуды гармоники в выходном спектре получается путём интерполяции.
В пункте 2 возможны вариации — например, получать значение скорости фаз также при помощи интерполяции, но суть остаётся неизменной: важную роль играет не фаза, а скорость её изменения, которая воспроизводится в выходном спектре.
8. Получив выходной спектр, зеркально восстанавливаем его правую часть. При этом мнимая часть берётся с противоположным знаком.
9. После этого можно выполнить обратное быстрое преобразование Фурье над полученным спектром, благодаря тому что спектр зеркален, мнимые части результата будут равны нулю, а выходной сигнал будет вещественным.
10. Затем выходной сигнал нужно домножить на оконный фильтр. Здесь я использую ту же формулу, что и в п.3, с последующим делением на половину степени перекрытия.
11. Полученный сигнал складывается с учётом степени перекрытия, с результатами предыдущих и последующих итераций
12. Выходной сигнал готов, можно его выводить на аудиоустройство, или сохранять в файл — кому что ближе.
Изменение тона
Для изменения тона в N раз выполняется растяжение (т.е. замедление) звука в эти же N раз, а затем ускорение воспроизведения в N раз, что приводит к изменению тона при сохранении исходного темпа воспроизведения.
FFT (быстрое преобразование Фурье)
При работе со звуком возможно ли вычислить только часть звукового спектра, не вычисляя всю остальную? Если это возможно, то был бы рад ссылкам на движки, желательно Open Source.
Отслеживать
8,677 18 18 золотых знаков 74 74 серебряных знака 181 181 бронзовый знак
задан 13 ноя 2011 в 15:04
Vladimir4152 Vladimir4152
70 1 1 серебряный знак 7 7 бронзовых знаков
Что значит часть спектра? Как это? А opensource fft — «Fastest Fourier Transform in the West.» www.fftw.org
13 ноя 2011 в 15:52
Про этот вариант реализации реализацию — мне известно. Возникла проблема с интеграцией. Возможно ли получить часть спектра, то есть не вычислять все значения, а получить необходимый диапазон, в целях оптимизации.
13 ноя 2011 в 21:10
Что значит «необходимый диапазон»? fft вычисляет значения в диапазоне от 0 до половины частоты дискретизации. Вам нужна часть значений? Берите часть. Если бы речь шла просто о DFT (дискретном преобразовании Фурье), то имело бы смысл заморачиваться на эту тему, а оптимизации fft, как мне кажется, сводят смысл операции на нет. Но разумеется, Вы можете сами посчитать DFT по нужной части спектра без всяких библиотек, просто я не уверен, что овчинка стоит выделки.
14 ноя 2011 в 3:03
В целях оптимизации можно только уменьшить разрешение спектра, т.е. уменьшить число гармоник. Спектр будет строиться очень быстро на 100-200 гармониках, хотя, при этом, конечно теряется его точность. Выбирайте: скорость или детализация спектра.