Столкнулся с такой задачей: Есть АЦП с частотой дискретизации 8 кГц - задан жестко. Есть микроконтроллер, работающий на частоте 120 МГц - тоже задан жестко. Частота сигнала - частота сети 50 Гц. Т.е. после АЦП на микроконтроллер идут 160 отсчетов на 1 период сигнала. Далее на 10 периодах применяется БПФ на 1024 отсчета (тут БПФ не рассматривается, уже реализовано).

Собственно в чем весь вопрос: Какой именно вид интерполяции применить лучше всего (есть линейная, квадратичная, кубическими сплайнами, полиномы Лагранжа и прочее)? Ведь необходимо из 1600 отсчетов сделать 1024 - по-любому необходима интерполяция. И еще необходимо обеспечить точность - в идеале около 1%.

И еще один момент - всё это делается на микроконтроллере, т.е. потом все эти алгоритмы будут преобразованы в код на C++.

Нет такой необходимости.

Берёте 1600 отсчетов, умножаете их на оконную функцию (исходя из N = 1600). Дополняете массив нулями до размера 2048. Дальше БПФ.

1. Зачем Вам интерполяция на 10 периодах? Делайте на одном.
2. Почему нужно переходить от 160 отсчетов к 128 (к примеру)?
3. Цель интерполяции?
5. Погрешности по видам интерполяции считаются или моделируются - это не большая сложность.

Всегда думал, что интерполяция - это вычисление новых отсчетов, но преобразование 1600 в 1024 отсчета - это вроде как уменьшение числа отсчетов, разве нет? Тогда это уже не интерполяция, а децимация.

В данном случае топикстартер желает изменить частоту дискретизации. Для этого потребуется и то и другое.
Тоже не вижу необходимости. Делать БПФ по 1024 отсчетам. А что частоты не "круглые" получаются, так не великая беда. Может, результаты интерполировать, чтобы красивые графики рисовать?

Да, но если добивать нулями, то сильно искажается спектр и тогда оценить уровень гармоник не представляется возможным. К тому же БПФ на 2048 точек уже становится слишком ресурсоемким, даже на Cortex-M4F. Но за идею спасибо.




1. Ибо БПФ согласно ГОСТ нужно считать именно на 10 периодах. На самом деле не БПФ, а гармоники нужно считать на 10 периодах, поэтому такое условие.
2. Вот и получается, что необходима передискретизация с 1600 отсчетов на 1024 отсчета.
3. Вот как раз и цель интерполяции - пункт 2.
4. А погрешность - согласен, можно легко промодилеровать.

Вот мне и хотелось узнать как это делается сейчас и какие есть варианты решения проблемы.




Если немного уйти от темы, то в общем случае на вход АЦП идет ослабленный сигнал электрической сети. АЦП на 8 кГц, как я уже говорил. Данные с АЦП приходят в микроконтроллер. И задача стоит в том, чтобы оценить уровень гармоник вплоть до 50. И вот одно из решение: Так как частота сети плавает, то применить интерполяцию и получить 1024 отсчета и далее БПФ и вот они гармоники.

Полифазный фильтр

Поправьте меня, если я что-то не правильно говорю:

Полифазный фильтр действует примерно так:

Сначала мы увеличиваем частоту дискретизации в М раз - т.е. интерполируем сигнал.
Затем прореживаем до нужного нам количества отсчетов, т.е. откидываем ненужные, т.е. производим децимацию сигнала?

Не искажается ли спектр после отбрасывания отсчетов?

Конечно искажается. Чтобы этих искажений не было в итоговом сигнале в полифазном фильтре (неожиданно) есть фильтр, который защищает от наложений.

Тогда совсем глупый вопрос: И всё это дело можно организовать программным образом на микроконтроллере? Если да, можете подсказать пару методов или алгоритмов. Пожалуйста.

лагранжем проще по сравнению с кубическими сплайнами, не надо по всем данным от начала до конца бегать туда/сюда чтобы граничные условия выполнить, можно на лету делать.
http://ru.dsplib.org/content/resampling_lagrange.html
но для 1% при 160 точках на период и линейная интерполяция возможно вполне сгодится, хотя для 20й гармоники останется только 8 точек на период что уже маловато будет.

