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

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

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

[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