Десь місяць тому (або трохи більше) зачепила мене тема визначення наявності частоти в зашумленому сигналі, що вирішується зокрема перетворенням Фур'є, і всякими кореляторними та іншими математичними прибамбасами... Звичайно, мені як людині далекої від вищої математики, захотілося простих шляхів . Методом гуглення було з'ясовано, що схоже завдання вирішують при детектруванні ДТМФ, і використовують алгоритм Герцеля (перетворення фур'є для однієї частоти). Вся красота цього алгоритму в тому, що не треба заганяти всі семпли в буфер (оперативка в мікроконтролері дефіцит), а можна виконувати розрахунок прямо в процесі отримання даних з АЦП, використовуючи в процесі всього три динамічних змінних для кожної частоти що шукається, і вже в кінці здійснити розрахунок для отримання ваги частоти в сигналі.
Після довгих пошуків був знайдений ексель файлик з прикладом алгоритму Герцеля, але він не працював на моєму комп'ютері (хоча і робив у гугл таблицях), до всього там використовувалася математика з числами з плаваючою комою, і деякі незрозумілі функції самого екселю, а я шукав алгоритм адаптований для мікроконтролерів. У результаті я взяв цей ексель за основу, оскільки в ньому була генерація масиву семплів з двох синусоїд. Довго перебираючі різні приклади коду алгоритму на гітхабі я все-таки народив робочу версію в екселі(буквально сьогодні, ділюся радістю). Причому версію яка не використовує навіть важкого поділу (тільки додавання, віднімання та множення), все брав з коду оптимізованого для мікроконтролера (в мікроконтроллері замість ділення використовується ссув на певну кількість біт, в екселі то замінив діленням на число кратне ступеню двійки). Хоча трохи брешу, фінальна операція це витяг кореня, але його робити не обов'язково (просто результатом буде амплітуда в квадраті).
Наприклад, я взяв дві частоти, 50Гц і 52Гц. 52Гц зробив максимальною амплітудою, а 50Гц амплітудою в МІЛЬЙОН РАЗ МЕНШЕ. На частоті 50Гц я отримував амплітуду 33 одиниці, на 52Гц я отримував амплітуду 32+мільйона одиниць. А на частоті 51Гц (яка між цими двома частотами) я отримував одну одиницю. Смуга детектування рівно 1Гц, при частоті семплів 1024Гц, та кількості обрахованих точок 1024 штуки (1 секундний інтервал). Тоб то, якщо нам потрібна роздільна здатність в 1Гц (детектувати зсув частоти на 1Гц, або дві різні частоти з сувом 1Гц), то потрібен інтервал спостереження в 1 секунду.
Потім я загрубив вихідну сиквенцію до рамок 12 бітного АЦП (+-2048), частоту 52Гц поставив макстимальної амплітуди, а частоту 50Гц зробив у 10000 (десять тисяч) разів слабше. І незважаючи на те, що вага молодшого розряду 12 бітного семплювання більша ніж амплітуда сигналу 50Гц, на цій частоті алгоритм показував тверді 5 одиниць, у той час як на пів герца в бік вже були звичайні від 0 до 1 одиниці тиші. Для мене особисто це магія якась із розряду незбагненного, тому з радістью ділюся. Тепер можна й мікроконтролер розчехляти.
ексель для погратися додаю(подробиці потім). На жовтому фоні введення частот, ліворуч задається жовта частота детектування, а праворуч дві жовті частоти генерації. На зеленому тлі задаються амплітуди частот що генеруються. На помаранчевому фоні результат амплітуди частоти детектування. Кількість семплів і частоту семпліровніе також можна змінювати.
Пару тижнів тому постала переді мною задача на детектування наявності двох частот в сигналі який міг бути періодично зашумлений ще двома частотами. Я спочатку зрозумів її не правильно і вирішив в лоб через швидке перетворення Фурьє з детектуванням всіх згаданих в завданні частот. Коли перечитав завдання вирішив вже нічого не переробляти, хоть і надлишково але задача була вирішена вірно. Мікроконтролера для даної задачі цілком вистачило, а додатково на нього чіпляти нічого не передбачалося. В процесі пошуку інформації для вирішення задачі натикався на алгоритм Герцеля, але пояснення було складною академічною мовою тому вирішив сильно не заглиблюватися і використав готову функцію на С знайдену в інтернеті для швидкого перетворення Фурьє. Цікаво було б глянути на реалізацію алгоритму Герцеля для мікроконтролера, можливо буде зрозуміліше ніж в тому викладенні на яке я натрапляв )
_________________ Знание приумножает скорбь и увеличивает размер исходного кода.
Тут 9615.0 це частота семплування (опитування) АЦП. Я в екселі рахую його по цій формулі, можна формулу глянути там (Koef). Коефіцієнти можна порахувати при ініціалізації, або не напрягати контроллер а порахувати руками та занести в контроллер як константи. Множення на 1 << 14 (16384) то щоб отримати число з фиксованною точкою де 14 біт після точки. Код алгоритму під спойлером.
Код:
//-------Goertzel function---------------------------------------// long int goertzel (int sample[], long int coeff, int N) //---------------------------------------------------------------// { //initialize variables to be used in the function int Q, Q_prev, Q_prev2, i; long prod1, prod2, prod3, power;
Q_prev = 0; //set delay element1 Q_prev as zero Q_prev2 = 0; //set delay element2 Q_prev2 as zero power = 0; //set power as zero
for (i = 0; i < N; i++) // loop N times and calculate Q, Q_prev, Q_prev2 at each iteration (цикл) { Q = (sample[i]) + ((coeff * Q_prev) >> 14) - (Q_prev2); // >>14 used as the coeff was used in Q15 format Q_prev2 = Q_prev; // shuffle delay elements Q_prev = Q; }
//calculate the three products used to calculate power (тут рахується вага частоти) prod1 = ((long) Q_prev * Q_prev); prod2 = ((long) Q_prev2 * Q_prev2); prod3 = ((long) Q_prev * coeff) >> 14; prod3 = (prod3 * Q_prev2);
power = ((prod1 + prod2 - prod3)) >> 8; //calculate power using the three products and scale the result down
">> 14" то бітовий зсув після множення щоб знову отримату фіксовану точку (яка при множенні зсувається), як що бачите то в коді то це характерний признак що алгоритм оптимізований для мікроконтроллера під числа з фіксованою точкою. В циклі вираховується Q з кожного семплу, і зберігається ще дві її копії, Q_prev - попередня та Q_prev2 - позапопередня. На початку в них записують нулі. Коли всі семпли прораховані то вираховується амплітуда частоти в сигналі трохи складнішою формулою.
Весь кайф в тому, що значення АЦП можна брати в циклі не з массиву (буферу в пам'яті), а прямо читати з АЦП привї'язавши виконання циклу до преривань АЦП. То б то не треба буфера для семплів. В такому циклі можна паралельно рахувати ще кілька частот (додасться ще парочка змінних на кожну частоту). Виходить можна те робити на дуже бідному до пам'яті мікроконтроллері.
Для відносно плавного безперервного детектування частоти роблять тандемний цикл (коли один дорахував до половини то другий вже закінчив і починає знов), такі цикли перекриваються або повністью, або частково. Ще є алгоритм Герцеля ковзного типу в якому закладено вікно чутивості, але я знаходив про нього англомовні публікації які не мали самого алгоритму який я міг бі зрозуміти (максимум якісь багатоповерхові формули з незрозумілими знаками).
_________________ 2 моноколеса, 2 электровела, 2 подводных буксира, надувной электро-каяк.
Зараз цей форум переглядають: немає зареєстрованих користувачів і гості: 56
Ви не можете створювати нові теми Ви не можете відповідати на повідомлення Ви не можете редагувати свої повідомлення Ви не можете видаляти свої повідомлення Ви не можете додавати додаткові файли
Вхід
Реєстрація
FAQ
Пошук
Користувачі
Топлист подяк
Поточний час: 27 квіт 2024, 21:05
286 
1009