ну и добивание сигнала нулями спектр вообще-то не портит, от того что 0 умножают на синус/косинус и интегрируют результат не меняется.
сетка частот съедет, но и так у вас 8кГц и 1024 отчёта = 7.8125Гц в 50Гц не особо попадают.
https://electronix.ru/forum/index.php?s=&am...t&p=1452224
если не хочется делать по 2048 точкам, берите 6 периодов и добивайте нулями до 1024.

про БПФ: вместо всего спектра и 1024 отсчётов, если нужны только гармоники 50Гц, то ДПФ только для этих нескольких частот может оказаться быстрее.
БПФ грубо говоря 100 тактов на отсчёт - 100 000 тактов для вычисления спектра. + оконная функция + передискретизация, но это мелочи по сравнению с самим БПФ, ещё %10 сверху.
а ДПФ, особенно если Герцелем без тригонометрии - несколько тактов на отсчёт. до ~20 гармоник выигрыша от БПФ нет, плюс считать можно налету без буферизации данных вообще.

Можно. Сам правда не делал (в смысле на МК не делал). Точно есть в Лайонсе и что-то было в книге Айчифер, Джервис "Цифровая обработка сигналов: практический подход".

Здесь, почитайте, для начала.
DSPLIB

Как раз на 10 периодах и получается максимум до 50 гармоники - я это уже моделировал и код писал в IAR - всё работает, БПФ на 1024 считает четко, погрешность очень маленькая и как раз до 50 гармоники.

Ведь частота сети плавает, хоть это и считается медленно меняющимся параметром, в связи с чем меняется и количество отсчетов на 10 периодов (Почему именно 10? ГОСТ 30804.4.30 - 2013). 1600 отсчетов - идеализированный случай. Будь частота 48, отсчетов будет уже 1667 на 10 периодов. К чему я это - получить 1024 из плавающих 1600 - сможет ли Герцель с таким справиться? И уложиться в 2.5 - 3 мс.

сверху гармоники ограничены частотой дисктеризации, а не временем измерения.
и еще раз, 8кГц и 1024 точки даёт разрешение по частоте 7.8125 Гц, что в 50Гц нацело не очень попадает и для первых гармоник, значительную часть от 1% ошибки можно в это месте легко собрать.


Герцель ничего не интерполирует, это и есть дискретное (обычное, не быстрое) преборазование Фурье (F(w)=\int{f(t)*cos[sin](wt)dt}), где вычисление тригонометрии(sin/cos) заменено цифровым осциллятором (одно умножение и сложение), но его, в отличии от БПФ который даёт весь спектр сразу, не обязательно считать для всех 1024 частот, можно посчитать только на тех 50ти частотах которые интересны, ну и длина выборки любая, степенью 2 не ограничена.

10 периодов это 200мс, а вычислить надо за 3мс?
120МГц и 3мс это 360000 тактов, или 350 тактов на отсчёт для 1024 точек.
или 7 тактов на отсчёт если считать для 50 частот.
Герцелю надо умножение, пару сложений, на отсчёт для одной частоты, в семь тактов если и уложится, то на пределе.
так что для 50 и больше частот БПФ, который выдаёт спектр целиком будет пожалуй побыстрее.

Разложение числа на простые множители:
1600 = 2 · 2 · 2 · 2 · 2 · 2 · 5 · 5
Микс БПФ по двойке и пятерке.

+1.
Подробнее можно посмотреть здесь: Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов.

О вариантах БПФ
Parallel.ru

Еще:
Морфизприбор - БПФ для ГАС:
Нажмите для просмотра прикрепленного файла
К слову - нетрадиционные БПФ, это сейчас - нетрадиционные. Тогда, в 70-х традиционными были скорее нетрадиционные сегодня, т.к. цена быстродействия за лишнюю точку была высока.

P.S.
Тогда, когда еще не все были повально увлечены БПФ по степеням двойки, достаточно интересным и более быстрым был алгоритм Винограда. Представляя одномерное БПФ как многомерное, с разложением по простым множителям, удавалось значительно повысить скорость обработки.
Нажмите для просмотра прикрепленного файла

Годная статья но маленькая, сейчас "нетрадиционное ффт" называется mixed radix fft. Вот например


