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

[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 гц основной гармоники) бина.

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

[Alexey Lukin]

Гармоники в спектре ищутся как локальные максимумы: a[k-1] < a[k] >= a[k+1].
Выбираете несколько главных максимумов, превышающих по амплитуде уровень шума, — это и будут ваши гармоники.

А зачем вам flat top окно? Если вы считаете суммы квадратов, то оно только транжирит ваше частотное разрешение. Берите что-нибудь попроще и пошире, типа Blackman-Harris.

[Zelepuk]

Я так и хотел делать - разбивать спектр на группы отсчётов и искать максимумы.
Дело как раз в том, что локальные максимумы находятся на разных расстояниях друг от друга. И при изменении частоты основной гармоники в широких пределах (45-55Гц) получается что группа для 8 гармоники (для 55Гц) "залезает" на группу для 9-й (для 45гц).

Остаётся только каким -то образом искать частоту основной гармоники и искать исходя из её значения. Или есть другой способ?


p.s. после применения окна Блэкмена-Харриса вижу, что точность определения гармоник ниже, чем прри использовании окна с плоской вершиной (flattop). Или нужно как-то просуммировать бины вблизи максимумов?

Сообщение отредактировал Zelepuk - Aug 11 2011, 05:41

[bahurin]

Вообще надо исходить из того с какой точностью вам надо померить уровни кратных гармоник и с каким динамическим диапазоном.

Вот тут написано как правильно выбирать окно для спектрального анализа. В вашем случае думаю стоит во первых увеличить выборку FFT, тогда широкие окна перестанут перекрываться. Во вторых вам надо брать локальные максимумы ваших гармоник (ничего вблизи этого максимума не складывать) и соотносить с максимальным значением на частоте 50 Гц, который тоже меряется как локальный максимум вблизи 50 Гц.

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

Сообщение отредактировал bahurin - Aug 11 2011, 08:59

[_4afc_]

Остаётся только каким -то образом искать частоту основной гармоники и искать исходя из её значения. Или есть другой способ?

Я делал так для звука... Может поможет?

Эскизы прикрепленных изображений

 

Прикрепленные файлы

 FFT1_new.m.txt ( 3.75 килобайт ) Кол-во скачиваний: 188

 

[Alexey Lukin]

Я так и хотел делать - разбивать спектр на группы отсчётов и искать максимумы.

Не надо разбивать на группы. Просто ищите максимумы.

p.s. после применения окна Блэкмена-Харриса вижу, что точность определения гармоник ниже, чем прри использовании окна с плоской вершиной (flattop). Или нужно как-то просуммировать бины вблизи максимумов?

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

[sup-sup]

Можно сделать FIR фильтр в частотной области. Для этого нужен шаблон с основной частотой и гармониками (чем больше - тем точнее), который нужно 'сжимать - разжимать' в пределах диапазона частот 45-55 Гц. Размер и тип окна выбрать, чтобы только старшие частоты не сливались.

[bahurin]

Не надо разбивать на группы. Просто ищите максимумы.


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

1. На группы нужно разбивать. Если частота первой основной гармоники может находится в пределах 45-55 Гц, то после взятия FFT необходимо искать локальный максимум в диапазоне 45-55 Гц. Никакие бины не складывать, просто методом перебора найти максимальный. Найдя этот локальный максимум надо запомнить его уровень A0 и его частоту f0. Для увеличения точности оценки максимума и частоты можно по трем точкам уточнить по полиному второй степени.

2. После того как мы оценили частоту основной гармоники, то предсказываем появление второй гармоники на частоте 2*f0 +- df, где df есть диапазон выбора второго локального максимума определяется точностью оценки частоты первой гармоники. Например если у вас частота дискретизации 4096 Гц, и вы берете 1024 точки, то точность оценки первой гармоники будет +-2 Гц. Тогда можно для второй гармоники взять df = 6 Гц (ошибка в 2 Гц также удвоится и будет +-4, с запасом можно взять +-6). Аналогично п.1 ищем локальный максимум, ничего не суммируя. Получили уровень второй гармоники равный A1 и частота второй гармоники f1. ВАЖНО частота второй гармоники также определена с точностью +- 2 Гц.

3. Предсказываем частоту третьей гармоники равную f2 = 3*f1/2. Поскольку мы вторую гармонику оценили с точностью +-2 Гц, то df для третьей гармоники можно также взять уже равной 4 Гц, поскольку df = 3*(+-2)/2 = +-3, а мы возьмем +-4, с запасом. Снова ищем локальный максимум запоминаем его уровень A3 и частоту f3

