Jump to content

    

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

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

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now