А вот такой вопрос топикстартеру, скажем первая гармоника (50 Гц) изменяется на 0.1 %, тогда 50я гармоника изменится примерно на пол бина FFT и чего Вы тогда намеряете?

Для ТС-ра, даже и это - слишком много, надеюсь, что - пока.

P.S.
Родной мой, ну ты хватил - давай еще CRAY и пр. сюда подтащим.
ТС-ру нужна реализация на микроконтроллере, какие уж тут ядра и параллельности.

В 4G, он же LTE, есть размерность FFT, равная 1536, так что там 3^1*2^9. И сейчас есть нетрадиционные размерности

Не пугайте юношу страшным словом "Метрология".

Ага для CRAY есть аналогичное . На самом деле ценность математики, в этих статьях - это возможность записать любое FFT в виде
одной формулы, а там CRAY или Intel I-7 или АRM это уже вторично.

Для тех, кто в такой теме давно - может и да.
Но, зачем ТС-ра пугать в его теме?

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

Задайте частоту дискретизации АЦП 6400 Гц. За 128 выборок наберете 20 мс. Вычислив БПФ, получите 64 гармоники 50 Гц сигнала, с полосой каждого бина 50 Гц.

Потом для большей точности измерения амплитуды усредняйте результаты БПФ.

Не можете уйти от 8 кВыб/с? Интерполируйте на 4, децимируйте на 5. Хотя я этого садомазо из-за упертости схемотехника (или кого?) не понимаю.

Да, АЧХ каждого бина не плоская, на краях загибается. Думаю, можно забить. Или взять 256 выборок, БПФ по 256, получить шаг 25 Гц. Лучше 512 или больше, чтобы знать, как относить палки спектра к нужным гармоникам сигнала.

10 раз усреднённый спектр по 20мс и один по 200мс совсем не одно и то же.

А то!

По-хорошему, я бы вообще делал десять замеров периодов, начиная каждый замер в момент прерывания по zero-detected.
Мало ли там, что - по ГОСТУ-у

Я бы вообще не считал (количество замеров)... Усреднял бы и усреднял непрерывно замеры.

Говорят, по ГОСТ-у положено 10. Так, что после 10-ти и после 11-го идет просто задержка выдачи результата на период.

Вообще алгоритмы измерения сетевых гармоник подробно изложены в ГОСТ 30804.4.7-2013 и там в качестве опорного метода рекомендуется БПФ.

http://meganorm.ru/Data2/1/4293778/4293778588.pdf



Проблему ТС я вижу в том, что у него частота дискретизации не синхронизируется с сетевой частотой, как рекомендует ГОСТ,
а является фиксированной 8 кГц. Скорее всего используется готовый AFE. Чтобы обеспечить работу в требуемом по ГОСТУ интервале измерений от 42 до 58 Гц сетевой частоты, на 10 периодах сети имеем возможный разброс от 1380 до 1900 отсчетов и отсюда проблему передискретизации в фиксированные 1024 из-за дробного и непостоянного коэффициента.

Интерполяцию можно делать через кубические сплайны, можно через интерполирующий фильтр. Здесь что будет быстрее не знаю, но я бы через сплайны делал.

Да, всё именно так. Я немного ступил и задал неверные начальные данные, из-за этого обсуждение свернуло немного в другую сторону. То, что Вы описали - именно это мне и надо получить. В ходе обсуждения и было предложено делать интерполяцию, чтобы получить фиксированные 1024 отсчета и как-то решить проблему с плавающими входными отсчетами. Мне как раз и хотелось узнать - какую интерполяцию применяют в данном случае? Точно не линейная - точность требуемую по ГОСТ не обеспечим. Сплайнами? Точность первоклассная, а вот ресурсов требует наравне с БПФ (примерно 70, а то и все 90 процентов). Делать сначала интерполяцию, а потом децимацию лишено всякого смысла.

Спасибо!

А учебники вы из принципа не читаете? Сначала интерполяция а потом децимация и никак иначе. В противном случае вылезут ложные частоты.

В вашем случае хватит линейной, но по правилам делаем Sinc. Хватит 3-х или 5-ти точечного ядра. Интерполяцию делаете интерполирующим фильтром.

