Основы цифровой обработки сигналов: Приближение сигнала функцией, Частотное представление сигнала, Преобразование Фурье
Рассмотрены 3 темы по основам цифровой обработки сигналов: приближение сигнала функцией, частотное представление сигнала, преобразование Фурье.
В данном посте освещены 3 темы по основам цифровой обработки сигналов:
Приближение сигнала функцией
Пришло время поговорить о приближении сигнала аналитической функцией. Сперва давайте вспомним математическую модель или представление реального сигнала из прошлых публикаций. Я говорю о модели типа сигнал + шум, или детерминированный сигнал, описанный функцией и случайный процесс. Отсчёты температуры воды, нагреваемой на горелке, хорошо аппроксимируются линейной зависимостью, или полиномом первого порядка.
В отличие от реального измеренного сигнала, который определён только в нескольких точках на конечном временном интервале, аналитическая функция определена в любой точке на временной оси. Мы можем использовать значения этой функции, чтобы оценить величину сигнала между реальными отсчётами или измерениями, а также благодаря ей мы можем выйти за пределы временного интервала измерения, и предсказать значение сигнала до или после наблюдения.
На самом деле оценить величину сигнала между отсчётами мы можем и при помощи простых методов интерполяции, например, таких, как линейная интерполяция. Она, кстати, используется при отрисовке графиков дискретных последовательностей. Вы видите, что точки на графике соединены линиями. Но для предсказания такие методы не очень годятся, и аналитическая функция, приближающая измеренные данные, видится здесь более универсальным решением.
Однако, задача выбора приближающей функции и её параметров – весьма нетривиальна. Какая функция лучше подойдёт для конкретных значений, как мы можем оценить качество приближения и выбрать лучший вариант? Поможет нам приложение Curve Fitting из состава Curve Fitting Toolbox, расширения MATLAB.
В этом приложении можно быстро попробовать различные способы интерполяции и приближения ваших данных, а также оценить численные метрики качества приближения.
Загрузим данные и откроем приложение под названием Curve Fitting во вкладке APPS в разделе математика, статистика и оптимизация. Выберем Untitled fit1 и загрузим данные по оси X и Y.
Он предлагает нам приближение полиномом, в данном случае полиномом первой степени, то есть линейная зависимость. Также нам доступны величины отклонения от нашего приближения и метрики, такие как средние квадраты отклонений или корни среднего квадрата отклонения. По этим метрикам мы можем оценивать качество нашего приближения. Также нам доступны самые разнообразные функции: экспоненциальная функция, приближение рядом фурье, Gaussian и так далее.
Перебирая различные функции можно найти ту, которая лучше всего подходит под наши данные. Конечно же также тут доступны различные методы интерполяции. Как вывести для них метрики не обозначается, потому что это не приближение аналитической зависимостью. К примеру можем выбрать сумму синусоид и посмотреть как сумма из 8 синусоид может приблизить наши данные.
Особо важным считаю отметить возможность выбора пользовательского уравнения, то есть пользователь сам может написать уравнение, которое, как он считает, лучше приближает эти данные. К примеру, можно поменять функцию экспоненты на функцию косинус и посмотреть, как синусоида приближает наши данные.
Помимо работы с одномерными данные иы также можем приближать двумерные данные. Загружаем значение X,Y, Z. По умолчанию инструмент предлагает нам интерполировать данные, но тут также доступна опция приближения полиномами, причём с разными степенями по оси X и Y.
Как вы видите здесь также мы можем посмотреть разницу между реальными значениями и приближающей плоскостью или поверхностью. Если мы поменяем порядок полинома по оси Y на тройку, то мы увидим приближение поверхности типа лист бумаги, можем поменять степень полинома по оси X и посмотреть насколько улучшились метрики приближения.
Приближение экспериментальных данных аналитической функцией используется во многих областях, в том числе в анализе данных.
Частотное представление сигнала
Пришло время поговорить о представлении сигнала в частотной области. В предыдущих публикациях мы рассматривали в основном временное представление сигнала, то есть зависимость значений сигнала от момента времени.
Частотное представление – это зависимость измеряемых параметров сигнала от частоты. Вы спросите, как это возможно, ведь частота является одним из измеряемых параметров сигнала. Для простых периодических сигналов это справедливо, но мы будем рассматривать более общий случай. Проведём аналогию с временным представлением. Мы можем отразить значения сигнала на временной оси.
У большинства сигналов мы можем также описать зависимость энергии от частоты, и подобное частотное представление именуется спектром сигнала. Как и сигнал во временной области, спектр может быть непрерывным или дискретным, конечным или периодическим.
Давайте рассмотрим пример, к оторый поможет понять, что же это такое - зависимость уровня сигнала от частоты.
Аудио-эквалайзер — это устройство, позволяющее усиливать или ослаблять аудио-сигнал в определённых полосах частот. Усиление сигнала в области низких частот добавит басов в проигрываемом треке, усиление в области высоких частот сделает высокие звуки громче и добавит обертона. Рассмотрим пример в MATLAB.
Загрузим в рабочее пространство музыкальный файл и прослушаем его командой sound. Теперь откроем Sinal Analyzer и оценим сигнал во временной и частотной областях. Напоминаю, для отражения спектра нам на вкладке Display надо выбрать кнопку Spectrum.
Рекомендую делать копию сигнала для наглядного сравнения. Изменения претерпит именно копия. Мы пропустим её через фильтр верхних частот, кнопка Highpass на вкладке аналайзер. Мы подавим все частоты ниже 1 кГц. И тут же наблюдаем изменения сигнала во временной и частотной областях.
Но нам интересно услышать результат. Для этого экспортируем значения копии сигнала в рабочее пространство и снова выполним команду Sound. Басов нет, мы слышим только средние и высокие частоты.
Музыкальный сигнал, который мы анализировали и обрабатывали, является сложным сигналом. В предыдущих публикациях мы говорили, что сложные сигналы можно представлять суммой простых сигналов. По факту это и лежит в основе частотного представления сигнала.
Синусоиды с разными частотами, помноженные на весовые коэффициенты, при сложении могут дать нам сложный сигнал. Но помимо амплитуд этих синусоид на форму суммы также влияет фаза каждой из составляющих. Если мы немного сдвинем сигналы относительно друг друга по временной оси – результат сложения будет заметно отличаться. Поэтому мы говорим, что весовые коэффициенты должны учитывать как изменения амплитуды, так и изменения фазы.
Убедимся в этом на знакомом примере формирования прямоугольного сигнала из нескольких синусоид.
Отличие этого скрипта от того, что мы рассматривали ранее – вектор начальных фаз для каждого из слагаемых. Вначале, когда все начальные фазы равны нулю, мы наблюдаем результат, близкий по форме к прямоугольному сигналу. Но если мы укажем ненулевые начальные фазы для наших синусоид, при этом не поменяв амплитуды, то форма из суммы меняется кардинально. Мы убедились, что информация о фазе важна не менее, чем информация об амплитуде элементарного сигнала.
В предыдущих публикациях мы также говорили об эквивалентом представлении сигнала набором весовых коэффициентов. Давайте применим это к частотному представлению. Каждая точка на частотной оси соответствует синусоиде со своей амплитудой и фазой. Поэтому каждый коэффициент должен содержать в себе две величины.
Для подобного представления весьма удобно пользоваться комплексными числами. Модуль коэффициента будет описывать амплитуду синусоиды, а аргумент – начальную фазу. В нашем случае в формулах переменные x и y соответствуют положению точки на комплексной плоскости в декартовых координатах, а формулы вычисления модуля и аргумента – переводу из декартовых в полярные координаты. Давайте рассмотрим комплексную плоскость.
Точкой на комплексной плоскости мы можем отразить значение комплексной величины. По горизонтальной оси откладывается значение действительной части, а по вертикальной – мнимой. Когда мы переходим в полярные координаты, то значение комплексной величины описывается вектором, откладываемым из начала координат до точки на плоскости. Длина вектора – это модуль или амплитуда комплексного числа, а угол поворота фи, отсчитываемый от горизонтальной оси – аргумент или фаза.
Из курса тригонометрии мы можем вспомнить формулу Эйлера, связывающую функции синуса и косинуса с комплексной экспонентой. Косинус – это проекция комплексной экспоненты, или вращения вектора против часовой стрелки на комплексной плоскости, на действительную ось, а синус – проекция на мнимую ось. Причём скорость вращения вектора определяет частоту синусоиды, длина вектора – амплитуду, а начальный угол – начальную фазу.
Если мы рассмотрим сумму синусоид, как сумму вращающихся векторов на комплексной плоскости, но проекция фигуры вращения, описываемой суммой на одну из осей даст нам похожую картину. Мы опять наблюдаем формирование прямоугольного сигнала. В сети вы можете найти примеры формирования пилообразного сигнала, это достигается за счёт изменения фаз и частот слагаемых, что в свою очередь приводит к изменению фигуры вращения.
Комплексные числа – это математическая абстракция, помогающая нам описывать некоторые явления физического мира более простой математикой. Но в реальном мире мы оперируем действительными сигналами, у них нет мнимой части.
Рассмотрим математическое описание действительной синусоиды, и разложим его по формуле Эйлера. Мы видим два слагаемых с зеркальными начальными фазами и частотами.
Если мы рассмотрим эти слагаемые на частотной оси, то мы увидим две гармоники с одинаковыми амплитудами, одинаковыми по значению, но разными по знаку начальными фазами, и находящиеся по разные стороны начала координат. Это – спектр действительной синусоиды. Он содержит информацию о фазе и частоте гармоник на выбранных частотах. Если нас не интересует информация о фазе, а важно только оценить распределение энергии сигнала в зависимости от частоты – то можно рассмотреть амплитудный спектр, то есть модуль полного спектра. Он будет вещественным.
Но что же это за отрицательная частота, на которой мы наблюдаем гармонику? Отрицательная начальная фаза для синусоиды вполне понятна, но как мы можем представить отрицательную частоту?
Отрицательная частота не имеет физического смысла, как например отрицательная амплитуда или фаза, но имеет математический смысл на комплексной плоскости. Вращение вектора по часовой стрелке соответствует отрицательному приращению угла, то есть отрицательной частоте.
Сумма двух векторов одинаковой амплитуды, вращающихся в противоположных направлениях с одинаковой скоростью, будет всегда лежать на действительной оси. При условии, конечно, что их начальные фазы равны по модулю и отличаются знаком, то есть вращающиеся вектора всегда зеркальны относительно действительной оси.
Подобным образом мы представляем вещественную синусоиду суммой двух комплексных вдвое меньшей амплитуды, что полностью соответствует формуле Эйлера. Спектр любого вещественного сигнала, состоящего из большого количества синусоид, всегда будет симметричен относительно нулевой частоты.
Итак, мы поняли, что представление сигнала в частотной области связано с разложением сложных сигналов на простые составляющие с учётом их амплитуд и фаз. Но какой же математический аппарат помогает нам с задачей декомпозиции? Об этом мы узнаем в следующей публикации.
Преобразование Фурье
В данной публикации мы поговорим о преобразовании Фурье. Это одно из важнейших преобразований, которое используется во многих областях науки, но для задач обработки сигналов мы используем преобразование Фурье, как способ разложения сигнала на частоты и амплитуды, то есть для перевода его из временного представления в частотное. Для непрерывных сигналов мы можем применять интегральную форму, но когда мы говорим о последовательностях чисел, или дискретных сигналах, то мы применяем дискретное преобразование Фурье.
По факту результатом дискретного преобразования Фурье будет набор комплексных отсчётов спектра. Но для вычисления каждого такого отсчёта нам требуется перемножить и просуммировать весь набор значений сигнала во временной области со значениями комплексной амплитуды соответствующей частоты. Отсчёты спектра на выбранной частотной сетке соответствуют элементам ряда Фурье, в который была разложена входная последовательность.
Существуют также обратные преобразования Фурье, переносящие сигнал из частотной области во временную, и подобным образом мы можем переводить сигналы в частотную область и обратно без потери точности, но только в том случае, если количество элементов в ряде Фурье бесконечно. Как вы понимаете, в реальных вычислителях подобное разложение выполнить невозможно, поэтому мы всегда ограничиваем количество элементов ряда Фурье. Это приводит к ошибкам при преобразовании. Давайте посмотрим, как они могут проявляться.
Рассмотрим уже знакомый нам пример формирования прямоугольного сигнала из суммы синусоид. Посмотрим на картинки справа налево.
При уменьшении количества гармоник, форма сигнала отдаляется от прямоугольной. Чем меньше элементов ряда Фурье мы используем – тем сложнее нам представить мгновенный переход сигнала из положительного в отрицательные значения. Подобный скачок может быть описан только синусоидой с бесконечно большой частотой, а мы в своём описании ограничены самой большой частотой в выбранном конечном ряде Фурье.
Искажения сигнала, связанные с ограничением ряда Фурье, называются эффектом Гиббса. При возникновении этого эффекта мы наблюдаем как размытие резких перепадов, так и нежелательные выбросы, и колебания в районе этих перепадов. Вы встречали эффект Гиббса в повседневной жизни, к примеру, в цифровых изображениях. Алгоритм сжатия jpeg основан на преобразовании Фурье, и на границах перепадов яркости пикселей мы можем наблюдать размытие, выбросы и колебания амплитуды. Вы могли не только видеть, но и слышать эффект Гиббса. При сжатии музыкальных файлов колебания, возникающие до резкого перепада, могут рассматриваться как пре-эхо, например при резком вступлении партии или ударе по тарелке.
Теперь давайте поговорим о быстром преобразовании Фурье (БПФ). БПФ – это эффективный алгоритм для вычисления дискретного преобразования Фурье (ДПФ). Оптимальная длина временной последовательности для БПФ – в степени двойки.
Есть несколько вариаций алгоритма, но мы познакомимся с наиболее распространённой его версией, алгоритмом Кули-Тьюки. Он основывается на том факте, что вычисление ДПФ из, к примеру, восьми точек может быть выполнено отбором чётных и нечётных значений входного сигнала и выполнением двух четырёхточечных ДПФ с дальнейшей операцией так называемого «прореживания по времени». Вывод формул алгоритма доступен в интернете, поэтому с ним я рекомендую ознакомиться самостоятельно. Мы же очертим только ключевые моменты.
Операция прореживания во времени требует умножения выходных отсчётов ДПФ на так называемые «поворачивающие коэффициенты», которые по сути являются значениями комплексной экспоненты. Причём этих значений на окружности берётся столько, сколько точек ДПФ мы хотим объединить.
Если мы рассмотрим ДПФ по двум точкам, то и поворачивающих коэффициентов, то есть точек на окружности комплексной экспоненты, нам нужно два. И их значения будут сугубо действительными – плюс один и минус один. В случае ДПФ по четырём или восьми точкам мы будем задействовать мнимые и комплексные коэффициенты, но простейший двуточечный вариант позволяет нам просто обойтись операциями сложения и вычитания отсчётов.
Операция эта получила название бабочки за внешнее сходство графа с летающим насекомым. И использование прореживания ДПФ по времени подобными бабочками позволяет восьмиточечный ДПФ вычислить из двух четырёхточечных, или же из четырёх двухточечных. Подобным образом можно разложить любой ДПФ длительности степени двойки.
Рассмотрим применение алгоритма БПФ в MATLAB.
Попробуем получить спектр аудиофайла при помощи встроенной функции fft. Аудио файл – это запись колебания гитарной струны, нота ля, 200 Гц. Послушаем его.
Выполним пераобразование Фурье функцией fft и попробуем отобразить результат командой плот.
Мы наблюдаем вот такую интересную картину. Связано это с тем, что функция плот строит нам график значений комплексной величины на комплексной плоскости. Построим амплитудный спектр сигнала. Воспользуемся функцией abs для вычисления модуля комплексного вектора.
Сигнал у нас действительный, выход функции эфэфти учитывает как область положительных частот, так и область отрицательных частот. Но для того, чтобы представить спектр в привычном нам виде, зеркально относительно нулевой частоты, можно воспользоваться функцией fftshift. Также ограничим пределы – от –2 до +2 кГц. Спектр симметричен, поэтому нам достаточно проанализировать его только в области положительных частот.
Мы наблюдаем гармоническую структуру в спектре сигнала, то есть помимо основной гармоники мы видим обертона.
Воспользуемся функцией findpeaks для определения положения локальных максимумов спектра. Значения частот мы немного округлили, и теперь явно прослеживаются гармоники – производные основной частоты в 220 Гц. Можно проверить, каким нотам соответствуют данные колебвания в таблице.
Функция ifft выполняет обратное преобразование Фурье. Мы можем перевести значения отсчётов спектра во временную область и послушать сигнал. Так как мы брали достаточно большое количество точек для представления сигнала в частотной области, особых искажений сигнала после преобразований мы не услышали.
БПФ и обратное БПФ применяются во многих задачах ЦОС, например в таких как фильтрация в частотной области, сжатие, осуществление быстрой дискретной свёртки, а также конечно же для спектрального анализа. Особенности применения БПФ для спектрального анализа мы рассмотрим в следующей публикации.
Комментарии