Jump to content

    
haker_fox

Коллеги, что порекомендуете для вычисления RMS по выборкам?

Recommended Posts

17 minutes ago, Arlleex said:

Плохо? Хорошо?

Я тоже так думал, по-сути это метод определения частоты методом "перехода через ноль". Но дело в том, что синусоида в сети нечистая, а вполне себе может содержать различные всплески, и как раз в точке экстремума. И тогда весь алгоритм - кисе под хвост))) Нужны именно методы, основанные на математическом анализе, которые "видят" всю картину целиком.

Share this post


Link to post
Share on other sites

Ну почему...
Вот, допустим. Черным - "идеальный" синус, красным - реальный.

image.gif

 

Вот это мы накопили АЦП.

Анализируем слева направо: точка А принимается в расчет как возможный экстремум. Но двигаемся дальше - точка B принимается теперь как возможный экстремум, потому что его амплитуда больше, чем у A. Точка C тоже показывает экстремум, но она меньше, чем текущая (B), значит, это не максимальное значение функции. Дальше встречаем переход сигнала через 0. С этого момента искать экстремум нет необходимости (в отрицательной полуволне). В этой точке также запоминаем количество выборок между началом и первым переходом через 0, ведь это - возможный полупериод сигнала.

Ищем снова переход через 0: следующий он будет через полупериод сигнала. По аналогии с точками A, B и C, обнаруживаем экстремум в точке E среди экстремумов D, E, F. Запоминаем его.

Теперь у нас есть: расстояние между точками B и E, а также значение, которое ранее мы приняли за возможный полупериод. Если этот возможный полупериод, умноженный на 2, совпадает с расстоянием между точками B и E с некоторой погрешностью, то частота определена верно. Если не совпадает - значит один из экстремумов или переход через 0 являлся всплеском и алгоритм начинается сначала. В установившемся режиме всплесков не будет, либо они будут периодичны, и, таким образом, входить в состав сигнала для расчета RMS. Ну либо надо покумекать, как непереодичные всплески фильтровать:bb:

Edited by Arlleex

Share this post


Link to post
Share on other sites

Вот и я о том же. В силу моего отдалённого знакомства с гильбертами, я бы попытался решать эту задачу по-простому, "по-деревенски". Сначала набрал бы количество сэмплов заведомо большее, чем нужно для самого длительного периода. Затем бы нашёл экстремумы и по ним прикинул период (и нужное количество сэмплов). Посчитал бы среднее. Следующим этапом подвигал бы количество сэмплов туда-сюда, наблюдая за близостью среднего к нулю, и определили бы точнее это количество, соответствующее фактическому значению периода. И только после этого считал бы RMS. По нескольким следующим периодам. В предположении, что частота меняется не столь быстро.

Понятно, что сеть грязновата, но Вы ведь сами сказали, что сигнал, в целом, достаточно близок к синусу. Поэтому фильтрация помогла бы в значительной степени адекватному определению периода.

Во-первых, перед АЦП всё равно нужен НЧ фильтр, хотя бы для antialiasing-а. Во-вторых, медианный фильтр убрал бы случайные всплески. И, в-третьих, ещё немного подфильтровать  совсем не вредно. Если не стоит задача найти энергетическое распределение по гармоникам, то этот путь вполне легитимен: общая энергия при фильтрации не потеряется.

Сильно ли я упрощаю?

6 hours ago, iiv said:

так я о том же, в алгоритме, что я привел, всего 13 существенных строчек, и ни каких Фурье и Гильбертов.

Кстати, примерно 5 лет назад я что-то похожее Вам рассказывал, но, к сожалению, не смог найти ссылку, поэтому запрограммировал заново :)

Да, спасибо. В тот раз мне не удалось разобраться в алгоритме, посмотрю внимательнее сейчас. Было бы здорово, если бы Вы смогли сделать некоторые комментарии к Вашему листингу. Не сочтите за придирку.

Share this post


Link to post
Share on other sites
11 hours ago, haker_fox said:

А не подскажете, как этот алгоритм называется? Или это просто здравый смысл?))))

Линейное предсказание по осциллирующей функции.

 

Вывод примерно такой:

 

пусть у Вас дана функция f(x), задискретизированная на f_i, i=1,...N и мы предполагаем, что она описывается одной функцией A*cos(Alpha*x)+B*sin(Alpha*x), (в общем случае это будет сумма таких функций с разными Alpha).