Поправьте меня, если я не прав. Все операции передискретизации и преобразование Фурье делается на микроконтроллере, CORTEX-M4F. Вы предлагаете исходные отсчеты с АЦП (а их от 1380 до 1900 ), приходящие на МК сначала интерполировать в N раз, а потом провести децимацию (отбросить каждый M отсчет) - это, как я понимаю и есть интерполирующий фильтр. А интерполяцию делать Sinc-функцией. Но ведь таким же образом можно сделать и сплайн-интерполяцию, подозреваю, что и примерно с такой же скоростью. Да, интерполяция Sinc-функцией очень точна, но вместе с тем делается гораздо дольше. Линейная интерполяция не обеспечит необходимой точности, заданной ГОСТ - уже моделировал в MathCAD - расхождения очень большие.

Прокомментируйте мое предложение:
https://electronix.ru/forum/index.php?showt...t&p=1459238
а то складывается впечатление, что вы его пропустили.

.

Восстановление значения сигнала между двумя отсчетами (интерполяция) автоматически означает повышение частоты
дискретизации в N раз. Чем ближе восстановленный отсчет к первому исходному отсчету - тем выше N. В пределе для непрерывного
сигнала частота стремиться к бесконечности.
Существуют важные частные случаи передискретизации, для которых вычисление новых отсчётов производится проще:
- децимация с целым коэффициентом (уменьшение частоты дискретизации в целое D число раз);
- интерполяция с целым коэффициентом (увеличение частоты дискретизации в целое N число раз);
- изменение частоты дискретизации в рациональное (N/D) число раз (этот случай можно рассматривать как комбинацию двух предыдущих).

https://ru.wikipedia.org/wiki/Передискретизация

https://blog.amartynov.ru/интерполяция-или-фильтр


Сетевая частота может сильно меняться и поэтому задавать фиксированную частоту нельзя. Ее надо синхронизировать с сетевой из-за требования обеспечить интервал измерений строго 10 периодов сетевой частоты (сетевая частота измеряется другим методом и заранее известна). Если нельзя уйти от 8 кВыб/с то коэффициент N/D (интерполяция/децимация) будет переменным и очень дробным, потребует фильтрации на очень высокой частоте. Думаю надо смотреть в сторону кубических сплайнов https://ru.wikipedia.org/wiki/интерполяция
ну, может быть, полифазного фильтра. Поправьте, если не прав.

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

Если говорить более общо, то меня интересовало то, каким способом сейчас делают интерполяцию применительно к моему случаю (задаче). Будь частота микроконтроллера выше и составляла хотя бы 200 МГц, то можно смело делать сплайн, лишь оптимизировав код под конкретный процессор и используя все его (микроконтроллера) функции и ресурсы, и не бояться, что пока будет обрабатываться сигнал мы не успеем провести другие измерения и какие-то другие действия.

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

Полоса каждого бина БПФ будет 50 Гц, влезет любое отклонение частоты от номинальной.

Если уж так надо привязаться к частоте сети, то... надо к ней привязаться. ФАПЧ!

и для 20й гармоники тоже?

не надо к ней привязываться ФАПЧем, её надо просто автокорелляцией измерять, а потом делать обычное преобразование Фурье для нужного количества гармоник по сетке частот, соответствующей измеренной частоте основной гармоники.

Впритык. Частота сети контролируется глобально. Шаг влево-вправо - расстрел. Так меня учили.
При отклонении основной гармоники на 1 Гц 20-я уйдет на 20 Гц, а 50-я на 50 Гц. Влазит.
А она кому-то нужна, аж 20-я гармоника?

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

Дополню. Интерполировать синус кубическим полиномом (все эти сплайны-лагранжи-фарроу-эрмиты-катмулы-с-роммами...), а потом вычислять в нем 20-ю гармонику? О_о. Их там будет немеряно всяких, в-основном, несуществующих в исходном сигнале.

А, по сетке частот Герцелем? Можно.

Как раз таки ГОСТ и требует измерение 50 гармоник сетевой частоты и от ГОСТа никуда не уйти. ГОСТ 30804.4.7-2013.

Интерполяция сплайном, по крайней мере как показывает моделирование в MathCad, не очень сильно искажает исходный сигнал. Порядка 0,5 В.

Герцелем вычислять все 50 гармоник? При условии, что количество отсчетов всегда плавает. Можно ли будет уложиться хотя бы в 5-7 мс? На CORTEX-M4F, 120 МГц.

