Теория
Для начала немного теории. Как известно все в подобных анализаторах используется быстрое преобразование Фурье и часто говорится, что ДПФ в подобных конструкция использовать нельзя, только БПФ да и то в на ассемблере. Я же использовал вместо этого дискретное преобразование Фурье(ДПФ) и преобразование по Уолшу. И в этой статье докажу, что можно использовать даже не только БПФ, а ДПФ написанный на С. Но сначала по порядку как из ДПФ получить простую функцию ДПФ и по Уолшу. ДПФ классически выглядит следующим образом:
Так как у мк мало ресурсов, то заменяют 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 Гц.
Во многих случаях задача получения (вычисления) спектра сигнала выглядит следующим образом. Имеется АЦП, который с частотой дискретизации Fd преобразует непрерывный сигнал, поступающий на его вход в течение времени Т, в цифровые отсчеты - N штук. Далее массив отсчетов подается в некую программку, которая выдает N/2 каких-то числовых значений (программист, который утянул из инета
написал программку, уверяет, что она делает преобразование Фурье).Чтобы проверить, правильно ли работает программа, сформируем массив отсчетов как сумму двух синусоид sin(10*2*pi*x)+0,5*sin(5*2*pi*x) и подсунем программке. Программа нарисовала следующее:
рис.1 График временной функции сигнала
рис.2 График спектра сигнала
На графике спектра имеется две палки (гармоники) 5 Гц с амплитудой 0.5 В и 10 Гц - с амплитудой 1 В, все как в формуле исходного сигнала. Все отлично, программист молодец! Программа работает правильно.
Это значит, что если мы подадим на вход АЦП реальный сигнал из смеси двух синусоид, то мы получим аналогичный спектр, состоящий из двух гармоник.
Итого, наш реальный
измеренный сигнал, длительностью 5 сек
, оцифрованный АЦП, то есть представленный дискретными
отсчетами, имеет дискретный непериодический
спектр.
С математической точки зрения - сколько ошибок в этой фразе?
Теперь начальство решило мы решили, что 5 секунд - это слишком долго, давай измерять сигнал за 0.5 сек.
рис.3 График функции sin(10*2*pi*x)+0,5*sin(5*2*pi*x) на периоде измерения 0.5 сек
рис.4 Спектр функции
Что-то как бы не то! Гармоника 10 Гц рисуется нормально, а вместо палки на 5 Гц появилось несколько каких-то непонятных гармоник. Смотрим в интернетах, что да как…
Во, говорят, что в конец выборки надо добавить нули и спектр будет рисоваться нормальный.
рис.5 Добили нулей до 5 сек
рис.6 Получили спектр
Все равно не то, что было на 5 секундах. Придется разбираться с теорией. Идем в Википедию
- источник знаний.
2. Непрерывная функция и представление её рядом Фурье
Математически наш сигнал длительностью T секунд является некоторой функцией f(x), заданной на отрезке {0, T} (X в данном случае - время). Такую функцию всегда можно представить в виде суммы гармонических функций (синусоид или косинусоид) вида:
(1), где:
K - номер тригонометрической функции (номер гармонической составляющей, номер гармоники)
T - отрезок, где функция определена (длительность сигнала)
Ak - амплитуда k-ой гармонической составляющей,
θk- начальная фаза k-ой гармонической составляющей
Что значит «представить функцию в виде суммы ряда»? Это значит, что, сложив в каждой точке значения гармонических составляющих ряда Фурье, мы получим значение нашей функции в этой точке.
(Более строго, среднеквадратичное отклонение ряда от функции f(x) будет стремиться к нулю, но несмотря на среднеквадратичную сходимость, ряд Фурье функции, вообще говоря, не обязан сходиться к ней поточечно. См. https://ru.wikipedia.org/wiki/Ряд_Фурье .)
Этот ряд может быть также записан в виде:
(2),
где , k-я комплексная амплитуда.
Связь между коэффициентами (1) и (3) выражается следующими формулами:
Отметим, что все эти три представления ряда Фурье совершенно равнозначны. Иногда при работе с рядами Фурье бывает удобнее использовать вместо синусов и косинусов экспоненты мнимого аргумента, то есть использовать преобразование Фурье в комплексной форме. Но нам удобно использовать формулу (1), где ряд Фурье представлен в виде суммы косинусоид с соответствующими амплитудами и фазами. В любом случае неправильно говорить, что результатом преобразования Фурье действительного сигнала будут комплексные амплитуды гармоник. Как правильно говорится в Вики «Преобразование Фурье (ℱ) - операция, сопоставляющая одной функции вещественной переменной другую функцию, также вещественной переменной.»
Итого:
Математической основой спектрального анализа сигналов является преобразование Фурье.
Преобразование Фурье позволяет представить непрерывную функцию f(x) (сигнал), определенную на отрезке {0, T} в виде суммы бесконечного числа (бесконечного ряда) тригонометрических функций (синусоид и\или косинусоид) с определёнными амплитудами и фазами, также рассматриваемых на отрезке {0, T}. Такой ряд называется рядом Фурье.
Отметим еще некоторые моменты, понимание которых требуется для правильного применения преобразования Фурье к анализу сигналов. Если рассмотреть ряд Фурье (сумму синусоид) на всей оси Х, то можно увидеть, что вне отрезка {0, T} функция представленная рядом Фурье будет будет периодически повторять нашу функцию.
Например, на графике рис.7 исходная функция определена на отрезке {-T\2, +T\2}, а ряд Фурье представляет периодическую функцию, определенную на всей оси х.
Это происходит потому, что синусоиды сами являются периодическими функциями, соответственно и их сумма будет периодической функцией.
рис.7 Представление непериодической исходной функции рядом Фурье
Таким образом:
Наша исходная функция - непрерывная, непериодическая, определена на некотором отрезке длиной T.
Спектр этой функции - дискретный, то есть представлен в виде бесконечного ряда гармонических составляющих - ряда Фурье.
По факту, рядом Фурье определяется некоторая периодическая функция, совпадающая с нашей на отрезке {0, T}, но для нас эта периодичность не существенна.
Периоды гармонических составляющих кратны величине отрезка {0, T}, на котором определена исходная функция f(x). Другими словами, периоды гармоник кратны длительности измерения сигнала. Например, период первой гармоники ряда Фурье равен интервалу Т, на котором определена функция f(x). Период второй гармоники ряда Фурье равен интервалу Т/2. И так далее (см. рис. 8).
рис.8 Периоды (частоты) гармонических составляющих ряда Фурье (здесь Т=2π)
Соответственно, частоты гармонических составляющих кратны величине 1/Т. То есть частоты гармонических составляющих Fk равны Fk= к\Т, где к пробегает значения от 0 до ∞, например к=0 F0=0; к=1 F1=1\T; к=2 F2=2\T; к=3 F3=3\T;… Fk= к\Т (при нулевой частоте - постоянная составляющая).
Пусть наша исходная функция, представляет собой сигнал, записанный в течение Т=1 сек. Тогда период первой гармоники будет равен длительности нашего сигнала Т1=Т=1 сек и частота гармоники равна 1 Гц. Период второй гармоники будет равен длительности сигнала, деленной на 2 (Т2=Т/2=0,5 сек) и частота равна 2 Гц. Для третьей гармоники Т3=Т/3 сек и частота равна 3 Гц. И так далее.
Шаг между гармониками в этом случае равен 1 Гц.
Таким образом сигнал длительностью 1 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 1 Гц.
Чтобы увеличить разрешение в 2 раза до 0,5 Гц - надо увеличить длительность измерения в 2 раза - до 2 сек. Сигнал длительностью 10 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 0,1 Гц. Других способов увеличить разрешение по частоте нет.
Существует способ искусственного увеличения длительности сигнала путем добавления нулей к массиву отсчетов. Но реальную разрешающую способность по частоте он не увеличивает.
3. Дискретные сигналы и дискретное преобразование Фурье
С развитием цифровой техники изменились и способы хранения данных измерений (сигналов). Если раньше сигнал мог записываться на магнитофон и храниться на ленте в аналоговом виде, то сейчас сигналы оцифровываются и хранятся в файлах в памяти компьютера в виде набора чисел (отсчетов).Обычная схема измерения и оцифровки сигнала выглядит следующим образом.
рис.9 Схема измерительного канала
Сигнал с измерительного преобразователя поступает на АЦП в течение периода времени Т. Полученные за время Т отсчеты сигнала (выборка) передаются в компьютер и сохраняются в памяти.
рис.10 Оцифрованный сигнал - N отсчетов полученных за время Т
Какие требования выдвигаются к параметрам оцифровки сигнала? Устройство, преобразующее входной аналоговый сигнал в дискретный код (цифровой сигнал) называется аналого-цифровой преобразователь (АЦП, англ. Analog-to-digital converter, ADC) (Wiki).
Одним из основных параметров АЦП является максимальная частота дискретизации (или частота семплирования, англ. sample rate) - частота взятия отсчетов непрерывного во времени сигнала при его дискретизации. Измеряется в герцах. ((Wiki))
Согласно теореме Котельникова, если непрерывный сигнал имеет спектр, ограниченный частотой Fмакс, то он может быть полностью и однозначно восстановлен по его дискретным отсчетам, взятым через интервалы времени 
, т.е. с частотой Fd ≥ 2*Fмакс, где Fd - частота дискретизации; Fмакс - максимальная частота спектра сигнала. Другими слова частота оцифровки сигнала (частота дискретизации АЦП) должна как минимум в 2 раза превышать максимальную частоту сигнала, который мы хотим измерить.
А что будет, если мы будем брать отсчеты с меньшей частотой, чем требуется по теореме Котельникова?
В этом случае возникает эффект «алиасинга» (он же стробоскопический эффект, муаровый эффект), при котором сигнал высокой частоты после оцифровки превращается в сигнал низкой частоты, которого на самом деле не существует. На рис. 11 красная синусоида высокой частоты - это реальный сигнал. Синяя синусоида более низкой частоты - фиктивный сигнал, возникающий вследствие того, за время взятия отсчета успевает пройти больше, чем пол-периода высокочастотного сигнала.
Рис. 11. Появление ложного сигнала низкой частоты при недостаточно высокой частоте дискретизации
Чтобы избежать эффекта алиасинга перед АЦП ставят специальный антиалиасинговый фильтр - ФНЧ (фильтр нижних частот), который пропускает частоты ниже половины частоты дискретизации АЦП, а более высокие частоты зарезает.
Для того, чтобы вычислить спектр сигнала по его дискретным отсчетам используется дискретное преобразование Фурье (ДПФ). Отметим еще раз, что спектр дискретного сигнала «по определению» ограничен частотой Fмакс, меньшей половине частоты дискретизации Fd. Поэтому спектр дискретного сигнала может быть представлен суммой конечного числа гармоник, в отличие от бесконечной суммы для ряда Фурье непрерывного сигнала, спектр которого может быть неограничен. Согласно теореме Котельникова максимальная частота гармоники должна быть такой, чтобы на нее приходилось как минимум два отсчета, поэтому число гармоник равно половине числа отсчетов дискретного сигнала. То есть если в выборке имется N отсчетов, то число гармоник в спектре будет равно N/2.
Рассмотрим теперь дискретное преобразование Фурье (ДПФ).
Сравнивая с рядом Фурье
Видим, что они совпадают, за исключением того, что время в ДПФ имеет дискретный характер и число гармоник ограничено величиной N/2 - половиной числа отсчетов.
Формулы ДПФ записываются в безразмерных целых переменных k, s, где k – номера отсчетов сигнала, s – номера спектральных составляющих.
Величина s показывает количество полных колебаний гармоники на периоде Т (длительности измерения сигнала). Дискретное преобразование Фурье используется для нахождения амплитуд и фаз гармоник численным методом, т.е. «на компьютере»
Возвращаясь к результатам, полученным в начале. Как уже было сказано выше, при разложении в ряд Фурье непериодической функции (нашего сигнала), полученный ряд Фурье фактически соответствует периодической функции с периодом Т. (рис.12).
рис.12 Периодическая функция f(x) с периодом Т0, с периодом измерения Т>T0
Как видно на рис.12 функция f(x) периодическая с периодом Т0. Однако из-за того, что длительность измерительной выборки Т не совпадает с периодом функции Т0, функция, получаемая как ряд Фурье, имеет разрыв в точке Т. В результате спектр данной функции будет содержать большое количество высокочастотных гармоник. Если бы длительность измерительной выборки Т совпадала с периодом функции Т0, то в полученном после преобразования Фурье спектре присутствовала бы только первая гармоника (синусоида с периодом равным длительности выборки), поскольку функция f(x) представляет собой синусоиду.
Другими словами, программа ДПФ «не знает», что наш сигнал представляет собой «кусок синусоиды», а пытается представить в виде ряда периодическую функцию, которая имеет разрыв из-за нестыковки отдельных кусков синусоиды.
В результате в спектре появляются гармоники, которые должны в сумме изобразить форму функции, включая этот разрыв.
Таким образом, чтобы получить «правильный» спектр сигнала, являющегося суммой нескольких синусоид с разными периодами, необходимо чтобы на периоде измерения сигнала укладывалось целое число периодов каждой синусоиды. На практике это условие можно выполнить при достаточно большой длительности измерения сигнала.
Рис.13 Пример функции и спектра сигнала кинематической погрешности редуктора
При меньшей длительности картина будет выглядеть «хуже»:
Рис.14 Пример функции и спектра сигнала вибрации ротора
На практике бывает сложно понять, где «реальные составляющие», а где «артефакты», вызванные некратностью периодов составляющих и длительности выборки сигнала или «скачками и разрывами» формы сигнала. Конечно слова «реальные составляющие» и «артефакты» не зря взяты в кавычки. Наличие на графике спектра множества гармоник не означает, что наш сигнал в реальности из них «состоит». Это все равно что считать, будто число 7 «состоит» из чисел 3 и 4. Число 7 можно представить в виде суммы чисел 3 и 4 - это правильно.
Так и наш сигнал… а вернее даже не «наш сигнал», а периодическую функцию, составленную путем повторения нашего сигнала (выборки) можно представить в виде суммы гармоник (синусоид) с определенными амплитудами и фазами. Но во многих важных для практики случаях (см. рисунки выше) действительно можно связать полученные в спектре гармоники и с реальными процессами, имеющими циклический характер и вносящими значительный вклад в форму сигнала.
Некоторые итоги
1. Реальный измеренный сигнал, длительностью T сек, оцифрованный АЦП, то есть представленный набором дискретных отсчетов (N штук), имеет дискретный непериодический спектр, представленный набором гармоник (N/2 штук).2. Сигнал представлен набором действительных значений и его спектр представлен набором действительных значений. Частоты гармоник положительны. То, что математикам бывает удобнее представить спектр в комплексной форме с использованием отрицательных частот не значит, что «так правильно» и «так всегда надо делать».
3. Сигнал, измеренный на отрезке времени Т определен только на отрезке времени Т. Что было до того, как мы начали измерять сигнал, и что будет после того - науке это неизвестно. И в нашем случае - неинтересно. ДПФ ограниченного во времени сигнала дает его «настоящий» спектр, в том смысле, что при определенных условиях позволяет вычислить амплитуду и частоту его составляющих.
Использованные материалы и другие полезные материалы.
Для цифровой обработки сигналов (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 раза количества информации происходит из-за того, что данные буфера имеют избыточность в виде дублирования значений.
Этот ряд может быть также записан в виде:
(2),
где , k-я комплексная амплитуда.
Связь между коэффициентами (1) и (3) выражается следующими формулами:
Отметим, что все эти три представления ряда Фурье совершенно равнозначны. Иногда при работе с рядами Фурье бывает удобнее использовать вместо синусов и косинусов экспоненты мнимого аргумента, то есть использовать преобразование Фурье в комплексной форме. Но нам удобно использовать формулу (1), где ряд Фурье представлен в виде суммы косинусоид с соответствующими амплитудами и фазами. В любом случае неправильно говорить, что результатом преобразования Фурье действительного сигнала будут комплексные амплитуды гармоник. Как правильно говорится в Вики «Преобразование Фурье (?) - операция, сопоставляющая одной функции вещественной переменной другую функцию, также вещественной переменной.»
Итого:
Математической основой спектрального анализа сигналов является преобразование Фурье.
Преобразование Фурье позволяет представить непрерывную функцию f(x) (сигнал), определенную на отрезке {0, T} в виде суммы бесконечного числа (бесконечного ряда) тригонометрических функций (синусоид и\или косинусоид) с определёнными амплитудами и фазами, также рассматриваемых на отрезке {0, T}. Такой ряд называется рядом Фурье.
Отметим еще некоторые моменты, понимание которых требуется для правильного применения преобразования Фурье к анализу сигналов. Если рассмотреть ряд Фурье (сумму синусоид) на всей оси Х, то можно увидеть, что вне отрезка {0, T} функция представленная рядом Фурье будет будет периодически повторять нашу функцию.
Например, на графике рис.7 исходная функция определена на отрезке {-T\2, +T\2}, а ряд Фурье представляет периодическую функцию, определенную на всей оси х.
Это происходит потому, что синусоиды сами являются периодическими функциями, соответственно и их сумма будет периодической функцией.
рис.7 Представление непериодической исходной функции рядом Фурье
Таким образом:
Наша исходная функция - непрерывная, непериодическая, определена на некотором отрезке длиной T.
Спектр этой функции - дискретный, то есть представлен в виде бесконечного ряда гармонических составляющих - ряда Фурье.
По факту, рядом Фурье определяется некоторая периодическая функция, совпадающая с нашей на отрезке {0, T}, но для нас эта периодичность не существенна.
Периоды гармонических составляющих кратны величине отрезка {0, T}, на котором определена исходная функция f(x). Другими словами, периоды гармоник кратны длительности измерения сигнала. Например, период первой гармоники ряда Фурье равен интервалу Т, на котором определена функция f(x). Период второй гармоники ряда Фурье равен интервалу Т/2. И так далее (см. рис. 8).
рис.8 Периоды (частоты) гармонических составляющих ряда Фурье (здесь Т=2?)
Соответственно, частоты гармонических составляющих кратны величине 1/Т. То есть частоты гармонических составляющих Fk равны Fk= к\Т, где к пробегает значения от 0 до?, например к=0 F0=0; к=1 F1=1\T; к=2 F2=2\T; к=3 F3=3\T;… Fk= к\Т (при нулевой частоте - постоянная составляющая).
Пусть наша исходная функция, представляет собой сигнал, записанный в течение Т=1 сек. Тогда период первой гармоники будет равен длительности нашего сигнала Т1=Т=1 сек и частота гармоники равна 1 Гц. Период второй гармоники будет равен длительности сигнала, деленной на 2 (Т2=Т/2=0,5 сек) и частота равна 2 Гц. Для третьей гармоники Т3=Т/3 сек и частота равна 3 Гц. И так далее.
Шаг между гармониками в этом случае равен 1 Гц.
Таким образом сигнал длительностью 1 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 1 Гц.
Чтобы увеличить разрешение в 2 раза до 0,5 Гц - надо увеличить длительность измерения в 2 раза - до 2 сек. Сигнал длительностью 10 сек можно разложить на гармонические составляющие (получить спектр) с разрешением по частоте 0,1 Гц. Других способов увеличить разрешение по частоте нет.
Существует способ искусственного увеличения длительности сигнала путем добавления нулей к массиву отсчетов. Но реальную разрешающую способность по частоте он не увеличивает.
3. Дискретные сигналы и дискретное преобразование Фурье
С развитием цифровой техники изменились и способы хранения данных измерений (сигналов). Если раньше сигнал мог записываться на магнитофон и храниться на ленте в аналоговом виде, то сейчас сигналы оцифровываются и хранятся в файлах в памяти компьютера в виде набора чисел (отсчетов).Обычная схема измерения и оцифровки сигнала выглядит следующим образом.
рис.9 Схема измерительного канала
Сигнал с измерительного преобразователя поступает на АЦП в течение периода времени Т. Полученные за время Т отсчеты сигнала (выборка) передаются в компьютер и сохраняются в памяти.
рис.10 Оцифрованный сигнал - N отсчетов полученных за время Т
Какие требования выдвигаются к параметрам оцифровки сигнала? Устройство, преобразующее входной аналоговый сигнал в дискретный код (цифровой сигнал) называется аналого-цифровой преобразователь (АЦП, англ. Analog-to-digital converter, ADC) (Wiki).
Одним из основных параметров АЦП является максимальная частота дискретизации (или частота семплирования, англ. sample rate) - частота взятия отсчетов непрерывного во времени сигнала при его дискретизации. Измеряется в герцах. ((Wiki))
Согласно теореме Котельникова, если непрерывный сигнал имеет спектр, ограниченный частотой Fмакс, то он может быть полностью и однозначно восстановлен по его дискретным отсчетам, взятым через интервалы времени , т.е. с частотой Fd ? 2*Fмакс, где Fd - частота дискретизации; Fмакс - максимальная частота спектра сигнала. Другими слова частота оцифровки сигнала (частота дискретизации АЦП) должна как минимум в 2 раза превышать максимальную частоту сигнала, который мы хотим измерить.
А что будет, если мы будем брать отсчеты с меньшей частотой, чем требуется по теореме Котельникова?
В этом случае возникает эффект «алиасинга» (он же стробоскопический эффект, муаровый эффект), при котором сигнал высокой частоты после оцифровки превращается в сигнал низкой частоты, которого на самом деле не существует. На рис. 5 красная синусоида высокой частоты - это реальный сигнал. Синяя синусоида более низкой частоты - фиктивный сигнал, возникающий вследствие того, за время взятия отсчета успевает пройти больше, чем пол-периода высокочастотного сигнала.
Рис. 11. Появление ложного сигнала низкой частоты при недостаточно высокой частоте дискретизации
Чтобы избежать эффекта алиасинга перед АЦП ставят специальный антиалиасинговый фильтр - ФНЧ (фильтр нижних частот), который пропускает частоты ниже половины частоты дискретизации АЦП, а более высокие частоты зарезает.
Для того, чтобы вычислить спектр сигнала по его дискретным отсчетам используется дискретное преобразование Фурье (ДПФ). Отметим еще раз, что спектр дискретного сигнала «по определению» ограничен частотой Fмакс, меньшей половине частоты дискретизации Fd. Поэтому спектр дискретного сигнала может быть представлен суммой конечного числа гармоник, в отличие от бесконечной суммы для ряда Фурье непрерывного сигнала, спектр которого может быть неограничен. Согласно теореме Котельникова максимальная частота гармоники должна быть такой, чтобы на нее приходилось как минимум два отсчета, поэтому число гармоник равно половине числа отсчетов дискретного сигнала. То есть если в выборке имется N отсчетов, то число гармоник в спектре будет равно N/2.
Рассмотрим теперь дискретное преобразование Фурье (ДПФ).
Сравнивая с рядом Фурье
Видим, что они совпадают, за исключением того, что время в ДПФ имеет дискретный характер и число гармоник ограничено величиной N/2 - половиной числа отсчетов.
Формулы ДПФ записываются в безразмерных целых переменных k, s, где k – номера отсчетов сигнала, s – номера спектральных составляющих.
Величина s показывает количество полных колебаний гармоники на периоде Т (длительности измерения сигнала). Дискретное преобразование Фурье используется для нахождения амплитуд и фаз гармоник численным методом, т.е. «на компьютере»
Возвращаясь к результатам, полученным в начале. Как уже было сказано выше, при разложении в ряд Фурье непериодической функции (нашего сигнала), полученный ряд Фурье фактически соответствует периодической функции с периодом Т. (рис.12).
рис.12 Периодическая функция f(x) с периодом Т0, с периодом измерения Т>T0
Как видно на рис.12 функция f(x) периодическая с периодом Т0. Однако из-за того, что длительность измерительной выборки Т не совпадает с периодом функции Т0, функция, получаемая как ряд Фурье, имеет разрыв в точке Т. В результате спектр данной функции будет содержать большое количество высокочастотных гармоник. Если бы длительность измерительной выборки Т совпадала с периодом функции Т0, то в полученном после преобразования Фурье спектре присутствовала бы только первая гармоника (синусоида с периодом равным длительности выборки), поскольку функция f(x) представляет собой синусоиду.
Другими словами, программа ДПФ «не знает», что наш сигнал представляет собой «кусок синусоиды», а пытается представить в виде ряда периодическую функцию, которая имеет разрыв из-за нестыковки отдельных кусков синусоиды.
В результате в спектре появляются гармоники, которые должны в сумме изобразить форму функции, включая этот разрыв.
Таким образом, чтобы получить «правильный» спектр сигнала, являющегося суммой нескольких синусоид с разными периодами, необходимо чтобы на периоде измерения сигнала укладывалось целое число периодов каждой синусоиды. На практике это условие можно выполнить при достаточно большой длительности измерения сигнала.
Рис.13 Пример функции и спектра сигнала кинематической погрешности редуктора
При меньшей длительности картина будет выглядеть «хуже»:
Рис.14 Пример функции и спектра сигнала вибрации ротора
На практике бывает сложно понять, где «реальные составляющие», а где «артефакты», вызванные некратностью периодов составляющих и длительности выборки сигнала или «скачками и разрывами» формы сигнала. Конечно слова «реальные составляющие» и «артефакты» не зря взяты в кавычки. Наличие на графике спектра множества гармоник не означает, что наш сигнал в реальности из них «состоит». Это все равно что считать, будто число 7 «состоит» из чисел 3 и 4. Число 7 можно представить в виде суммы чисел 3 и 4 - это правильно.
Так и наш сигнал… а вернее даже не «наш сигнал», а периодическую функцию, составленную путем повторения нашего сигнала (выборки) можно представить в виде суммы гармоник (синусоид) с определенными амплитудами и фазами. Но во многих важных для практики случаях (см. рисунки выше) действительно можно связать полученные в спектре гармоники и с реальными процессами, имеющими циклический характер и вносящими значительный вклад в форму сигнала.
Некоторые итоги
1. Реальный измеренный сигнал, длительностью T сек, оцифрованный АЦП, то есть представленный набором дискретных отсчетов (N штук), имеет дискретный непериодический спектр, представленный набором гармоник (N/2 штук).2. Сигнал представлен набором действительных значений и его спектр представлен набором действительных значений. Частоты гармоник положительны. То, что математикам бывает удобнее представить спектр в комплексной форме с использованием отрицательных частот не значит, что «так правильно» и «так всегда надо делать».
3. Сигнал, измеренный на отрезке времени Т определен только на отрезке времени Т. Что было до того, как мы начали измерять сигнал, и что будет после того - науке это неизвестно. И в нашем случае - неинтересно. ДПФ ограниченного во времени сигнала дает его «настоящий» спектр, в том смысле, что при определенных условиях позволяет вычислить амплитуду и частоту его составляющих.
Использованные материалы и другие полезные материалы.