Если рассмотреть подпространство f(x), f(x+d), f(x+2*d), то легко заметить, что оно не полное и имеет ранг 2, в комплексном случае ранг 1, а в случае нескольких Alpha, надо взять больше таких сдвигов и все равно мы наткнемся на случай, когда подпространство будет меньшей размерности, чем число таких функций. Очевидно, что разумно взять d таким, что он соответствует скорости семплирования, тогда достаточно построить подпространство в виде f_i, f_{i+1}, f_{i+2}

 

Итак, это подпространство имеет нулевое (в случае отсутствия шума) сингулярное число, и правый его вектор соответствует коэффициентам полинома, корни которого суть Alpha.

То есть алгоритм строит довольно хорошее приближение к этому минимальному сингулярному числу и соответствующему вектор, в моем случае правый сингулярный вектор имеет вид (1, a2, a3). Далее - квадратное уравнение и его корень и, искомая Alpha.

 

Понятно, что если у вас куча гармоник, то надо брать больше f_i, f_{i+1}, ... f_{i+k} и можно все их найти.

В алгоритме, что я привел, предполагается, что минимальное сингулярное число равно нулю, но это может плохо работать в случае если есть еще гармоники близкие по уровню к основной частоте. Есть аналогичный алгоритм, с честным сингулярным числом, он требует немного больше вычислений, но ИМХО, с Вашим семплингом он даже в атмегу влезет, могу тоже привести, там на пару десятков строк больше из-за честного вычисления минимального сингулярного числа с его правым вектором.

 

Пользую с 1997 (на каком-то аспирантском курсе было), с тех пор использовал очень много раз, всегда работало. Ссылку дать не могу, ибо не знаю на кого точно сослаться - в таком виде это уже классика, скорей всего первый это сформулировал Колмогоров в 1941, но оригинал его не читал (не нашел), да и может кто до него это тоже делал, посему врать не буду. Аналог в многомерном случае сам недавно тиснул, ссылку на многомерную задачу могу дать. Если что-то непонятно объяснил - спрашивайте.

Share this post


Link to post
Share on other sites
2 hours ago, iiv said:

Итак, это подпространство имеет нулевое (в случае отсутствия шума) сингулярное число...

Вы, похоже, математик.. ;)

 

Напомнило анекдот:

Quote

Воздушный шар сбился с курса, и воздухоплаватель срочно опустился с ним вниз. Увидев внизу человека, он спросил:
- Извините, где я нахожусь?
- Вы находитесь на воздушном шаре, в 15м над землей. Ваши координаты - 5°28'17" N и 100°40'19" E.
- Похоже, вы математик, - вздохнул воздухоплаватель.
- Да, я математик, - согласился прохожий. - Как вы догадались?
- Ваш ответ, по-видимому, точный и полный, но для меня совершенно бесполезный. Я по-прежнему не знаю, где я нахожусь, и что мне делать. Вы мне нисколько не помогли, только напрасно отняли время. 

 

Share this post


Link to post
Share on other sites
1 minute ago, blackfin said:

Вы, очевидно, математик.. ;)

да, верно, только стоит уточнить, что я - вычислительный математик, посему я предложил простое решение в виде готовой, оттестированной и рабочей программы на основе проверенной десятилетиями теории - бери и пользуйся, не так ли?

Share this post


Link to post
Share on other sites
5 minutes ago, iiv said:

на основе проверенной десятилетиями теории - бери и пользуйся, не так ли?

Не аргумент.. Теория рядов Фурье проверена уже столетиями - "бери и пользуйся, не так ли?" ;)

Share this post


Link to post
Share on other sites
Just now, blackfin said:

Не аргумент.. Теория рядов Фурье проверена уже столетиями - "бери и пользуйся, не так ли?" ;) 

Если так просто, взяли бы и накропали бы в 9 операторов Фурье с поиском максимума, а мы бы тут все вместе поржали бы сколько килобайт кода у Вас получилось бы и во сколько сотен раз он бы проигрывал бы мной предложенному алгоритму.

 

