Оконное взвешивание, алгоритм применения Опции

[Zelepuk]

Есть задача: нужно посчитать THD в промышленной сети 50Гц.
THD считается по 10 гармоникам (кратным оновной частоте). Основная частота меняется в пределах 45...55Гц.
Я хочу применить алгоритм:

-оцифровываю сигнал с частотой 4096Гц (это априорная величина задана жёстко),
-затем оконное взвешивание (использую flattopwin в MATLAB - окно с плоской вершиной),
-делаю БПФ на 512 точек (окно в 0,125с.)

Теперь самое интересное, если гармоник 10, то нужно просто найти 10 пиков в спектре и вокруг взять корень квадратный из суммы квадратов окружающих эти пики бинов. Но если частота основной гармоники плавает от 45 до 55 Гц, то получается например для 9-й гармоники бины на которых основная энергия меняются от 52 до 59 (при частоте основной гармоники 45-55Гц), а для 10-й Гармоники основная энергия распределяется на 58-70 бинах (в зависимости от изменения основной гармоники), 11-я гармоника уже бегает с 63 (при 45 гц основной гармоники) до 77 (при 55 гц основной гармоники) бина.

Иначе говоря трудно ожидать пики в одних и тех же местах спектра при изменении частоты основной гармоники, посему ищется способ однозначного определения местоположения пиков (амплитуд) гармоник в спектре. Речь идёт именно о поиске, так как на графиках точно видно что амплитуды всех гармоник имеют место быть с достаточно выскоой точностью.

[SPACUM]

Ошибка SPACUM

Был не прав. Прошу прощения. Амплитуды выбросов АЧХ окна флаттоп на синусах и косинусах равны.
При оцифровке 12кГц и числе точек 2000 на висящем куске провода на логарифмическом спектре раздельно видны 50 гармоник.
Нечетные стоят, четные меняются.
Идеи:
1. Продифференцировать сигнал до АЦП чтобы совсем мелких гармоник не было.
2. На логарифмическом спектре использовать кепстр для определения 1/ fundamental frequenсy

[Alex11]

Коли спрашивали, признаюсь: оцифровка на частоте 25600 Гц, основная измеряемая гармоника - 50 Гц. При вычислении RMS используются все до Fs/2, при расчете гармонического состава - 40 гармоник. Но в последнем случае очень тяжело смотреть погрешности, т.к. все существующие эталоны дают отвратительную точность. То, что приводилось (0.01%) - это точность измерения RMS и амплитуды основной гармоники. Это на существующем приборе. При моделировании точность рассчетов (смесь float и double) на уровне 1e-8.

[SPACUM]

Коли спрашивали, признаюсь: оцифровка на частоте 25600 Гц, основная измеряемая гармоника - 50 Гц. При вычислении RMS используются все до Fs/2, при расчете гармонического состава - 40 гармоник. Но в последнем случае очень тяжело смотреть погрешности, т.к. все существующие эталоны дают отвратительную точность. То, что приводилось (0.01%) - это точность измерения RMS и амплитуды основной гармоники. Это на существующем приборе. При моделировании точность рассчетов (смесь float и double) на уровне 1e-8.

Большое спасибо за ответ. Мне такую задачу приходится решать при калибровке приборов и именно из-за очень плохих эталонов. В виброметрии других не бывает. Попробую Ваш метод. Но только в конце сентября. Работы много и отпуск. В реальном приборе точность расчетов (int32) на уровне 1e-8 для числа точек 2^n (на уровне 1.5e-7 для числа точек 2^n*10^n).

[Signal]

Есть задача: нужно посчитать THD в промышленной сети 50Гц.
THD считается по 10 гармоникам (кратным оновной частоте). Основная частота меняется в пределах 45...55Гц.

Здесь было предложено уже достаточно грамотных решений на основе преобразования Фурье. Я не буду повторять их. Также не буду говорить про затратные альтернативы типа MUSIC (мне самому стало интересно, как его упростить или построить новый алгоритм для случая кратных гармоник, может получится и не очень затратно). Позволю себе пофантазировать.
В описании сигнала ничего не сказано о шуме. Слукавлю и буду считать его несущественным (он таки может быть несущественным по сравнению с нелинейностью). Тогда достаточно сравнить мощность первой гармоники с мощностью остального сигнала без учета шума. Заметьте, знать точно частоту "основного тона" и мощности каждой гармоники в этом случае не нужно.

