Теория
Для начала немного теории. Как известно все в подобных анализаторах используется быстрое преобразование Фурье и часто говорится, что ДПФ в подобных конструкция использовать нельзя, только БПФ да и то в на ассемблере. Я же использовал вместо этого дискретное преобразование Фурье(ДПФ) и преобразование по Уолшу. И в этой статье докажу, что можно использовать даже не только БПФ, а ДПФ написанный на С. Но сначала по порядку как из ДПФ получить простую функцию ДПФ и по Уолшу. ДПФ классически выглядит следующим образом:
Так как у мк мало ресурсов, то заменяют cos и sin на массивы размерностью N. Кроме того мк 8 разрядный и целесообразнее массивы хранить в виде 8 разрядных значений. Так как cos и sin меняется от -1 до 1, то лучше всего это диапазон увеличить в 127 раз, так как переменная 8 разрядная знаковая может хранить в себе значения от -127 до 127. Таким образом с учетом преобразований формулы будет:
где m меняется от 0 до N-1 с шагом равный k, когда m становится больше N, m уменьшают на N-1. Всего испльзуется 12 каналов, так что мк по силу ДПФ на такое маленькое количество каналов.
Например имеем 512 отсчетов АЦП нужно посчитать мнимую и действительную части для 150Гц при частоте дискретизации 19200 Гц:
Таким образом реальная и мнимая части получаются гораздо быстрее нежели традиционным способом, но в 127 раз больше. Для того, чтобы получить их реальный значения нужно поделить на 127, но у мк нет деления, поэму гораздо рациональнее будет не делить, а сдвинуть! Один сдвиг эквивалентен делению на 2. То есть если сдвинуть7 раз число то по сути поделили на 128! Так как потери в точности уже были неизбежны, то деление на 128 картины не изменит.
Дискретное преоразование Фурье для 150 Гц при частоте дискретизации 19200 Гц тогда выглядит следующим образом:
Для Уолша заменяем синусоиду и косинусоиду на меанды соответствующих периодов. То есть для sin от 0 до 180 градусов будет 1 а от >180 до 360 будет -1. Соответственно для синуса от 0 до 90 это 1, от 90 до 270 это -1 и от 270 до 360 это 1. Тем самым все вычисление мнимой и действительной части будут простым накапливанием сумм и разностей значения АЦП. То есть когда например синус равен 1, то значение АЦП прибавляется, а когда -1 отнимается. Недостаток такого решения заключается опять же в погрешности, которая неизбежно увеличивается и достигает 20%. Но так как в моей конструкции всего 8 значений то опять же существенно разницы мало кто заметит.
Пример реализации расчета мнимой и действительной части для 150 Гц при частоте дискретизации 19200 и 512 отсчетов:
Таким образом получаем довольно быстро мнимые и действительные части без процедур умножения.
И так получив мнимую и действительную части необходимо найти амплитуду спектра. Для этого необходимо найти корень из сумм квадратов мнимой и действительной части. Но если воспользоваться функцией из библиотеки math извлечение получится долгим и функция к тому же съест не хилый кусок от ПЗУ. Немного покопавшись в интернете я нашел элегантную функцию которую потом еще немного упростил в силу того, что она оперирует маленькими значениями. Вот это функция:
Сравнив эту функцию и функцию из библиотеки math пришел к выводу, что ее точности вполне хватает, чтобы результат был одинаков. Сама функция весит 2% против 12% от ПЗУ мк. Кроме того вычисляет гораздо быстрее.
Но как же получилось, что мк успевает расчитать 12 каналов да еще и в ДПФ. Кроме всех ухищрений со сдвигом вместо деления и быстрой функции квадрата есть еще одна уловка. Про которую я сейчас раскажу. Дело в том, что чем выше частота выделения тем уже полоса пропускания фильтра, так как переход cos и sin убыстряется и число периодов растет. А чем больше таких проходов cos и sin тем уже полоса пропускания. Например для частоты 150 Гц cos и sin повторяются 4 раза, а для 1,2 кГц cos и sin повторяются уже 32 раза. Отсюда видно, что для того чтобы на всех диапазонах полоса пропускания была равномерной и охватывал всех диапазон частот число отсчетов с ростом частоты фильтрации надо уменьшать. Например для 150 Гц бурутся все 512 отсчетов, для 600 Гц 256 отсчетов, а для 2,4 кГц 32 отсчета и так далее. Не трудно заменитить, что уменьшая число отсчетов с ростом частоты круто увеливается скорость ДПФ, так как умножений и сумм уже нужно делать гораздо меньше.
Практическая реализация
И так теоретическая часть подготовлена можно приступать к описанию конструкции. Вся конструкция состоит из одного микроконтроллера, 4-х транзисторов, нескольких конденсаторов и много резисторов. Резисторов лучше поставить много, хотя можно ограничиться только резисторами по горизонтали, т.е. по одному на каждый вывод порта. Схема классическая кроме единственного, что я использовал по 3 порта за 1 проход динамической развертки вместо 1 как везде делают. Это позволило уменьшить частоту развертки и уменьшить число транзистров до 4. Получилась фактически шкала на 24х4.
Анализатор спектра работает на частоте дискретизации 19,2 кГц от кварца на 16 МГц.
Анализатор спектра рассчитывает спектральные амплитуды следующих частот:
9,6 кГц, 4,8 кГц, 2,4 кГц,1,6 кГц, 1,2 кГц, 800 Гц, 600 Гц, 500 Гц, 400 Гц, 300 Гц, 150 Гц, 75 Гц. Программа проверялась и для 33 Гц и ДПФ успевал при тома что размерность cos и sin становится равный 512, но решил ограничится 75 Гц.
Здесь имеются частоты которые не кратны 2 в n-й степени, но тем не менее вычисляются. Например 400 Гц при делении на 19200 получаем 48 которое не кратно 2 в степени n. Выход из положение я выбрал взяв близкое число к числу 2 в степени n. Наиболее близкое это 240 оно близко к 256. То есть из 512 мы взяли только 240 отсчетов. Кроме того нельзя просто взять просто близкое. Например мы могли взять и 480 которое близко к 512, но тем не менее взяли близкое к 256. Объяснение этому в том, что на разных частотах число отсчетов влияет на полосу пропускания. Чем больше число отсчетов тем уже полоса пропускания. Связано это с тем что на высокой частоте косинус проходит гораздо быстрее период нежели на низкой и амплитуда вычисляется на столько точно, что соседние частоты просто выбрасываются и между частотами образуются слепые зоны частот которые анализатором не воспринимаются. Для того, чтобы анализатор воспринимал все частоты и охватывал весь спектр необходимо на высоких частотах расширить полосу взяв меньше отсчетов, а на низких как можно больше сузить взяв отсчетов соответственно больше. Таким образом на путем практического подбора числа отсчетов я подобрал такие:
9,6 кГц 16 отсчетов, 4,8 кГц 32 отсчета, 2,4 кГц 32 отсчета, 1,6 кГц 60 отсчетов, 1,2 кГц 64 отсчета, 800 Гц 240 отсчетов, 600 Гц 256 отсчетов, 500 Гц 252 отсчета, 400 Гц 240 отсчетов, 300 Гц 512 отсчетов, 150 Гц 512 отсчетов, 75 Гц 512 отсчетов.
Таким выбором числа отсчетов удалось полосу равномерной по всему диапазону частот.
Еще один подводный камень получился на частоте 9,6 кГц. Так как мнимой части нет(это легко проверить подставив в формулу выше при 512 отсчетах 256 номер спектра и синус будет всегда равен 0), то реальная часть может достаточно сильно изменяться за счет того, что значение косинуса будет вычисляться через раз в противофазе к основному сигналу. То есть будет вычисляться раз. Для того, чтобы этого не было необходимо вычислить хотя бы 2 значения реальной части сдвинутой на 90 градусов и выбрать максимальный из двух значений.
Алгоритм программы накопление 512 отсчетов в промежутке перевод мк в режим сна и пробуждение когда очередной отсчет готов. Кроме того 150 Гц происходит развертка светодиодов - это раз в 128 от частоты дискретизации в 19200. То есть до того как ацп снимет все отсчеты он успеет полностью провести одну полную развертку. Как только все отсчеты готовы в основном цикле программы происходит вычисление всех амплитуд спектра. В это время развертка продолжается, но мк не впадает в сон, а считает амплитуды. Как только амплитуды посчитаны мк переводится в сон и программа повторяется заново. Амплитуды рассчитываются исходя из 20 дб диапазона, то есть прологарифмированы.
Исходя из времени на получение всех отсчетов и время на расчет всех амплитуд частота обновления находится в районе 10-15 Гц.
В статье описан небольшой анализатор аудиоспектра (0 - 10 кГц), состоящий из ЖК-дисплея 16x2 и микроконтроллера ATmega32. Используется простой алгоритм ДПФ (Дискретное Преобразование Фурье). БПФ (Быстрое Преобразование Фурье) отличается от ДПФ только большей скоростью но и более сложным алгоритмом.
ДПФ медленный по сравнению с БПФ. Данный ЖК анализатор спектра не требует большой скорости, которую может обеспечить БПФ, и если изображение на экране будет меняться с частотой около 30 кадров / сек, то это более чем достаточно для визуализации звукового спектра. Для ЖК-дисплея не рекомендуется слишком высокая частота обновления. Звук с частотой дискретизации 20 кГц даёт 32 точки ДПФ. Поскольку результат преобразования симметричен, мне нужно использовать только первые 16 результатов. Соответственно максимальная частота 10 кГц. Таким образом, 10кГц/16 = 625Гц.
Можно увеличить скорость вычисления ДПФ. Если есть точка N ДПФ, то необходимо найти синус и косинус (N ^ 2) / 2. Для 32-точечного ДПФ, необходимо найти синус и косинус 512. Прежде чем искать синус и косинус, нам нужно найти угол (градусы), который занимает некоторое время процессора. Для ускорения этого сделаны таблицы для синуса и косинуса. Синус и косинус сделаны 16-битными переменными, умножив значения синуса и косинуса на 10000. После преобразования нужно разделить каждый результат на 10000. Теперь можно рассчитать 120 32-точечных ДПФ в секунду, что более чем достаточно для анализатора спектра.
Дисплей
Для отображения столбиков использованы пользовательские символы для ЖК-дисплея, загруженные в 64 Байт встроенной памяти дисплея.
Аудио вход
Одной из наиболее важных частей анализатора спектра является получение сигнала с электретного микрофона. Особое внимание должно быть уделено разработке предварительного усилителя для микрофона. Нам нужно установить нулевой уровень на входе АЦП и максимальный уровень равный половине напряжения питания, т.е. 2,5В. На него может подаваться напряжение от -2,5В до +2,5В. Предусилитель должен быть настроен так, чтобы не превышать этих границ. В схеме использован операционный усилитель LM324 в качестве предварительного усилителя для микрофона.
В качестве устройства отображения используется двухстрочный символьный ЖК индикатор. Основным моментом при реализации данного проекта является не аппаратная часть, а программная, точнее реализация дискретного преобразования Фурье (ДПФ) на 8-разрядном микроконтроллере. Сразу следует отметить, что автор не является экспертом в этой области и поэтому начал с основ - с простого дискретного преобразования Фурье. Алгоритм быстрого преобразования Фурье является не только быстрым, но и достаточно сложным.
Дискретное преобразование Фурье (в англоязычной литературе DFT, Discrete Fourier Transform) - это одно из преобразований Фурье, широко применяемых в алгоритмах цифровой обработки сигналов (его модификации применяются в сжатии звука в MP3, сжатии изображений в JPEG и др.), а также в других областях, связанных с анализом частот в дискретном (к примеру, оцифрованном аналоговом) сигнале. Дискретное преобразование Фурье требует в качестве входа дискретную функцию. Такие функции часто создаются путем дискретизации (выборки значений из непрерывных функций).
Принципиальная схема анализатора спектра звукового сигнала очень проста и условно ее можно разделить на цифровую часть и аналоговую.
Цифровая часть образована микроконтроллером и подключенным к нему ЖК индикатором. Микроконтроллер тактируется от кварцевого резонатора 16 МГц, в качестве опорного напряжения для АЦП микроконтроллера используется напряжение питания +5 В.
Шина данных ЖК индикатора подключена к порту C микроконтроллера (линии ввода/вывода PC0-PC3), шина управления подключена к порту D(PD5, PD6) микроконтроллера. Индикатор работает в 4-битном режиме. Переменный резистор номиналом 4.7 кОм используется для регулировки контрастности. Для работы с индикатором были созданы пользовательские символы для отображения 8 горизонтальных столбиков анализатора, эти пользовательские символы занимают все 64 Байта ОЗУ ЖК индикатора.
Микроконтроллер работает от внешнего кварцевого резонатора 16 МГц.
Аналоговая часть устройства - это самая важная часть и представляет собой предварительный усилитель сигнала электретного микрофона, выход которого подключается к каналу ADC0 встроенного в микроконтроллер АЦП. Уровень нуля на входе АЦП нам необходимо установить равным точно половине опорного напряжения, т.е. 2.5 В. В этом случае мы сможем использовать положительную и отрицательную полуволну сигнала, но его амплитуда не должна превышать установленный предел, т.е. коэффициент усиления должен быть точно настроен для предотвращения перегрузки. Всем вышеуказанным условиям отвечает распространенная микросхема низкопотребляющего операционного усилителя .
Алгоритм ДПФ несколько медленнее в сравнении с быстрым преобразованием Фурье. Но наш анализатор спектра не требует высокой скорости, и если он способен обеспечить скорость обновления около 30 кадров в секунду, этого будет более чем достаточно для визуализации спектра звукового сигнала. В любом случае, в нашем варианте возможно достичь скорости 100 кадров в секунду, но это уже слишком высокое значение параметра для двухстрочного символьного ЖК индикатора и оно не рекомендуется. Частота дискретизации равна 20 кГц для 32 точечного дискретного преобразования Фурье и поскольку результат преобразования симметричен, нам нужно использовать только первую половину, т.е. первые 16 результатов. Следовательно, мы можем отображать частотный спектр в диапазоне до 10 кГц и разрешение анализатора составляет 10 кГц/16 = 625 Гц.
Автором конструкции были предприняты попытки увеличения скорости вычисления ДПФ. Если это преобразование имеет N точек, то мы должны найти N2/2 значений синуса и косинуса. Для нашего 32 точечного преобразования необходимо найти 512 значений синуса и косинуса. Но, прежде чем найти их нам необходимо вычислить угол (градусы), что займет некторое процессорное время, поэтому было решено использовать для этих вычислений таблицы значений. При расчетах в программе микроконтроллера не используются вычисления с плавающей точкой и числа двойной точности (double), так как это займет больше времени на обработку на 8-разрядном микроконтроллере. Вместо этого значения в таблицах поиска используются 16-разрядные данные целочисленного типа (integer), умноженные на 10000. Затем после выполнения преобразования результаты делятся на 10000. При таком подходе имеется возможность выполнять 120 32-точечных преобразований в секунду, что более чем достаточно для нашего устройства.
Демонстрация работы анализатора спектра на микроконтроллере ATmega32
Загрузки
Исходный код (программа микроконтроллера, таблицы данных синуса, косинуса и угла) -
- Понятно, что на АВР-ке дальше светомузыки сложно уехать, не те параметры. Но 120 32-точечных преобразований в секунду для большинства задач может быть достаточно. А выборку 625Гц, можно конечно и другую взять, по точнее потеряв частоту обновления. Стоит отметить, что МК будет себя плохо чувствовать, в плане производительности мало что на него еще навешаешь. Но тут можно же организовать выдачу результата по аппаратным методам передачи данных. Тогда это будет вспомогательный МК, а основной будет только принимать с него данный и обрабатывать совместимо с другими процессами. По большому счету все же в частоту проца упирается. Когда-то получалось разгонять мегу выше 20 МГц, но для этих задач наверно получим только глюки на высоких частотах. Идея хороша, только бы больше мат части расписано было бы... именно ее реализация на МК
- я и поинтересней анализаторы делал: You Tube или вариант на цветном ЖКИ: You Tube в основе знаменитая библиотека Чена:)
- "нам необходимо вычислить угол (градусы)" А может кто-нибудь подробнее рассказать как рассчитываются значения для этих таблиц?
- С таблицей синусов и косинусов все понятно. Не понятно как рассчитываются значения в таблице degree_lookup?