4. Для четвертой гармоники предсказываем частоту ее появления на f4 = 4*f3/3. далее везде полосу для поиска локального максимума можно брать +-4 Гц.

Итак ничего не надо суммировать и т.д., нужно лишь запусить итерационную процедуру поиска локальных максимумов в узкой полосе предсказывая каждый раз очередную гармонику. Соотвественно потом надо поделить полученные уровни An на A0.

[Zelepuk]

Я эспериментирую в CBuilder над FFT и окнами.
Нашёл реализуцию fixed-point FFT http://www.jjj.de/fft/fix_fft.tar.gz

Генерирую в Матлаб:
- сигнал (чистая синусоида частотой 53 гц - чтобы увидеть растекание спектра) длиной 1024;
- окно с плоской вершиной длиной 1024;

Все массивы в матлаб перевожу в int16 и копирую в CBuilder.

Делаю комплексное фурье(с использование подпрограммы на ссылке выше). Нахожу амплитуду (корень квадратный из суммы квадратов мнимой и действительной части поэлементно соответственно).

В итоге получаю через printf нечто вроде

0
0
0
16
1256
3256
1245
145
11
0
1
0
1
0
0
0
1
0

(все числа вымышлены и только для наглядности).

Когда суммирую числа из кучки выше (16+1256+3256+1245+145+11) и умножаю на 2 - получается в точности то, что и в матлаб!
Даже удивился точности!
Значит всё работает подумал я!

Да не тут то было.

Как только к сигналу примешиваю 5-ю гармонику с амплитудой 0.5 от основной, то сразу же спектр, выдаваемый целочисленным алгоритмом в CBuilder сильно искажаться по всем гармоникам, в итоге получаю кучу помех в спектре, причём там где спектр должен быть нулевым.

например

0
10
1
1
16
32
1114
123
0
1
0
1
10
32
556
124
11
0
1
0
1
15
123
112
144
114
10
1
1



Но в матлаб всё работает превосходно.

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

Почему так может происходить?
Может дело в реализации FFT? Может кто поделится другой реализацией (только для int16). На плавучку очень не хочется переходить (код будет перенесйн на микроконтроллер без FPU)

Сообщение отредактировал Zelepuk - Aug 12 2011, 09:18

[bahurin]

Когда суммирую числа из кучки выше (16+1256+3256+1245+145+11) и умножаю на 2 - получается в точности то, что и в матлаб!
Даже удивился точности!
Значит всё работает подумал я!

.

меня сильно настораживает тот факт что вы так спокойно решили складывать и умножать подгоняя свои результаты под матлабовские. Миром правят нормировки. По счастливой случайности складывая "кучку" вы получили одну и туже нормировку с матлабом. В другом случае у вас не проканало. Может стоит сначала отмоделировать все в матлабе и потом пошагово искать и устранять разночтения в свое программе

Сообщение отредактировал bahurin - Aug 12 2011, 10:08

[Zelepuk]

к сожалению с целыми числами полноценно работать в матлабе не научился.

Можно подробнее про нормировку?

[bahurin]

к сожалению с целыми числами полноценно работать в матлабе не научился.

Можно подробнее про нормировку?

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

[ivan219]

Zelepuk а как вы получаете гармоники? Из нулевой гармоники.

[Zelepuk]

Zelepuk а как вы получаете гармоники? Из нулевой гармоники.

всмысле как? Я просто умножаю на окно, делаю БПФ и затем смотрю что получилось. Нулевая гармоника - это постоянная составляющая.

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

Дело как раз получается в том, что я делаю всегда одно и тоже с сигналом: сигнал+окно-> БПФ -> суммирование на некотором интервале (именно суммирую, я просто складываю "палки" и всё)

Если сигнал есть одна гармоника (пробовал разные частоты от 45 до 55Гц) - всё красиво и точность по амплитуде 0,01%, если сигнал представляет собой сумму гармоник (пусть есть 7-я и 5-я гармоники с амплитудами 0,1 от основной) - наблюдается погрешность.

По вашим словам получается когда сигнал чистая синусоида, то мой алгоритм учитывает ослабление от окна, когда сумма гармоник - не учитывает. Но я всегда делаю с сигналом одни и теже операции.

Сообщение отредактировал Zelepuk - Aug 12 2011, 12:02

[bahurin]

По вашим словам получается когда сигнал чистая синусоида, то мой алгоритм учитывает ослабление от окна, когда сумма гармоник - не учитывает. Но я всегда делаю с сигналом одни и теже операции.

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