Потребуется три фильтра (не считая ФНЧ перед/после АЦП). Первый (простейший рекурсивный без умножителей) вырезает "постоянку" (смещение АЦП), не трогая первую гармонику и выше. Второй может быть вычислительно совмещен с третим (зависит от реализации), основное назначение - разделить сигнал на две полосы: ниже 55 и выше 90 Гц. Как вариант, при независимой реализации ФНЧ/ФВЧ, выход КИХ ФНЧ можно считать сразу децимированным (опять же для экономии). Далее усредняем мощность с выходов двух фильтров и делим. Либо, что экономнее, итеративно находим усреднение отношения двух мощностей без деления. Ценой одного умножения этим же методом можно получить и квадратный корень отношения двух мощностей. Либо, если хочется иметь значения в dB, то вычислить разность двух логарифмов. Возможно, потребуется согласовать задержки фильтров, чтобы скомпенсировать ошибки при изменении амплитуд. Естественно, точность измерения будет зависеть от избирательности фильтров (это верно и при выборе оконной функции при использовании преобразования Фурье), неравномерности в полосе пропускания, длины окна усреднения, точности арифметики. Ну и конечно же от шума, пренебрежение которым было разрешающим для подхода. При непрерывной (поотсчетной) обработке использование рекурсивных фильтров должно быть значительно экономнее по затратам, однако... однако не буду отвлекаться на вопросы реализации рекурсивных фильтров на малобитной арифметике (MCU в сопутствующей теме назван слабым).

Фантазирую дальше. Первый ФВЧ, удаляющий постоянную составляющую, можно совместить с последующим ФНЧ, что открывает возможность для добавления Гильбертовой пары для обоих полос, и от непрерывной обработки можно перейти к блочной. Для неперекрывающихся окон затраты, не считая накладных расходов, окажутся: 4 умножения и 4 сложения на отсчет для фильтров, 4 умножения, 2 сложения и одно деление (или два логарифма) на блок для вычисления мощностей и их отношения. При этом затраты памяти - одна ячейка для текущего отсчета и 4 аккумулятора по одному на каждый фильтр. Если же сначала сохранять в память отсчеты на всю длину фильтров, скажем N, то потребуется 2 умножения вместо 4 на каждый отсчет (за счет симметрии коэффициентов фильтров), но больше накладных расходов на организацию вычислений - сдается мне, экономии не получится, наоборот, зависит от процессора. При равной длине фильтров N, их коэффициенты займут в ROM 2N ячеек.

Как ни странно, в итоге получилось именно оконное взвешивание, как в заголовке темы. Такое взвешивание оказывается достаточно самостоятельным методом, а не только сопутствующим ДПФ. Для тех, кому все-таки странно, могу сказать даже, что это ДПФ сопутствует оконному взвешиванию для более эффективного с точки зрения экономии вычислений (и унифицированного!) формирования требуемых окон. Не имея необходимости во множестве унифицированных окон (часто так и бывает), достаточно сформировать необходимые окна заранее (а не в процессе обработки) и "взвешивать" ими без использования ДПФ. На чем и сэкономить... если получится

[Alexey Lukin]

Таким образом вычисляется THD+N. Если нет шума, то оно будет равно THD.

[Signal]

Фантазирую дальше. Первый ФВЧ, удаляющий постоянную составляющую, можно совместить с последующим ФНЧ, что открывает возможность для добавления Гильбертовой пары для обоих полос, и от непрерывной обработки можно перейти к блочной. Для неперекрывающихся окон...

Вот тут я зафантазировался. Такой номер пройдет только с основным тоном - он один в полосе. Для суммы гармоник во второй полосе огибающая не будет равна сумме огибающих (чего бы хотелось). А потому такое кардинальное децимирование, увы, даст некорректный результат. Прошу прощения, "Остапа несло" .