Теория
Для начала немного теории. Как известно все в подобных анализаторах используется быстрое преобразование Фурье и часто говорится, что ДПФ в подобных конструкция использовать нельзя, только БПФ да и то в на ассемблере. Я же использовал вместо этого дискретное преобразование Фурье(ДПФ) и преобразование по Уолшу. И в этой статье докажу, что можно использовать даже не только БПФ, а ДПФ написанный на С. Но сначала по порядку как из ДПФ получить простую функцию ДПФ и по Уолшу. ДПФ классически выглядит следующим образом:
Так как у мк мало ресурсов, то заменяют 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 кадров / сек, то это более чем достаточно для визуализации звукового спектра. Но я итак могу достичь частоты около 100 кадров / сек, однако для ЖК-дисплея не рекомендуется слишком высокая частота обновления. Звук с частотой дискретизации 20 кГц даёт 32 точки ДПФ. Поскольку результат преобразования симметричен, мне нужно использовать только первые 16 результатов. Соответственно максимальная частота 10 кГц. Таким образом, 10кГц/16 = 625Гц.
Я пытался увеличить скорость вычисления ДПФ. Если есть точка N ДПФ, то необходимо найти синус и косинус (N ^ 2) / 2. Для 32-точечного ДПФ, необходимо найти синус и косинус 512. Прежде чем искать синус и косинус, нам нужно найти угол (градусы), который занимает некоторое время процессора. Для этого я сделал таблицы для синуса и косинуса. Я сделал синус и косинус 16-битными переменными, умножив значения синуса и косинуса на 10000. После преобразования я должен разделить каждый результат на 10000. Теперь я могу рассчитать 120 32-точечных ДПФ в секунду, что более чем достаточно для анализатора спектра.
Дисплей
Я использовал пользовательские символы для ЖК-дисплея загруженные в 64 Байт встроенной памяти ЖК-дисплея. В интернете я увидел видео, где ЖК-дисплей 16х2 используется в качестве дисплея анализатора спектра и использовал эту идею.
Аудио вход
Одной из наиболее важных частей анализатора спектра является получение сигнала с электретного микрофона. Особое внимание должно быть уделено разработке предварительного усилителя для микрофона. Нам нужно установить нулевой уровень на входе АЦП и максимальный уровень равный половине напряжения питания, т.е. 2,5В. На него может подаваться напряжение от -2,5В до +2,5В. Предусилитель должен быть настроен так, чтобы не превышать этих границ. Я использовал операционный усилитель LM324 в качестве предварительного усилителя для микрофона.
Список радиоэлементов
| Обозначение | Тип | Номинал | Количество | Примечание | Магазин | Мой блокнот | |
|---|---|---|---|---|---|---|---|
| Дисплей | |||||||
| МК AVR 8-бит | ATmega32 | 1 | В блокнот | ||||
| Конденсатор | 22 пФ | 2 | В блокнот | ||||
| Конденсатор | 0.1 мкФ | 1 | В блокнот | ||||
| Электролитический конденсатор | 100 мкФ | 1 | В блокнот | ||||
| Резистор | 100 Ом | 1 | В блокнот | ||||
| Подстроечный резистор | 4.7 кОм | 1 | В блокнот | ||||
| Кварцевый резонатор | 16 МГц | 1 | В блокнот | ||||
| LCD-дисплей | 16х2 | 1 | В блокнот | ||||
| Блок питания | 5 В | 1 | В блокнот | ||||
| Аудио вход | |||||||
| U1 | Операционный усилитель | LM324 | 1 | В блокнот | |||
| С1 | Конденсатор | 1 мкФ | 1 | В блокнот | |||
| С8 | Конденсатор | 0.01 мкФ | 1 | В блокнот | |||
| R1 | Резистор | 220 кОм | 1 | В блокнот | |||
| R2, R3 | Резистор | 10 кОм | 2 | В блокнот | |||
| R4, R9 | Резистор | 1 кОм | 2 | В блокнот | |||
| R5 | Резистор | ||||||
Для цифровой обработки сигналов (DSP) существует множество специализированных процессоров, это DSP от Texas Instruments серии TMS320, которая включает в себя и простые целочисленные ядра, и такие монстры как подсемейство семейство C6000, обрабатывающие данные с плавающей точкой. Есть целая серия ADSP от Analog Devices (в которую входит более-менее универсальный BlackFin), есть и более простые решения от MicroChip - dsPIC.
Однако, специализированный DSP - это хорошо, но вот всегда ли он так необходим? Да, при огромном потоке информации он просто незаменим, однако есть и более простые задачи обработки. Конкретно меня интересовала задача двойного преобразования - делается свертка аудиосигнала, тем самым получается спектр, затем над спектром можно провести любые операции и выполнить обратное преобразование, получив тем самым обработанный сигнал. Все это нужно провести в реальном времени и получить качество не ниже телефонного.
Сейчас не 2000-ый год, существуют и существенно упали в цене однокристальные решения на производительных ядрах ARM7/Cortex-M3, они являются 32 битными, имеют аппаратную реализацию операции 32-битного умножения (мало того, почти DSP-операцию умножения с накоплением и 64 битным результатом), а Cortex-M3 включает еще и аппаратное деление.
Сразу хочу предупредить, обработка сигналов не является моим профилем, почти все знания (а точнее понимания принципов) сохранились еще с института, сейчас же захотелось просто проверить и реализовать их. Это я к тому, что могут встречаться неточности в описании, подмены понятий и прочее. Собственно, академическая точность меня не очень сильно волновала.
Почти для любого DSP, озвученная выше задача является простой и подъемной. Но как же поведет себя на ней RISC-ядро общего применения? Если рассматривать AVR или PIC, то их вряд ли хватит. Сказывается 8-битность и низкая тактовая частота. Хотя, у Elm-Chan -а есть конструкции, в которых он на AVR проводит БПФ и рисует спектр сигнала. Однако, в данном случае в реальном времени делается просто визуализация (с минимальной точностью обработки), а не полная обработка сигнала с приемлемым аудио-качеством.
В качестве подопытного чипа был выбран LPC2146, основанный на ядре ARM7TDMI-S и имеющий максимальную тактовую частоту 60МГц (на практике, на работает на 72 и даже на 84МГц). Почему? Во-первых, у меня есть отладочная плата для него, во-вторых, на борту имеется АЦП и ЦАП, т.е. требуется минимальная внешняя обвеска.
Немного теории
В первую очередь, интересно было оценить производительность БПФ (Быстрого Преобразования Фурье) на ARM-микроконтроллерах. По этой оценке можно сделать вывод: хватит ли у него скорости для обработки потока аудиоданных, и сигнал с какой частотой дискретизации и на сколько каналов можно будет обработать на таком микроконтроллере.
На основе преобразования Фурье можно сооружать хитрые фильтры (при этом с очень заманчивыми характеристиками). Меня же в первую очередь интересовали задачи изменения тональности сигнала (поднятие и опускание спектра) и "отражение" спектра. Последнее требуется в SDR-радиоприемниках для прослушивания радиопередачи в нижней боковой полосе LSB .
Не буду загружать теорией и объяснять что такое Преобразование Фурье, материалов на эту тему довольно много, из того, что использовал я: вики и глава из очень хорошей и информативной книги .
Программная реализация
Программных реализаций БПФ существует очень много, однако, я писал собственную. Основная цель, которую я преследовал: оптимизация кода на конкретную архитектуру. Во-первых, я сразу ориентировался на 32-битность, во-вторых требовались только целочисленные вычисления и желательно было избежать операции деления. Под эти требования подобрать нечто готовое уже проблемно.
Все константы, которые можно было рассчитать заранее, были рассчитаны и помещены в таблицы (в основном это значения тригонометрических функций). Это основная оптимизация алгоритма, в остальном он почти полностью повторяет описанный алгоритм .
Самым же важным является требование к целочисленным вычислениям. В процессе реализации, была даже допущена ошибка, из-за которой происходило переполнение в одной из 32 битных переменных цикла. Мало того, она возникала не на всех тестовых данных, поэтому доставляла приличную головную боль, пока не была найдена.
Все исходные тексты я собрал в одном архиве . Реализация не окончательная, есть дублирующие вычисления (не учитывается симметрия спектра и фаз), и требуется оптимизация использования буферов, так как в текущий момент для расчетов используется слишком много оперативной памяти (почти 12к для преобразования из 1024 точек).
Первые тесты
Барабанная дробь: тестирую скорость работы преобразования для выборки из 1024 точек, частота процессорного ядра 60МГц. Тестирование проводилось в эмуляторе, так что не претендует на 100% точность, однако этот результат можно использовать в качестве показателя (по моему предыдущему опыту эмулятор хоть и врал, но не сильно). Тест первой версии кода, компилятор RealView MDK, опция оптимизации O3, ARM-режим генерации кода.
Итак, что мы видим:
По 6мс на каждое преобразование, итого чуть более 12мс на преобразование "туда и обратно". Получается, что при частоте дискретизации 44100Гц (стандартная для аудио) и выборках разрешением до 16 бит, чистые вычисления будут занимать ~44*12мс = 528мс. И это на промежуточной версии микропрограммы, когда еще не закончены некоторые оптимизации кода (по прикидкам алгоритм может быть ускорен еще чуть ли не в 2 раза)! По-моему просто отличный показатель.
Итого, загрузка ядра предполагается в районе 50%, еще 50% остается на преобразования над спектром и накладные расходы при работе с АЦП, ЦАП и прочие передачи данных. Если же понизить частоту выборок до "телефонного" уровня (около 4800-9600Гц), то загрузка ядра будет еще ниже (порядка 15-30%).
Итак, с математической частью более-менее понятно. Можно приступить к конкретной реализации.
Железо
Для тестовой платформы взята отладочная плата Keil MCB2140 с динамиком. Напаян шнурок Mini-Jack для подключения к линейному выходу устройства и собрана простейшая входная цепочка. Как уже было сказано, на плате уже установлен динамик, подключенный к аналоговому выходу микроконтроллера и есть элементы управления (кнопка и потенциометр).
Вот набросок входной цепи:
Отладка софта происходила поэтапно:
- Отладка всей необходимой периферии: АЦП, ЦАП, таймеры, светодиодная индикация.
- Тест с оцифровкой сигнала: оцифровываю данные с требуемой скоростью и кладу в буфер, затем извлекаю данные из буфера и воспроизвожу сигнал. Т.е. простой сдвиг сигнала во времени, без каких-либо преобразований. На этом этапе тестируется механизм работы с 2 буферами, необходимый для дальнейшей работы.
- К предыдущей версии добавляется прямое и обратное преобразование Фурье. Этот тест окончательно проверяет корректность работы кода БПФ, а так же проверка, что код укладывается в доступную производительность.
- После этого основной костяк приложения сделан, можно перейти к практическим применениям.
Проблема возникла после добавления в код БПФ: в сигнале появились посторонние шумы и свисты. Вообще, данное поведение показалось мне довольно странным, т.к. без преобразования сигнал проходя цифровой тракт оставался достаточно чистым. Первая причина этого: после аналоговой цепи, амплитуда сигнала на АЦП была не полной (0-3.3В), а только в пределах 0.5-2В на максимальной громкости плеера, вторая: достаточно сильный шум из-за целочисленных вычислений (+-1 единица, что оказалось достаточным для появления слышимых помех).
Для борьбы с первой проблемой, было решено заняться подстройкой аналоговой части. А для решения вопроса с шумами - попробовать использовать lowpass-фильтр.
Прикладное применение 1: изменяем тональность сигнала
На плате предусмотрен потенциометр (переменные резистор), его можно использовать для управления. В данном случае он задает смещение спектра сигнала вверх-вниз, вполне достаточно чтобы "преобразить" любимые композиции.
Вот что происходит в частотной области:
При этом результат преобразования содержится в 2 буферах. Один буфер - действительная часть, а другой - мнимая. Физический смысл полученных в них чисел: в действительной части содержатся значения гармоник, в мнимой - сдвиг фаз для данных гармоник. Мало того, как видно, изначальный сигнал описывается N-значениями, а после преобразования получается 2N-значений. Количество информации не меняется, а увеличение в 2 раза количества информации происходит из-за того, что данные буфера имеют избыточность в виде дублирования значений.
В предыдущей мы подключали дешевый китайский LCD экран к плате STM32L4 Discovery . Теперь мы попробуем реализовать на этой комбинации что-то выходящее за рамки традиционного моргания светодиодом, а именно анализатор звукового спектра, который использует имеющийся на плате микрофон. Заодно я расскажу, как пользоваться операционной системой FreeRTOS, и зачем она нужна, а также почему в нотной октаве 12 нот, и чем 53 ноты лучше, чем 12.
Оцифровка звука
Мы хотим получать сигнал с микрофона, вычислять его спектр с помощью быстрого преобразования Фурье (FPU нам в помощь) и показывать результат на LCD в виде "цветного водопада". Силу звука будем кодировать цветом. Будем рисовать с краю дисплея строку пикселей, где самый левый пиксель будет соответствовать минимальной частоте, а самый правый - максимальной, при этом предыдущая картинка будет смещаться на одну строку, освобождая место для новой строки. Наш микроконтроллер слишком сложен, чтобы начать с нуля, поэтому начнем с примера из комплекта STM32Cube, который называется DFSDM_AudioRecord. Что такой DFSDM? Это Digital Filter for Sigma-Delta Modulation. Дело в том, что в отличие от старых добрых аналоговых микрофонов, тот, что стоит на плате Discovery, выдает сигнал не в виде напряжения, пропорционального звуковому давлению, а в виде последовательности нулей и единиц с тактовой частотой в несколько мегагерц. Если пропустить эту последовательность через фильтр низких частот, то получится тот самый аналоговый сигнал. В предыдущих моделях микроконтроллеров приходилось делать цифровой фильтр, чтобы получить звуковой сигнал в цифровом виде. Теперь в микроконтроллере есть специальный модуль для этого, и все, что требуется, - это настроить его на старте программы. Для этого можно или углубиться в чтение документации, или воспользоваться готовым примером. Я пошел по второму пути. Следующая картинка иллюстрирует внутреннюю структуру программы DFSDM_AudioRecord.
Оцифрованный звук с помощью DMA попадает в кольцевой буфер. DMA вызывает прерывание дважды: один раз - когда буфер заполнен наполовину, второй раз - когда он заполнен полностью. Процедура обработки прерываний просто выставляет соответствующий флажок. Функция main() после инициализации исполняет бесконечный цикл, где проверяются эти флажки и, если флажок выставлен, копируется соответствующая половина буфера. Пример копирует данные в другой буфер, откуда они, опять-таки с помощью DMA, отправляются на усилитель наушников. Я оставил эту функциональность, добавив вычисление спектра звукового сигнала.
Когда задач много
Прямолинейный способ добавить новую функциональность в наш код - добавить еще флажков и написать функции, которые будут вызываться, если эти флажки выставлены. В результате обычно получается каша из флажков, функций-обработчиков и глобального контекста, который вынужден быть глобальным, поскольку решение одной задачи разбивается на множество мелких шагов, реализованных отдельными функциями - обработчиками событий. Альтернативный способ - поручить управление задачами операционной системе, например FreeRTOS. Это позволяет значительно упростить логику за счет того, что каждая задача решается в рамках своего цикла обработки событий, которые взаимодействуют друг с другом посредством функций операционной системы. Например, мы можем добавить задачу обработки данных в виде отдельного цикла, который будет ждать готовности данных на синхронизационном примитиве - семафоре. Семафор устроен очень просто: вы можете пройти его, если флажок поднят, при этом флажок автоматически опускается. Поднимет флажок в нашем случае источник данных, когда подготовит данные для другой задачи. Подобным образом можно создавать произвольные цепочки из задач-источников данных и задач-потребителей данных подобно тому, как это происходит, например, в операционной системе линукс.
Конечно, одновременность исполнения задач - это иллюзия, особенно, когда вычислительное ядро всего одно. В этом случае мы можем говорить о том, что у нас есть единственный поток исполнения программы процессором. Семафоры, как и другие синхронизационные примитивы, играют роль волшебной кроличьей норы, в которую проваливается поток исполнения, чтобы вынырнуть в другой задаче.
Подключить FreeRTOS к своему проекту достаточно просто. Нужно лишь заменить бесконечный цикл, которым обычно заканчивается функция main() в микроконтроллере, на вызов osKernelStart(). После этого компилятор объяснит вам, чего именно ему не хватает для компиляции. Все действия, которые вы до этого выполняли в цикле, нужно перенести в отдельную задачу и зарегистрировать ее с помощью вызова xTaskCreate. После этого вы сможете добавить еще столько задач, сколько захотите. Нужно иметь ввиду, что между вызовами xTaskCreate и osKernelStart лучше не размещать никакого кода, работающего с железом, поскольку здесь системный таймер может работать неправильно. Вызов обработчика таймера операционной системы osSystickHandler() нужно добавить в SysTick_Handler(), а две функции SVC_Handler и PendSV_Handler убрать из своего кода, поскольку они реализованы в коде ОС. При регистрации задач важно не ошибиться с размером стека. Если он окажется слишком мал, вы получите краши в самых неожиданных местах. Первым при переполнении стека страдает сама структура, описывающая задачу. В IAR есть возможность посмотреть список задач. Если вы видите в нем задачу с измененным именем, значит нужно увеличить размер стека.
Вычисляем спектр
Для вычисления спектра мы воспользуемся быстрым преобразованием Фурье. Соответствующая функция уже есть в библиотеке. Она получает буфер, заполненный комплексными данными, и формирует результат там же. Соответственно, на входе ей нужен буфер, где оцифрованный звук чередуется с нулями (комплексная часть 0). На выходе мы получаем комплексные числа, для которых сразу вычисляем квадрат модуля, сложив квадраты действительной и мнимой части. Мы делаем это только для половины буфера, поскольку спектр симметричен. Вторая половина нам понадобилась бы, если бы мы захотели сделать обратное преобразование, но для простого показа спектра она не нужна. Некоторые дополнительные усилия необходимы для того, чтобы иметь возможность вычислять спектр в разных спектральных диапазонах. Чтобы получить спектр для низких частот, я аккумулирую данные за несколько циклов чтения буфера, эффективно снижая частоту дискретизации звука, которая изначально составляет 44.1kHz. В итоге получается 6 диапазонов - 20kHz, 10kHz, 5kHz, 2600Hz, 1300Hz, 650Hz. Для переключения диапазонов используется джойстик и отдельная задача. Джойстик также выполняет функции запуска / останова "водопада", а также регулировки чувствительности. Показывать спектр удобнее в логарифмических единицах (децибелах), поскольку его динамический диапазон обычно весьма велик, и в линейном масштабе мы сможем различить лишь самые сильные составляющие спектра. Логарифм считается довольно долго даже на FPU, поэтому я заменил реальный логарифм кусочно-линейной аппроксимацией, которую легко получить, зная формат представления числа в float32 . Старший бит - это знак. Следующие 8 бит - двоичная экспонента плюс 127. Оставшиеся биты - это дробная часть мантиссы при том, что целая часть равна 1 (нюансы денормализованных чисел для простоты опустим). Значит, выделив из float32 экспоненту и прихватив несколько старших бит мантиссы, можно получить неплохую аппроксимацию логарифма. Полученное число мы с помощью предварительно заготовленной таблицы преобразуем в RGB код для показа на LCD. Получается цветовая шкала на 90 или 60 децибел. Уровень громкости, соответствующий нулю этой шкалы, можно настраивать, нажимая джойстик вверх и вниз.Выводим картинку - о пользе чтения даташитов
Теперь нам осталось вывести картинку и оживить наш "водопад". Прямолинейный способ сделать это - хранить картинку со всего экрана в буфере, обновлять ее там и перерисовывать каждый раз, когда появляются новые данные. Мало того, что это решение крайне неэффективное, у нас еще и недостаточно памяти, чтобы хранить всю картинку. Казалось бы, у самой LCD достаточно памяти для этого, и она должна уметь делать с ней что-то интересное. Действительно, изучение даташита позволило обнаружить доселе никем не использованную команду скроллинга, которая позволяет динамически менять способ отображения памяти контроллера LCD на экран. Представим себе, что память - это замкнутая в кольцо лента, которую вы видите под стеклом экрана. Команда Vertical Scrolling Start Address (0x37) позволяет задать позицию на ленте, соответствующую верхнему краю экрана. Значит, все, что нам нужно, чтобы оживить "водопад" - это записать в эту позицию новый спектр и прокрутить ленту памяти. Соответствующий код был добавлен в драйвер LCD, позаимствованный у уважаемого Peter Drescher , и адаптированный как описано . Единственный недостаток подобного подхода: скроллинг работает только вдоль длинной стороны экрана. Соответственно, для вывода спектра доступна только короткая сторона.Почему в октаве 12 нот?
Перейдем к практическим применениям нашего устройства. Первое, что легко увидеть на спектре, это гармоники, то есть частоты, кратные частоте основного тона. Особенно много их в голосе. Есть они и в звуках, которые издают музыкальные инструменты. Легко понять, почему ноты соседних октав различаются по частоте в 2 раза: тогда ноты более высокой октавы совпадают по частоте со второй гармоникой нот низкой октавы. Говорят, что при этом они звучат «в унисон». Чуть сложнее разобраться в том, почему в октаве 12 нот - семь основных (белые клавиши на клавиатуре фортепьяно) плюс 5 дополнительных (черные клавиши). Дополнительные ноты обозначаются через основные с диезными и бемольными знаками, хотя по сути никакой разницы между ними и основными нотами нет - все 12 нот образуют геометрическую прогрессию так, что отношение частот между соседними нотами равно корню 12-й степени из 2. Смысл такого деления октавы на ноты в том, чтобы для любой ноты нашлись другие ноты, отличающиеся от нее по частоте в полтора раза - такая комбинация называется квинтой. Ноты, образующие квинту, звучат в унисон потому, что вторая гармоника одной ноты совпадает по частоте с третьей гармоникой другой ноты. На фото ниже показаны спектры нот До и Соль, образующих квинту, совпадающие гармоники обведены желтым.
Как же получилось, что нот 12? Поскольку ноты образуют геометрическую прогрессию, перейдем к логарифмам. ln(1.5)/ln(2) = 0.58496… Близкое значение получается у дроби 7/12 = 0.583… То есть, семь полутонов (интервалов между соседними нотами) оказываются весьма близки к квинте - 1.498. Интересно, что гораздо большую точность дает дробь 31/53 = 0.58491.., так что квинта отличается от 1.5 только в пятом знаке после запятой. Этот факт не остался незамеченным, но музыкальные инструменты с 53 нотами в октаве не получили распространения. Их сложно настраивать, на них сложно играть, а процент людей, способных почувствовать разницу с обычными инструментами, исчезающе мал.
Введение
Книги и публикации по цифровой обработке сигналов пишут авторы зачастую не догадывающиеся и не понимающие задач, стоящих перед разработчиками. Особенно это касается систем, работающих в реальном времени. Эти авторы отводят себе скромную роль бога, существующего вне времени и пространства, что вызывает некоторое недоумение у читателей подобной литературы. Данная публикация имеет целью развеять недоумения, возникающие у большинства разработчиков, и помочь им преодолеть «порог вхождения», для этих целей в тексте сознательно используется аналогии и терминология сферы программирования.
Данный опус не претендует на полноту и связность изложения.
Добавлено после прочтения комментариев.
Публикаций о том как делать БПФ немеряно, а о том как сделать БПФ, преобразовать спектр, и собрать сигнал заново, да еще и в реальном времени, явно не хватает. Автор пытается восполнить этот пробел.
Часть первая, обзорная
Существуют два основных способа построения дискретных линейных динамических систем. В литературе, такие системы принято называть цифровыми фильтрами, которые подразделяются на два основных типа: фильтры с конечной импульсной характеристикой (КИХ) и фильтры с бесконечной импульсной характеристикой (БИХ).
Алгоритмическая сущность фильтра с КИХ заключается в дискретном вычислении интеграла свертки:
Где x(t) – входной сигнал
y(t) – выходной сигнал
h(t) – импульсная характеристика фильтра или реакция фильтра на дельта функцию. Импульсная характеристика является обратным преобразованием Фурье комплексной частотной характеристики фильтра K(f).
Для формирования ясной картины у читателя, приведем пример дискретного вычисления интеграла свертки на языке С в реальном времени.
#define L (4) //длинна фильтра
int FIR(int a)
{
static int i=0; //текущая позиция
static int reg[L]; //массив входных значений
static const int h[L]={1,1,1,1};//импульсная характеристика
int b=0;//выходное значение
reg[i]=a; //копируем входное значение в массив входных значений
for(int j=0;j Вызывая данную функцию через определенные интервалы времени T и передавая ей в качестве аргумента входной сигнал, на выходе мы получим выходной сигнал, соответствующий реакции фильтра с импульсной характеристикой вида: H(t)=1 при 0 Всем собравшимся фильтр с такой импульсной характеристикой более известен под названием «фильтр скользящего среднего», и, соответственно, реализуется он гораздо проще. В данном случае такая импульсная характеристика используется для примера. Синтезу импульсных характеристик КИХ фильтров посвящена масса литературы, также имеются готовые программные продукты для получения фильтров с заданными свойствами. Автор предпочитает глючный инструмент Filter Design из пакета Matlab, но это дело вкуса. Используя фильтр с конечной импульсной характеристикой, удается немножечко воспарить над привычной реальностью, так как, в природе не существует динамических систем, имеющих конечную импульсную характеристику. Фильтр КИХ - попытка зайти в частотно-временную область с другого конца, не так как ходит природа, поэтому частотные характеристики этих фильтров зачастую обладают неожиданными свойствами. Намного ближе к природе фильтры с бесконечной импульсной характеристикой. Алгоритмическая сущность фильтров с бесконечной импульсной характеристикой сводится к рекуррентному (не путать с рекурсивным!) решению дифференциального уравнения, описывающего фильтр. То есть, каждое последующие значение выходного сигнала фильтра вычисляется на основании предыдущего значения. Именно так протекают процессы в реальном мире. Камень, падая с небоскреба каждую секунду, увеличивает свою скорость на 9.8м/с Speed=Speed+9.8, и пройденный путь каждую секунду увеличивается Distance=Distance+Speed. Кто скажет, что это не рекуррентный алгоритм, пусть первый бросит в меня камень. Только в нашей Матрице
временной интервал вызова функции возвращающей положение камня много меньше цены деления доступных нам средств измерения. Отдельно хотелось бы определить понятие «порядок фильтра». Это количество переменных которые подвергаются рекуррентным операциям. В приведенном примере функция, возвращающая скорость камня - первого порядка, функция, возвращающая пройденный путь - второго порядка. Для окончательного просветления читателя приведем пример на языке С самого простого фильтра низких частот, широко известного как фильтр «фильтр экспоненциального сглаживания» #define alfa (2) //параметр сглаживания
int filter(int a)
{
static int out_alfa=0;
out_alfa=out_alfa - (out_alfa >>alfa) + a;
return (out_alfa >> alfa);
}
Вызывая данную функцию с частотой F и передавая ей в качестве аргумента входной сигнал, на выходе мы получим выходной сигнал, соответствующий реакции фильтра низких частот первого порядка с частотой среза: Приведенный пример исходного кода совершенно неудобоварим с точки зрения понимания сути алгоритма. С точки зрения рекуррентной сути (смотри «падение камня») алгоритма, правильнее y=y+((x-y)>>alfa);, но в этом случае происходит потеря alfa значащих разрядов. Рекуррентное выражение фильтра, из примера кода, построено таким образом, чтобы избежать потери значащих разрядов. Именно конечная точность вычислений может испортить всю прелесть цифрового фильтра с бесконечной импульсной характеристикой. Особенно это заметно на фильтрах высоких порядков, отличающихся высокой добротностью. В реальных динамических системах такая проблема не возникает, наша Матрица
производит вычисления с невероятной для нас точностью. Синтезу подобных фильтров посвящена масса литературы, также имеются готовые программные продукты (см. выше). Часть вторая. Фурье – фильтр
Из вузовских курсов (у вашего покорного слуги это был курс ОТЭЦ) многие собравшие помнят два основных подхода к анализу линейных динамических систем: анализ во временной области и анализ в частотной области. Анализ во временной области - это решение дифференциальных уравнений, интегралы свертки и Дюамеля. Эти методы анализа дискретно воплотились в цифровых фильтрах БИХ и КИХ. Но существует частотный подход к анализу линейных динамических систем. Иногда его называют операторным. В качестве операторов используются преобразование Фурье, Лапласа и т.п. Далее мы будем говорить только о преобразовании Фурье. Данный метод анализа не получил широкого распространения при построении цифровых фильтров. Автору не удалось найти вменяемых практических рекомендаций по построению подобных фильтров на русском языке. Единственное краткое упоминание такого фильтра в практической литературе [Рабинер Л., Гоулд Б., Теория и применение цифровой обработки сигналов 1978], но в данной книге рассмотрение подобного фильтра очень поверхностно. В указанной книге данная схема построения фильтра называется: «свертка в реальном времени методом БПФ», что, по моему скромному мнению, совершенно не отражает сути, название должно быть коротким, иначе времени на отдых не останется. Реакция линейной динамической системы есть обратное преобразование Фурье от произведения изображения по Фурье входного сигнала x(t) на комплексный коэффициент передачи K(f): В практическом плане, данное аналитическое выражение предполагает следующий порядок действий: берем преобразование Фурье от входного сигнала, умножаем результат на комплексный коэффициент передачи, выполняем обратное преобразование Фурье, результатом которого является выходной сигнал. В реальном дискретном времени такой порядок действий выполнить невозможно. Как брать интеграл по времени от минус до плюс бесконечности?! Его можно взять только находясь вне времени… В дискретном мире для выполнения преобразования Фурье существует инструмент - алгоритм быстрого преобразования Фурье (БПФ). Именно его мы и будем использовать при реализации нашего Фурье-фильтра. Аргументом функции БПФ является массив временных отсчетов из 2^n элементов, результатом два массива длинной 2^n элементов соответствующие действительной и мнимой части преобразования Фурье. Дискретной особенностью алгоритма БПФ является то, что входной сигнал считается периодичным с интервалом 2^n. Это накладывает некоторые ограничения на алгоритм Фурье-фильтра. Если взять последовательность выборок входного сигнала, провести от них БПФ, умножить результат БПФ на комплексный коэффициент передачи фильтра и выполнить обратное преобразование …ничего получится!
Выходной сигнал будет иметь огромные нелинейные искажения в окрестности стыков выборок. Для решения этой проблемы необходимо применить два приема: Такие функции широко применяются в технике цифровой обработки сигналов, и называть их принято - окнами. По скромному мнению автора лучшим, с практической точки зрения, является окно имени Хана: На рисунке приведены графики иллюстрирующие свойства окна Хана длинной 2^n=256. Экземпляры окна построены с половинным перекрытием k=128. Как видно все оговоренные выше свойства имеются в наличии. По просьбам трудящихся, на следующем рисунке приведена схема вычислений Фурье-фильтра, при длине выборки 2^n=8, количество выборок 3. На подобных рисунках очень сложно отобразить процесс вычислений, особенно тяжело показать его цикличность, поэтому мы и ограничились количеством выборок равным трем. Входной сигнал разбивается на блоки длинной 2^n=8 с перекрытием 50%, от каждого блока берется БПФ, результаты БПФ подвергаются нужной трансформации, берется обратное БПФ, результат обратного БПФ скалярно умножается на окно, после умножения блоки складываются с перекрытием. При выполнение трансформаций спектра, не стоит забывать о главном свойстве массива БПФ действительных сигналов, первая половина массива БПФ комплексно сопряжена со второй половиной, т.е Re[i]=Re[(1< Теперь мы знаем все, что необходимо для написания алгоритма Фурье-фильтра. Опишем алгоритм на языке С. #include Вызывая функцию FureFilter() с частотой FSempl и передавая ей в качестве аргумента входной сигнал, результатом получим выходной сигнал. В данном примере входной сигнал обрабатывается следующим образом: сигнал пропускается через полосовой фильтр с частотами среза FsrLow, FsrHi, подавляются все спектральные составляющие выше и ниже указанных частот, сдвигается спектр сигнала (для звуковых сигналов это воспринимается как эффект Буратино-Карабаса), амплитудный спектр сигнала подвергается сглаживанию фильтром низких частот (для звука это эффект гулкого помещения). Данные действия с сигналом выполнены в качестве примера, для того чтобы показать технические приемы обработки сигнала в частотной области, такие как: соблюдение комплексно-сопряженности коэффициентов, восстановление комплексного спектра по амплитудному, не используя тригонометрических функций и т.п. Заключение
Стоит отметить, что, скорее всего, данная функция Фурье-фильтра, на практике окажется неработоспособна. При вызове данной функции даже с невысокой частотой 8000Гц, она не успеет выполнится к моменту следующего вызова, не хватить быстродействия. Данный программный код Фурье-фильтра приведен в качестве описания алгоритма, без привязки к конкретным аппаратным ресурсам, и имеет чисто образовательные цели
(см. Введение). При практической реализации следует распараллелить выполнение функции заполнения-опорожнения буфера BufInOut (лучше сразу ПДП и т. п.) и функции обработки буфера ObrBuf(), но это уже совсем другая история.
