Щоб перевірити, чи правильно працює програма, сформуємо масив відліків як суму двох синусоїд 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://ua.wikipedia.org/ wiki/Ряд_Фур'є .)
Цей ряд може бути записаний у вигляді:
(2),
де, k-я комплексна амплітуда.
Зв'язок між коефіцієнтами (1) та (3) виражається такими формулами:
Зазначимо, всі ці три уявлення низки Фур'є цілком рівнозначні. Іноді під час роботи з рядами Фур'є буває зручніше використовувати замість синусів і косінусів експоненти уявного аргументу, тобто використовувати перетворення Фур'є у комплексній формі. Але нам зручно використовувати формулу (1), де ряд Фур'є представлений у вигляді суми косінусоїд з відповідними амплітудами та фазами. У будь-якому разі неправильно говорити, що результатом перетворення Фур'є дійсного сигналу будуть комплексні амплітуди гармонік. Як правильно говориться в Вікі «Перетворення Фур'є (ℱ) - операція, яка зіставляє однієї функції речової змінної іншу функцію, також речової змінної.»
Разом:
Математичною основою спектрального аналізу сигналів є перетворення Фур'є.
Перетворення Фур'є дозволяє представити безперервну функцію f(x) (сигнал), визначену на відрізку (0, T) у вигляді суми нескінченного числа (нескінченного ряду) тригонометричних функцій (синусоїд та косинусоід) з певними амплітудами і фазами, також розглядаються на відрізку (0, T). Такий ряд називається поряд Фур'є.
Зазначимо ще деякі моменти, розуміння яких потрібно правильного застосування перетворення Фур'є до аналізу сигналів. Якщо розглянути ряд Фур'є (суму синусоїд) на всій осі Х, то можна побачити, що поза відрізком (0, T) функція представлена поруч Фур'є буде періодично повторювати нашу функцію.
Наприклад, на графіці рис.7 вихідна функція визначена на відрізку (-T2, + T2), а ряд Фур'є представляє періодичну функцію, визначену на всій осі х.
Це тому, що синусоїди самі є періодичними функціями, відповідно і їх сума буде періодичною функцією.
рис.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. Сигнал, виміряний на відрізку часу Т визначено тільки на відрізку часу Т. Що було до того, як ми почали вимірювати сигнал, і що буде після того – науці це невідомо. І в нашому випадку – нецікаво. ДПФ обмеженого у часі сигналу дає його «справжній» спектр, у тому сенсі, що за певних умов дозволяє обчислити амплітуду та частоту його складових.
Використані матеріали та інші корисні матеріали.
Теорія
Спочатку трохи теорії. Як відомо все в подібних аналізаторах використовується швидке перетворення Фур'є і часто говориться, що ДПФ в подібних конструкціях використовувати не можна, тільки БПФ та й то на асемблері. Я ж використав натомість дискретне перетворення Фур'є (ДПФ) та перетворення за Уолшем. І в цій статті доведу, що можна використовувати навіть не тільки БПФ, а ДПФ написаний на С. Але спочатку по порядку як з ДПФ отримати просту функцію ДПФ і Волшу. ДПФ класично виглядає так:
Так як у мк мало ресурсів, то замінюють 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.
Аналізатор спектру працює на частоті дискретизації 192 кГц від кварцу на 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 Гц 255 відліків, 400 Гц 240 відліків, 300 Гц 512 відліків, 150 Гц 512 відліків, 75 Гц 512 відліків.
Таким вибором числа відліків вдалося смугу рівномірної по всьому діапазону частот.
Ще один підводний камінь вийшов на частоті 96 кГц. Так як уявної частини немає (це легко перевірити підставивши у формулу вище при 512 відліках 256 номер спектру і синус буде завжди дорівнює 0), то реальна частина може сильно змінюватися за рахунок того, що значення косинуса буде обчислюватися через раз в протифазі до основного сигналу . Тобто обчислюватиметься раз. Для того, щоб цього не було, необхідно обчислити хоча б 2 значення реальної частини зсунутої на 90 градусів і вибрати максимальний з двох значень.
Алгоритм програми накопичення 512 відліків у проміжку переведення мк у режим сну та пробудження коли черговий відлік готовий. Крім того 150 Гц відбувається розгортка світлодіодів - це раз на 128 від частоти дискретизації в 19200. Тобто до того як ацп зніме всі відліки він встигне повністю провести одну повну розгортку. Як тільки всі відліки готові в основному циклі програми відбувається обчислення всіх амплітуд спектру. У цей час розгортка продовжується, але мк не впадає в сон, а вважає амплітуди. Щойно амплітуди пораховані мк перетворюється на сон і програма повторюється заново. Амплітуди розраховуються виходячи з 20 дБ діапазону, тобто прологарифмовані.
Виходячи з часу отримання всіх відліків і час на розрахунок всіх амплітуд частота оновлення знаходиться в районі 10-15 Гц.
Для цифрової обробки сигналів (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 для підключення до лінійного виходу пристрою і зібраний найпростіший вхідний ланцюжок. Як було зазначено, на платі вже встановлено динамік, підключений до аналогового виходу мікроконтролера і є елементи управління (кнопка і потенціометр).
Ось малюнок вхідного ланцюга:
Налагодження софту відбувалося поетапно:
- Налагодження всієї необхідної периферії: АЦП, ЦАП, таймери, світлодіодна індикація.
- Тест з оцифруванням сигналу: оцифровую дані з необхідною швидкістю і кладу в буфер, потім витягую дані з буфера і відтворюю сигнал. Тобто. простий зсув сигналу у часі, без будь-яких перетворень. На цьому етапі тестується механізм роботи з двома буферами, необхідний для подальшої роботи.
- До попередньої версії додається пряме та зворотне перетворення Фур'є. Цей тест остаточно перевіряє коректність роботи коду БПФ, а також перевірка, що код вкладається в доступну продуктивність.
- Після цього основний кістяк програми зроблено, можна перейти до практичних застосувань.
Проблема виникла після додавання до коду БПФ: у сигналі з'явилися сторонні шуми та свисти. Взагалі, ця поведінка видалася мені досить дивною, т.к. без перетворення сигнал проходячи цифровий тракт залишався досить чистим. Перша причина цього: після аналогового ланцюга, амплітуда сигналу на АЦП була не повною (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 нотами в октаві не набули поширення. Їх складно налаштовувати, на них складно грати, а відсоток людей, здатних відчути різницю із звичайними інструментами, зникаючий.
Всі сигнали, незалежно від того, ви їх придумали або спостерігали у Всесвіті, насправді просто сума простих синусоїд різних частот.
Я зробив невеликий аудіо аналізатор спектру (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 | Резистор | ||||||