А почему "уложиться хотя бы в 5-7 мс"? Почему не в 0,5-0,7 мс?
И почему именно "CORTEX-M4F, 120 МГц"?
В природе есть много процессоров, хороших и разных..

FFT надо ~100 таков на отсчёт, Герцелю - одно умножение и пара сложений, соответственно для ~20 гармоник Герцель будет быстрее чем FFT для всего спектра.
ну пусть даже будет 10 тактов на отсчёт для Герцеля, для 10 периодов 50Гц на 8кГц это 1600 точек * 10 тактов на отсчёт * 50 гармоник / 120МГц = 6мс.
плавающее количество отсчётов Герцелю никак не мешает и не надо ничего интерполировать.
то что последняя точка может не попасть в целое количество периодов на пол периода АЦП, +-0.5/8кГц = 60мкс, для 10ти периодов думаю никак на точность не повлияет.

Ну, это зависит от процессора. Для 16-ти битного FFT на 2048 точек blackfin'у требуется ~16 тактов на отсчет..
Кроме того, есть процессоры, в которые встроен аппаратный ускоритель FFT.

ну так то блэкфин, там и два алу, и аппаратные циклы, и куча шин, что позволяет за один такт и два МАСа сделать и ещё и в память слазить за следующими данными.

а у ТС - кортекс М4.

1024-FFT (Complex in Q15 Format)
Cortex-M0 : 855 733 cycles
Cortex-M0+ : 664 531 cycles
Cortex-M3 : 204 244 cycles
Cortex-M4 : 89 839 cycles

хотя для М4F попадались цифры в 40 тактов на отсчёт с плавающей запятой.
FFTN M4 A8 A9 A15 BF5xx BF70x SHARC21489
64 3709 3773 3358 2264 2200 1526 783
128 9811 6384 5682 3830 5249 3431 1334
256 21575 11114 9891 6668 11744 7611 2542
512 37813 21852 19448 13111 27385 17084 5189
1024 96630 50738 45157 30443 60216 37568 10972

так что 100 тактов на отсчёт для оценки вполне пойдёт.

ну и для анализа сетевых 50Гц это всё-таки постараться надо чтобы в вычислительную мощность процессора упиралось.

Кстати не факт что 20 "Гертцелей" будут быстрее одного FFT . Скажем есть 1024 отсчета для FFT требуется порядка N * log2 ( N ) = 10240 комплексных операций. Считаем что для Гертцеля требуется одна комплексная операция на отсчет , итого 20 Гертцелей будут стоить 20*1024 = 20480 операций.
Тут кто-то про ФАПЧ выше писал... Такая идея... Хватаемся ФАПЧ за первую гармонику, из генератора фапч делаем опорные комплексные экспоненты с частотами 2*f1, 3*f1.... 20*f1, считаем их корреляцию со входным сигналом.

-----------------------------------
p.s. Кажется лет 6 назад была точно такая же тема здесь. И вроде было все подробно разобрано, и кажется пришли к мысли что передискретизация(интерполяция) особо не нужна чтобы измерить частоту (амплитуду) отдельно стоящей гармоники.

Это дежавю каждый год, если не чаще.
Вообще есть специальные flat top окна, читай FIR фильтр с плоской АЧХ в полосе пропускания, чтобы не вносить ошибку в амплитуду из-за частотного сдвига, точного значения частоты не нужно, только грубое, чтобы выбрать нужный фильтр в банке, куда гармоника попала, громадное FFT тоже не нужно, поскольку АЧХ фильтров в банке будут сильно перекрываться, поэтому можно проредить и использовать эффективную реализацию банка фильтров в виде полифазник+небольшое FFT.

не все комплексные операции одинаково полезны (с).
вон на 1024FFT кортексу М0 надо в десять раз больше тактов чем М4.
а вот для того чтобы умножить и сложить пару чисел ему точно в десять раз больше тактов не надо.


и чем это от Герцеля отличается?
только частоту основной гармоники просто измерить, автокорреляцией или через среднеквадратичное значение производной, а не из фапча. и потом считать Герцелем те же интегралы Фурье для частот 2*f1, 3*f1.... 20*f1,.