Да, примерно так. На окне наблюдения делал ресемплинг, чтобы фильтр Фурье работал с тем окном, которое надо (условно, 1/50=0.02 с для 50 Гц и 1/51=0.019608
для 51 Гц).
Ресемплинг можно делать кусочно-линейный (это самый быстрый), либо можно заморочиться с монотонными сплайнами.

Я брал реальные осциллограммы с КЗ и строил для них частоту от времени. В осциллограммах была какая-то апериодика у токов, какие-то небольшие махры в напряжениях... более-менее нормально определялась частота 50 Гц +- какие-то мелкие дроби (что-нибудь типа 49.95 Гц).
По напряжениям частота хорошо считалась (особенно если обработать данные по всем трем фазам). По токам хуже (апериодика мешала и всякие искажения синусоиды токов).

Вот
comtradematlabparserv091.zip

Запускать TEST_OSCILL_HANDLE.m

(также можно прочесть комтрейд файлом TEST_READ_COMTRADE и записать что-то в комтрейд файлом TEST_WRITE_COMTRADE)

Мне кажется идейно это похоже на ФАПЧ. И там и тут нужно так подстроить частоту фильтра Фурье и период интегрирования чтобы вектор остановился (равно дисперсия стала минимальной). В ФАПЧ это делается в контуре обратной связи по разностной частоте, а тут в лоб вычисляется подбором нужной частоты (или более быстрым методом). Вычислений, конечно, много (если правильно понял). Надо посчитать Фурье на каждой точке да еще и несколько раз, пока процесс не сойдется. Не все современные терминалы всё успевают посчитать даже один раз между соседними точками, поэтому переходят на расчеты раз в 2.5, 5, 10 мс.

А можно чуть подробнее по алгоритму? Ищем абсциссу минимума функции D(f)? Эта функция примерно похожа на параболу? Получается что при идеальных сигналах и идеальной настройке и дисперсия и производная дисперсии равны нулю, что приводит к соотношению 0/0 в выражении для обновления f.

Добавлено: 2021-03-31 21:55:14

А какие?

Мне кажется несоответствие немного другое.
Процесс вида:
x(n) = x(n-1) - f ( x(n-1) ) / f' ( x(n-1) )
ищет не минимум функции f(x), а ноль (корень). И это логично, сама функция стремится к нулю, а производная нет и на каждом шаге добавка всё уменьшается и уменьшается. А в районе минимума функции данный процесс может начать колебания рядом с минимумом.
Для поиска минимума функции f(x) ищут корень функции f'(x), значит процесс для поиска минимума функции дисперсии должен выглядеть как
x(n) = x(n-1) - f' ( x(n-1) ) / f'' ( x(n-1) )

Дисперсия сигнала D(f) есть функция от f
мы фильтруем массив значений сигнала x фильтром Ф(f) несколько раз,
получаем набор отфильтрованных отсчетов Xdft(f)=filter(x,f)
далее находим дисперсию D(x_all, f)=D(filter(x1,f)...filter(x2,f)....filter(xn,f)), x_all=[x1,x2....xn],  сами x1,x2... это тоже массивы, выделяемые из окна данных, со сдвигом.

Далее нам надо подобрать такую частоту f, чтобы D(x_all,f) была равна нулю, т.е. получается уравнение D(f)=0.
f получается корень. df на каждой итерации по Ньютону уменьшается, поэтому, если все правильно сделать, колебаний быть не должно (разве что какие-то мелкие, из-за округления).

Зачем там минимум и вторая производная? Тем более, что все равно численно производную придется искать.

А, понял.
Да, наверное так лучше.

Добавлено: 2021-04-01 12:49:02

Или действительно попробовать другие способы поиска минимума.

ОК,

Вы знаете условия, при которых проводят сертификацию ?,

в них и должно быть задано отношение сигнал/шум

Смысл в том, что у нас не америка и вместо 50 Гц не может оказаться 60 Гц или там 30 Гц или сигнал с какими-то очень мощными ВЧ-всплесками. Там везде синусоида примерно 50 Гц частотой, не очень загрязненная, для нее вполне реально сделать определение частоты за малый промежуток времени.
Если бы у нас был какой-то радиосигнал, где есть куча разных частот, тогда да, надо было бы брать большое окно данных.

Я потестил этот алгоритм, который выше. Там есть определенные затыки с Ньютоном, т.е. всякие кривые значения, которые иногда он выдает, связаны с кривизной алгоритма поиска минимума дисперсии (точка минимума пролетает, и далее начинаются кривые итерации "не туда").
Но для контроля сделал поиск минимума по методу золотого сечения. Работает хорошо (брал сигнал с рандомным шумом в несколько процентов). Но долго. Надо еще думать над производительностью и тестировать.

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

Я так понимаю, весь смысл всех этих математических "изощрений" заключается в том, что если нельзя просто определить частоту сигнала как число колебаний в единицу времени путём измерения периода (или полупериода), то вычисление частоты производится путём перебора или подстановки значений частоты, пока на и большим образом не будет удовлетворять точности. Действительно, если частота синусоиды 50 Гц и взяли это значение, то дисперсия будет ноль, более того если она скажем 49,97747812731 и взять такую же, то будет также...

Где-то ещё лет 25 назад таким способом вычислял частоты двух эквивалентных эдс при АР, путём точного измерения периода АР, всё было замечательно если он был постоянным, но как только брал натурные сигналы АР, всё разваливалось, поэтому и забросил (да и компьютера не было хорошего, считал на MatchCad 2.0 ещё под DOS на Мозовия с тактовой 2,38 МГц - один опыт считался 3...4 мин.).