Основы цифровой обработки сигналов: Спектральный анализ, Дискретная свёртка, Линейные стационарные системы
Рассмотрены 3 темы по основам цифровой обработки сигналов: спектральный анализ, дискретная свёртка, линейные стационарные системы.
В данном посте освещены 3 темы по основам цифровой обработки сигналов:
Спектральный анализ
Спектральный анализ сигнала можно осуществлять различными методами, но в данной публикации мы рассмотрим особенности спектрального анализа на основе БПФ.
Первая особенность – это векторный вход алгоритма БПФ. Алгоритм всегда работает с массивами значений входного сигнала, и на выходе получается вектор комплексных отсчётов спектра сигнала. Алгоритм требует векторный вход, но отсчёты сигнала зачастую поступают скалярно, то есть по одному. Поэтому перед непосредственным БПФ мы осуществляем буферизацию, или накопление отсчётов в вектор, а также определённые операции предобработки, о которых мы сейчас поговорим.
Входной вектор временных отсчётов мы называем окном, а выходной вектор мы называем точками или отсчётами БПФ. Размер входного окна тем или иным способом приравнивается к длине БПФ. Если во входном окне отсчётов больше, чем нужно – лишние отбрасываются, если меньше – входной вектор дополняется нулями.
От того, сколько точек БПФ мы получаем зависит разрешение по частоте, то есть та частотная сетка, в которой мы будем устанавливать амплитуды и фазы гармоник. Разрешение также зависит от частоты дискретизации сигнала. По факту мы берём весь частотный диапазон (от минус половины частоты дискретизации до половины частоты дискретизации) и заполняем его точками, в которых мы оцениваем спектр.
Представим, что мы анализируем действительный сигнал, и поэтому рассмотрим область положительных частот. Она ограничивается пределами от нуля до половины частоты дискретизации. В этом диапазоне мы рассматриваем nx2 точек БПФ. Расстояние между соседними точками равно разрешению по частоте. При конечном числе точек БПФ, спектр всегда будет выглядеть дискретным, но это не значит, что у реальных сигналов есть колебание на частоте fk, есть колебание большей амплитуды на частоте fk+1, а между ними ничего нет. Спектр непрерывен, но мы вынуждены рассматривать его в окрестностях конкретных частотных отсчётов.
На самом деле мы не можем подсчитать реальный спектр, но мы можем оценить его с конечной точностью. Часто мы оперируем понятием спектральной плотности мощности, то есть мощности сигнала, приходящейся на единичный интервал частоты, или попавшей в определённую полосу на оси частот.
Часто мы не сможем определить точное положение той или иной гармоники в спектре, так как она попадает между точками БПФ. Энергия этой гармоники может проявится в соседних точках БПФ. Этот эффект называется утечкой спектра. Избежать его полностью невозможно, но сгладить его влияние за счёт потери в точности в определении других параметров помогут оконные функции!
Оконная функция – это набор весовых коэффициентов, изменяющих значения амплитуды отсчётов сигнала во временной области. Когда мы осуществяем спектральный анализ большого сигнала, мы делим его на небольшие отрезки, которые мы взвешиваем оконной функцией и отправляем на вход алгоритма БПФ. Различные оконные функции позволяют выполнить оценки различных параметров спектра с большей точностью, но сами внесут искажения.
Идеальная оконная функция в частотной области имела бы форму прямоугольника, и получала только энергию в окрестностях своей частотной засечки. Но реальные функции в частотной области будут собирать энергию из соседних полос своими боковыми лепестками, или же не позволят различить близко расположенные гармоники за счёт слишком широкого основного лепестка.
Давайте выполним спектральный анализ в MATLAB при помощи встроенных функций. Входные данные со времен прошлой публикации не поменялись, это запись колебания струны 220 Гц. Но теперь мы воспользуемся встроенной функцией для оценки спектральной плотности мощности – pwelch.
В качестве входных аргументов эта функция может принимать размер входного окна, длину БПФ, количество перекрывающихся отсчётов (эта опция позволяет получить более сглаженный вид спектра, но пока что мы её не задействуем), и также мы указываем тип оконной функции. В нашем случае выбрана функция хэмминга, но вы можете поменять её на любую другую самостоятельно. Закомментированно здесь создание прямоугольного окна. Мы можем оценить вид выбранной оконной функции во временной области.
Затем мы запускаем основную функцию pwelch для оценки спектральной плотности мощности, а также подсчитываем текущее разрешение по частоте.
Пока что разрешение в районе 2,7 Гц, значение положения пиков мы рассматриваем с точностью до 1 Гц, так как округляем их функцией раунд. Ну и наблюдаем либо проблемы с настройкой гитары, либо же отклонения связанные с нехваткой точности и неровной частотной сеткой.
Количество точек входного сигнала меньше, чем количество точек БПФ, и дополнение нулями не поднимает фактического разрешения, а лишь интерполирует спектр. Так как мы не меняем nFFT – частотная сетка также не менялась, и положения пиков не изменились даже после того, как мы увеличили размер окна, но теперь мы наблюдаем большее количество гармоник, энергетики которых не хватало в коротком отрезке входного сигнала.
Если же мы увеличим количество точек БПФ, то частотная сетка теперь стала более плотной, и мы оцениваем положение пиков более точно. Теперь можно с уверенностью утверждать, что проблема не в разрешении, а в том, что я не умею настраивать гитару на слух. Спасибо, спектральный анализ!
Теперь давайте поговорим о том, что же такое нелинейные искажения с точки зрения спектрального анализа.
- Линейным считается изменения сигналов, приводящие к изменению амплитуд и/или фаз его спектральных компонент.
- Нелинейным называется изменения сигналов, приводящие к появлению новых спектральных компонент.
Подобное определение позволяет нам описать линейные и нелинейные системы с точки зрения того как они влияют на спектр сигнала.
Линейная система изменит амплитуда и фаза спектральных отсчетов, в то время как нелинейная система добавит те спектральные компоненты, которых не было в спектре входного сигнала.
Рассмотрим линейное и нелинейное искажение синусоиды. Во-первых, давайте посмотрим на нашу синусоиду, исходную синусоида 220 Гц во временной и частотной областях, также послушаем ее.
Можем наблюдать только одну спектральную компоненту. Затем попробуем выполнить линейное искажение сигнала, то есть изменение амплитуды. Давайте посмотрим, как эти сигналы отличаются на слух. Мы не услышали разницы в тоне, только услышали разницу в громкости.
Теперь давайте посмотрим как выглядит нелинейное искажения сигнала. В качестве примера нелинейного искажения мы используем здесь ограничение по уровню. Давайте посмотрим, что произошло со спектром ограниченного сигнала.
В нём появились дополнительные частотные компоненты, гармоники, которые мы точно также услышим и в окрасе или тоне нашего сигнала.
В завершении давайте поговорим об оконном преобразовании Фурье, которое позволяет нам получить спектрограмму сигнала. Оконное преобразования Фурье (по-английски short time Fouriern transform, STFT) – это метод оценки спектра сигнала, изменяющегося во времени. Во многих методах оценки стационарного сигнала мы также использовали оконные функции и разделение сигнала на окошки, но после этого мы сигнал склеивали, осредняли и наблюдали спектр всего сигнала. В то время как оконное быстрое преобразование Фурье позволяет нам получать кусочный спектр, который мы затем накладываем на временную ось. Спектрограммы называются, по сути, трёхмерная картина изменения спектра во времени.
Рассмотрим пример классического нестационарного во времени сигнала, то есть речи человека. Мы будем наблюдать спектр и спектрограмму сигнала во времени при помощи инструмента Spectrum Analyzer. В качестве источника аудиоданных я буду использовать микрофон моего компьютера. Для обработки и чтения я буду использовать системные объекты. Отдельные системные объекты я буду использовать для захвата аудиопотока с микрофона, отдельные системные объекты я буду использовать для построения спектра и спектрограммы. Эти системные объекты входят в состав расширения DSP System Toolbox. Когда я их подготовил, я могу начать обращаться к ним в основном цикле. Цикл мы запускаем стандартной конструкцией for, запускаем его на исполнение 10 тысяч раз, и мы будем наблюдать в живую, как на спектре и спектрограмме отражается моя речь. Запустим.
Вот так на спектре и спектрограмме выглядит речь человека. Если мы хотим попробовать что-то интересное, можем посмотреть, как на спектре и спектрограмме выглядит, к примеру, гласные звуки. Мы наблюдали гармоническую структуру. Если мы попробуем издать согласные звуки, например “Ш”. Мы увидим шумоподобный сигнал, который занимает большую часть спектра. Звук “С” чем-то напоминает синий или фиолетовый шум. А если я попробую посвистеть, то я могу увидеть явный пик. Подобный пример, наверное, будет очень интересно проделать самостоятельно.
Дискретная свёртка
В этой публикации мы поговорим об одном из ключевых понятий в цифровой обработке сигналов – дискретной свёртке. Как и большинство операций, которые мы рассматривали в предыдущих публикациях, свёртка может быть как непрерывной, так и дискретной. По определению, свёртка может вычисляться для двух функций или последовательностей, и на выходе мы получим третью, результирующую, функцию или последовательность.
Непрерывная свёртка является одним из видов интегрального преобразования, и по своей сути весьма близка к корреляционной функции. Напомню, в случае с корреляционной функцией, мы оценивали степень схожести одной функции с другой в зависимости от временного сдвига. Проще говоря, проезжались одной функцией по другой, и считали площадь перекрытия. В случае же со свёрткой, мы оцениваем степень схожести одной функции с зеркально отражённой и сдвигаемой во времени второй функцией.
Давайте рассмотрим пример вычисления дискретной свёртки, то есть свёртки двух последовательностей. Как вы уже знаете, при переходе от непрерывных сигналов к дискретным, операции интегрирования заменяются более простыми перемножениями и суммами. Рассмотрим две дискретные последовательности по пять отсчётов каждая: х со значениями 4 2 1 2 3, h - 1 2 4 2 3. Для того, чтобы показать временной сдвиг – дополним первую последовательность нулями. Вторую последовательность отражаем зеркально, и поэлементно перемножаем с первыми пятью отсчётами, которые у нас на входе. Элементы результата перемножения складываем между собой. Получаем число 4.
Затем сдвигаем первую последовательность на один отсчёт, и повторяем операции перемножения и суммирования.
Таким образом мы вычисляем элементы выходной последовательности игрик!
Как вы заметили, мы повернули зеркально вторую последовательность h, а двигаем первую x. На самом деле принципиальной разницы нет, так как важен именно относительный сдвиг, но подобное представление мы взяли неспроста. Первую последовательность мы можем рассматривать как входной сигнал, а вторую – как некую систему, через которую этот сигнал проходит. Результат линейной дискретной свёртки, а именно её мы сейчас вычисляли, содержит больше отсчётов, нежели каждая из сворачиваемых последовательностей.
Для того, чтобы быстро посчитать свёртку в MATLAB, мы можем воспользоваться встроенной функцией conf. Давайте создадим эти две последовательности в командном окне, и посчитаем их свёртку в одно действие:
Как мы видим, результат соответствует рассчитанному ранее.
Давайте перейдём от математических абстрактных вычислений к практическому применению операции свёртки для цифровой обработки сигналов. Но для этого надо будет на минутку вернуться назад в аналоговый мир. Вспомним теорию цепей, базовые вещи – интегрирующая RC-цепь. Допустим, мы включаем источник напряжения на некий промежуток времени, а затем выключаем его. Во временной области зависимость напряжения источника будет представлять собой прямоугольник. А цепь мы описываем её импульсной характеристикой, то есть реакцией цепи на дельта функцию.
Напоминаю, дельта-функция, или функция Дирака, это функция, не равная нулю только в начале координат, или в момент времени t равный нулю. То есть бесконечно короткий и бесконечно мощный импульс. Так вот, реакция системы на подобный импульс и называется импульсной характеристикой. С её помощью мы можем описывать динамику системы.
И для того, чтобы вычислить выходной сигнал, в нашем случае напряжение на конденсаторе – нам необходимо осуществить свёртку входного сигнала (напряжения источника) с импульсной характеристикой цепи. Свёртка используется для вычисления выхода цепи/системы при известном входном сигнале и известной характеристике цепи/системы! Запомните, это справедливо и для дискретных систем и сигналов, в чём мы убедимся в последующих публикациях.
А в завершении этой публикации мы рассмотрим пример применения свёртки для поиска сигналов. Если система имеет импульсную характеристику, обратную значениям входной последовательности, то в момент совпадения отсчётов выход операции свёртки будет максимальным. На этом основывается принцип согласованной фильтрации. Согласованные фильтры применяются в тех случаях, когда нам не важно сохранить форму принимаемого сигнала, а важно определить его наличие в эфире и узнать время его прибытия. Это одна из типичных задач радиолокации.
В том случае, если мы заранее знаем форму сигнала, который мы ожидаем принять, нетрудно подобрать коэффициенты фильтра, то есть значения второй последовательности для операции свёртки, такими, чтобы при зеркальном отражении они равнялись отсчётам сигнала. В этом случае в момент прихода на вход системы зашумлённого сигнала значение на выходе фильтра будет максимальным.
Подобными образом можно обнаруживать сигналы на фоне шумов. Рассмотрим пример в MATLAB.
В примере мы будем искать некую случайную последовательность X при помощи согласованного фильтра с импульсной характеристикой h. Обратите внимание, что импульсную характеристику h мы формируем при помощи отражения последовательности X. Функция fliplr отражает элементы массива. Затем мы добавление случайные элементы до и после нашей последовательности, а также добавляем ещё дополнительный шум. Мы знаем, куда мы поместили нашу последовательность, поэтому запишем положение конца нашей последователи в переменную expected_num. Затем вычисляем свертку нашего зашумленного сигнала с импульсной характеристикой h. Максимальное значение свёртки должно показать положение конца нашей последовательности. Давайте проверим.
На данном графике представлен зашумлённый сигнал и исходная последовательность.
На данном графике представлена наша свертка. Как видите она имеет явный максимум. Можно попробовать понажимать кнопку Run ещё несколько раз, чтобы убедиться что характер шума не влияет на положение искомой последовательности. Если мы сдвинем нашу последовательность на несколько отчетов, то мы точно также увидим, что наша функция свертки успешно отслеживает её положение даже на фоне шумов.
Линейные стационарные системы
В данной публикации мы поговорим о линейных стационарных системах. Для начала давайте вспомним, что мы называем системой. Система – это некий объект, который каким-то образом преобразует или анализирует входной сигнал, или воздействие, и порождает выходной сигнал, или отклик системы.
Что же такое линейная стационарная система? Это система, которая удовлетворяет условиям линейности и стационарности. Линейные стационарные системы бывают непрерывными и дискретными, но в рамках данного курса мы будем рассматривать только дискретные системы.
Давайте разберёмся с определениями. Система называется линейной, если выход или отклик системы получается из линейной комбинации отсчётов входного и выходного сигналов, в том числе задержанных во времени. Под линейной комбинацией здесь подразумевается взвешенная сумма. Отсчёты могут умножаться на различные коэффициенты, потом складываются между собой – и мы получаем значение отсчёта выходного сигнала. Потом переходим к следующему дискрету времени.
Стационарность – это свойство системы выдавать один и тот же отклик на одинаковое входное воздействие вне зависимости от времени. То есть мы не можем получить разные сигналы на выходе, если подали один и тот же импульс вчера и сегодня. Результат будет одинаковым.
Ну о дискретности я ради экономии времени лишний раз говорить не буду. Давайте лучше познакомимся со способами описания линейных стационарных систем.
Первый способ описания, с которым мы познакомимся – это разностное уравнение. Этот способ оперирует непосредственно отсчётами сигнала во временной области. По сути, разностное уравнение показывает, как из задержанных и текущих отсчётов входа, а также задержанных отсчётов выхода посчитать текущее значение на выходе системы.
Как я и упоминал, значения отсчётов могут быть умножены на различные весовые коэффициенты. И при записи значений отсчётов входа (x) и выхода (y) в виде взвешенных сумм задержанных отсчётов, нам на самом деле нужно только знать значения этих самых коэффициентов. Например, в нашем уравнении нет слагаемого x с задержкой на 3 отсчёта. Значит ему соответствет весовой коэффициент ноль. Подобным образом мы накладываем вектор коэффициентов на равномерно нарастающие задержки сигналов, и можем всегда восстановить разностное уравнение. За входные отсчёты отвечают коэффициенты b, за выходные – а. Давайте рассмотрим конкретное разностное уравнение.
Текущее значение выхода системы – это y(n). Текущее значение входа – x(n). У нас также присутствуют предыдущие значения входа и выхода. Решим это уравнение отсносительно текущего значения выхода. Оно получается из суммы входа, входа задержанного на один такт, и выхода, задержанного на один и два такта.
Соответствующие коэффициенты при входных отсчётах, вектор b – 0.75, 0.5
Выходные коэффициенты – 1, 0.2, 0.1
Разностное уравнение можно представить схематически! На вход схемы один за одним поступают отсчёты входного сигнала. В каждый момент времени для того, чтобы вычислить икрик эн нам нужно взять текущее значение входа, умножить на 0.75, заглянуть в регистр (это такая маленькая ячейка памяти, она может хранить значение сигнала до следующего такта), взять оттуда значение входа с предыдущего такта, умножить его на 0.5. То же самое происходит и с сохранёнными предыдущими выходными значениями, они умножаются на соответствующие коэффициенты, и потом это всё идёт на сумматор. Результат, сумма, и есть наш выход y(n).
По разностным уравнениям можно составлять схемы линейных дискретных систем и реализовывать их в железе.
Мы продолжаем рассматривать линейные стационарные системы. И, как вы помните из предыдущей публикации, динамику системы мы можем описать импульсной характеристикой. К дискретным системам это так же относится. Но у дискретной системы импульсная характеристика получается в результате реакции её на единичный дискретный импульс. Это подобие функции Дирака для дискретных сигналов, у неё есть всего один отсчёт равный единице в момент времени t0.
Если нам известна импульсная характеристика дискретной системы, мы можем определить отклик для любого входного воздействия при помощи операции дискретной свёркти.
Что если у нас есть коэффициенты разностного уравнения, и нам нужно быстро вычислить импульсную характеристику системы? Переходим в MATLAB и пользуемся встроенной функцией. Для начала задаим вектора a и b. А затем запишем h=impz(b,a). Можно визуализировать характеристику функцией stem.
Но давайте вернёмся к слайдам и рассмотрим ещё один способ описания ЛДСС. Для этого нам нужно будет познакомиться с z-преобразованием. Также, как преобразование Фурье является основой анализа сигналов, z-преобразование – это основа анализа систем. Оно ставит в соответствие значениям сигнала во временной области отсчёты в комплексной частотной области. Задержка сигнала на k-отсчётов становится умножением на z в степени k. Умножение на коэффициент отправляется в знаменатель, и самое главное – операция свёртки заменяется простым умножением.
Как вы видите, вычисления становятся гораздо проще и эффективнее. Поэтому часто анализ сложных систем осуществляют именно в z-области, и затем проводят обратное преобразование во временную область.
Выход системы, который мы рассчитывали при помощи свёркти входа с импульсной характеристикой в z-области высчитывается как перемножение с z-формой импульсной характеристики. Отношение входа и выхода в зэд-облатси называется передаточной функцией системы. Это дробная функция с полиномами в числителе и знаменателе.
Давайте для примера переведём наше разностное уравнение в z-область. Теперь мы можем записать передаточную функцию, коэффициенты при x пойдут в числитель, а при y – в знаменатель. Как вы заметили, в числителе и знаменателе полиномы, которые так же полностью определяются знакомыми нам коэффициентами. Как и в случае разностного уравнения, нам достаточно знать только вектора a и b.
И в завершении поговрим об ещё одном интересном описании линейных стационарных систем, так называемом нуль полюсном описании. Переменная z – это некоторая комплексная величина, которую можно отобразить на комплексной плоскости. А передаточная функция определена для всех значений z, то есть может быть представлена трёхмерной поверхностью на комплексной плоскости. Мы сейчас не будем рисовать эту поверхность, но определённые ключевые точки можем определить. А именно точки, в которых наша передаточная функция равна нулю, и точки, в которых она равна бесконечности. Как вы знаете, дробь равна нулю, когда её числитель равен нулю, и равна бесконечности, когда её знаменатель равен нулю. Мы можем решить уравнения числителя и знаменателя, определить корни. Корни числителя – это так называемые нули функции, а корни знаменателя – полюса.
Для того, чтобы было проще решать – домножим числитель и знаменатель на зэд в квадрате, так как минимальная степерь у нас – минус вторая. И получаем корни для числителя и знаменателя. Их мы можем отразить на комплексной плоскости. Нули помечаются кружками, полюса – крестиками. Это отображение называется нуль-полюсной диаграммой.
Построим её в MATLAB. Корни числителя и знаменателя передаточной функции можно вычислить в одно действие – для этого воспользуемся функцией tf2zpk. Нули помещаются в вектор z, полюса – в p. Теперь отобразим их на диаграмме функцией zplane.
А о том, как связана нуль-полюсная диаграмма с амплитудно частотной характеристикой системы – в следующей публикации.
Комментарии