ТС решает задачу в реальном времени на очень дохлом процессоре, и вычислительная сложность 5 N Log_2 N против 3N + подбор окна под Фурье означает в сотни раз вычислительно более сложную задачу. В моей программе всего-то 9 операторов, Фурье в такое число операторов с поиском максимума и получением правильного отскалированного результата Вы не уложитесь, чуть-чуть не хватит :) Чисто Фурье - да, можно уложиться, только оно работать очень медленно будет, а со всем обвесом - не уложитесь. Не говоря уж о том, что для Фурье гарантированно надобно окно в 4N чисел (а лучше больше), а в мною предложенном алгоритме надо только несколько переменных.

 

А если вспомнить, что Фурье даст дискретный максимум и Вы будете приплясывать аппроксимируя дельта функцию, чтоб уточнить результат, а в мной предложенном алгоритме ответ получается изначально не дискретным, а точным с точностью до шумов, Фурье опять проигрывает.

 

Не желание понять и прочитать хотя бы Колмогорова или любую классику по линейному предсказанию не означает, что этих методов нет или они хуже :(

Share this post


Link to post
Share on other sites
1 minute ago, iiv said:

В моей программе всего-то 9 операторов..

Не нужно так кричать. Я вас услышал.. ;)

 

В моей программе всего-то 1 оператор - умножение на функцию exp(-a*t^2).. взятую из таблицы.. ;)

Share this post


Link to post
Share on other sites
1 minute ago, blackfin said:

В моей программе всего-то 1 оператор - умножение на функцию exp(-a*t^2).. взятую из таблицы.. ;) 

О как весело, Вы предлагаете ТС умножать в лоб на матрицу Фурье ввергнув его в N*N арифметическую сложность вместо линейной что у меня и N log N что в БПФ? Самому-то не стыдно такую глупость писать и советовать?

 

Или это упертая зависть, что Вы даже умножение на матрицу Фурье на раз запрограммировать не сможете не запутавшись в комплексных коэффициентах, а тут кто-то взял да и поделился алгоритмом, причем работающим.

Share this post


Link to post
Share on other sites
5 hours ago, iiv said:

О как весело, Вы предлагаете ТС умножать в лоб на матрицу Фурье

В моей программе Фурье не вычисляется. Входной сигнал умножается на exp(-a*t^2), полученное произведение возводится в квадрат, после чего суммируется по всем точкам на интервале измерения. Затем полученная сумма делится на вес оконной функции Гаусса и вычисляется квадратный корень.

 

Вуаля! RMS готов!

 

А теперь "тут все вместе можете поржать".. :)

 

Share this post


Link to post
Share on other sites

Коллеги, предлагаю сохранять вежливое отношение к друг другу:acute: ибо за вашими эмоциями и меряниями количеством операторов, трудно отследить истину!

Share this post


Link to post
Share on other sites

Интересно почитать посты математически подкованных людей:wink:

Жаль только, что я на матанализе в универе, хоть и учился более-менее хорошо, не мог "приткнуться" с практической точки зрения применения абсолютного большинства задач, которые там проходили.

Но, как мне кажется, нужно не забывать о реальном конечном быстродействии вычислительного устройства (МК), который, помимо подсчета RMS, уверен, у автора должен выполнять еще ряд задач.

Edited by Arlleex

Share this post


Link to post
Share on other sites
42 minutes ago, Arlleex said:

нужно не забывать о реальном конечном быстродействии вычислительного устройства (МК)

Те, кому не хватает быстродействия МК, обычно просто берут стандартный MATLAB'овский алгоритм Exponential Weighting Method и отливают его в железе.. ;)

 

Точность, правда, не гарантирована.. Какая получится, такая получится.. ;)

 

 

Share this post


Link to post
Share on other sites
21 hours ago, blackfin said:

Те, кому не хватает быстродействия МК

Трудно оценить запас быстродействия... но... стандартный набор: работа с сетью, microSD, флешками, графика 128 на 64... всё скромненько, как у всех)))) С наступающим Новым Годом!!!

Share this post


Link to post
Share on other sites

Join the conversation

You can post now and register later. If you have an account, sign in now to post with your account.

Guest
Reply to this topic...

×   Pasted as rich text.   Paste as plain text instead

  Only 75 emoji are allowed.

×   Your link has been automatically embedded.   Display as a link instead

×   Your previous content has been restored.   Clear editor

×   You cannot paste images directly. Upload or insert images from URL.