Хочу минимизировать эфект "растекания" спектра. Сигнал содержит 40 гармоник кратных основной частоте.
Известен интервал изменения основной гармоники (50+ 1Гц)
Частота дискретизации задана жёстко 16000Гц.
По расчётам кол-во точек сигнала равно N = 8192.
Затем по некоторому алгоритму получаем частоту основной гармоники (пусть частота равна fосн. = 50.55Гц).
Далее, по идее, для минимизации "растекания" нужно интерполировать сигнал с учётом полученной частоты основной гармоники.

Расчётные данные показали, что нужно пересчитать на частоту fd = (327.68*fосн) = 16564,224 Гц.
________________________________________________
Вот здесь самое интересное начинается!!!

Чтобы идеализировать интерполяцию просто изначально задаю частоту 16564,224 Гц на том же интервале времени, который равен интервалу времени выборки 8192 отчётов на частоте 16000Гц (т.е. (1/16000)*8192 = 0,512 с.).
Другими словами окно временное оставляю тем же, что и для 16000Гц.

Конечно же в спектре сигнала наблюдается сильное искажение...

Вот такой вопрос:

Я сделал вывод что при интерполяции каким-то образом меняется величина временного окна... иначе никак((
КАК ЭТО ДЕЛАЕТСЯ????

Встал на этом вопросе прирасчётах, очень нужна помощь...

Заранее спасибо всем ответившим!

Такого эффекта на сегодняшний день наука не знает...



угу. Вы хотите сделать передискретизацию так, чтоб гармоника была кратна бину дпф.



Гы. Вы повышаете частоту дискретизации. На интервале 0.512 c получится больше отсчетов.
В Вашем случае их станет 8480,882688. Никакой метафизики тут нету...



Что именно? Как вычислить дпф на 8к точек по 8480,882688-и точкам? Дык вариантов много.
Можно вычислить несколько дпф на 8к точек со смещением по выборке из 8480, а потом их усреднить. Есть еще более сложные способы.
Искажения спектра скорее всего связаны с горбатой интерполяцией.

ДПФ это банк фильтров с прямоугольной ИХ, АЧХ каждого фильтра медленно спадает и фильтры сильно перекрываются между собой, можно использовать банк фильтров с другой импульсной характеристикой фильтра прототипа так чтобы фильтры не перекрывались и имели большое подавление за пределами своей подполосы и тогда никакого "растекания" не будет.

А зачем так сложно, что требуется получить (измерить). Может быть просто взять окно какое-нибудь, наложить на выборку посчитать FFT и ничего не интерполировать?

Здесь предыистория:
http://electronix.ru/forum/index.php?showtopic=66966

Я так понимаю TigerSHARC хочет точно померить амплитуды всех 40-ка гармоник.

... ну и как это звучит по научному?



да



совершенно верно



так я же пишу, что здесь не интерполяция, пока, а просто частоту дискретизации сменил на том же интервале...

Вопрос в том что при интерполяции мы имеем дело с с некоторым количеством отчётов на некотором жестко заданом интервале (в частности 0.512с). А я думаю что при интерполяции длина окна должна меняться.
вобщем нужно менять длинну окна... как-то...

что бы не быть голословным приведу вырезку из диссертации...
на этой странице говорится об особенностях применения ДПФ для анализа 40 гармоник промышленной частоты (50Гц)...

Знает, в книге А.Б. Сергиенко "Цифровая обработка сигналов" посмотрите.

По теме: может, использовать общепринятые окна или метод MUSIC?

всем кто смотрит тему: "пожалуйста посмотрите файл с вырезкой из диссертации в моём посте выше"

Стандартное окно (Хемминга и др...) не применимо из-за погрешности порядка процентов.
Здесь строится прецизионная система.

Нужно использовать именно интерполяцию (предположительно фильтром Фарроу 3-го порядка)

Не читайте сергиенку. Книшки надо правильные читать, чтоб спектр не растекался.

ps
Прошу прощения за оффтоп...

как по научному-то??? у Сергиенко тоже "растекание" как "растекание" упоминается... нипойму вашей ирронии..

Где окно погрешности вносит ? В измерении амплитуды гармоник? Насчет интерполятора, Вы смотрели какая будет АЧХ у полиномиального интерполятора 3го порядка? Наверное теже несколько процентов получится завал АЧХ на верхних частотах...

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

Для прецизионности системы нужна интерполяция синком

А что спектр надо оценивать не только для гармоник основной частоты, но и между ними?
Если нет - то зачем БПФ? Достаточно ДПФ и быстрее. Существует для любого окна N
Как говорил petrov - банк фильтров на кратных частотах.

Интерполяция меняет размер окна???

Разве на 40 гармоник ДПФ будет быстрее??

40 гармоник ДПФ - это примерно 80 * N действительных операций (40 для синуса, 40 для косинуса MAC)
FFT - это 2*N*logN = 26*N комплексных операций. Вы же вычисляете не 40, а 4096 гармоник,
из которых Вас интересует, возможно, только каждая 10-я. Зависит от постановки задачи.
Конечно, если Вам нужен красивый спектроанализатор тогда другое дело...

Если не быстрее, то примерно одинаково



Какая интерполяция?
Если Вы пересчитаете свои 8192 отсчета на другую дискретную сетку, то понятно, что размер окна в отсчетах изменится. А в физических единицах - нет. Ширина полосы фильтра, который можно получить по этим точкам в герцах не изменится

Ну так делайте частоту дискретизации такую чтобы на на 8K выборке помещалось целое число периодов основного тона. Хотите тактовую частоту АЦП меняйте, хотите ресамплинг делайте. Только как-то довольно сомнительно это. ИМХО туже информацию можно извлечь без передискретизации.

p.s. А в дисертациях не обязательно умные вещи пишут...

Укажите какие.
О квалификации Сергиеко А.Б. не Вам судить. Лично с ним знаком.

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

По теме:
Ознакомился с бумажкой. Потом ознакомился с гостом и приложением. И не очень понял, в чем же проблема? Сделаете Вы передискретизацию, выполните дпф на нужное число точек. Все нужные частоты будут точно совпадать с частотами бинов. По сути вы и измените длину окна. Окно 0.512 c на 16000Гц станет 0.4945...с на 16564,224Гц.

А как объяснить такое.
Читаю:
"... С учётом близкого расположения гармоник на частотной оси ( через каждые 50 Гц) и малых отклонений частоты основной: гармоники (до сотых долей Гц) от номинала требуется высокое разрешение БПФ (вплоть до 0.03 Гц). Таким образом, размер БПФ (длина выборки) превышает 533 тысячи отсчетов прни f=16 кГц. но при выполнении условна отсутствия искажений размер БПФ можно существенно уменьшить - до порвлка 10 тысяч отсчетов при той жe частоте дискретизации)"

Получается не соответствует действительности?

P.S. условие отсутсвия искажений fd*w = N*fосн.
fd - частота дискретизации;
w - некоторое целое число;
N - размер БПФ;
fосн. - частота основной гармоники (49-51Гц)



наконец-то!!! спасибо это и хотел услышать... просто в MATLAB не додумаюсь как это сделать... - изменить длинну окна...
как-то нужно время пересчитать наверное...
мы же изначально на интерполятор пускаем сигнал из отсчётов в окне 0.512 с частотой 16000, преобразую в сигнал с частотой 16564,224 - получается больше отчётов на этом же окне(((
наверное сам принцип "интерполяции" не понимаю...

Рабинер Голд Теория и применение цос. Годится не только как базовый учебник, но и как справочник.
Правда, матлабских скриптов там нет. Эт, канешна, серьезный минус...




Нет. О квалификации Сергиенко А.Б. я судить не могу и не хочу. Но как читатель, вполне могу покритиковать его книгу.



Завидую.

Например, меняете частоту АЦП в соответствии и берёте всегда 8192 отсчёта. Тогда "условие отсутствия искажений"
будет выполняться и все гармоники точно лягут на бины ДПФ.

В случае интерполяции прийдётся или отбросить часть отсчётов или добавить нулей, чтобы получить примерно то же самое. Но трудно будет сделать хороший интерполятор. Адаптивный интерполятор Фарроу сам будет давать искажения, а чтобы сделать прецизионный ресемплер с фильтрацией Вам нужно будет пересчитывать фильтры под основную частоту в реальном времени...

Хотя может Фарроу и будет работать за счёт большого оверсамплинга, это вообще-то вопрос. 50гц*40=2 кгц
8 отсчетов период может мало для "прецизионной" интерполяции, но можно пробовать ресамплинг сделать на повышеной частоте, потом прорежать и делать ДПФ (или БПФ)

Да нет все наоборот, нужно получить 8K (== размер FFT) отчетов, из k*8K входных отсчетов (k- коэффициент вроде k=0.99.. 1.01). После ресамплера длина окна постоянна, а фактически, во времени длина окна будет кратна периоду основного тона.

всё верно! Так и надо делать... никакого аппаратного изменения частоты.
мне сам принцип изменения окна временнного при интерполяции интересует.
thermit сказал что окно меняется. Теперь ищу информацию, где доступно написано как именно.

Интерполяция в контексте цос - повышение частоты дискретизации. те отсчеты сигнала берутся чаще, и как следствие - их становится больше в единице времени.

Да не нужен Вам никакой ресемплер, наложите Гауссовское окно и посчитайте правильно. Получите очень даже прецизионно. И чему равно Ваше прецизионно? У нас прибор уже работает на этом принципе. Точность по напряжению 0.01% - это полностью, а не только точность вычислений, по частоте и фазе еще существенно выше.

мой алгоритм предполагает реализацию для приборов контроля качества электрической трёхфазной сети в реальном времени.
Интересно:
1) у вас - это у кого? (если не секрет)
2) 0.01% Это с чётом погрешности первичных датчиков?
3) какие параметры системы(чатсота дискретизации, длинна БПФ...)?

Единственное что приходит на ум при применении сглаживающих окон это введение како-то поправки...


Но как вы ответите на такой тезис:
"... несмотря на рекомендации, приведенные в ГОСТ 13109-97 (обязательное приложение Е6), по выбору типа и длины окна, четких правил выбора этих параметров не существует, А поскольку любое окно искажает спектр полезного сигнала» так как произведение функций во временной области соответствует свертке их образов в частотной области, то погрешность вычисления спектра зависит от длины окна и вида спектра анализируемого сигнала. Выбор типа окна должен основываться на предполагаемом спектральном составе сигналов, но спектральный состав сигналов от датчиков тока и напряжения может изменяться во времени, что определяется характером включенной нагрузки. Поэтому единственный правильный вариант - обеспечение условия отсутствия искажений при применении прямоугольного окна. Тем более что это не противоречит ГОСТ 13109-97"

Привет SHARK. По Твоей тематике можешь посмотреть работы VU Jiekang (Ziejiang Univ.) . Сейчас работает в Америке в энергетической компании. Может и прибор уже создал. Никакого БПФ! Я использую его метод, но в другой области.

Книга бесусловно хорошая, но есть проблемы с примерами.

Мы знакомы???)

Что за метод?
Где найти его материалы??? поиск ничего не дал(мне по крайней мере)...

Прошу прощения за фамильярность. Навалом ссылок.
Вот журнал: IEEE transactions on power delivery ISSN 0885-8977, 2005, vol. 20, no1, pp. 366-374 [9 page(s) (article)] (22 ref.)
Еще: http://www.springerlink.com/content/r113673647v40251/
Там и почта его есть.
Похоже это та статья. В бесплатном виде пока не нашел.
Мой экземпляр погиб вместе с жестким, на бумаге осталось. Может подскажете где эти transactions почитать подешевле можно?

а я откуда знаю... я с этими методами даже не знаком

Так попробуйте познакомиться со статейкой из IEEE. Кому нужно Вам или нам?
http://electronix.ru/forum/index.php?showt...&start=1215

Там как раз интерполяция по Лагранжу
Статейке в Springerlink из китайского журнала кажется защищенная от копирования, если я правильно понял, скопировать бесплатно её не получится

Вот эта?

В каком смысле - "Потом удалю"?

Вот та статья VU Jiekang из IEEE более соответствующая контексту оценки частоты силовой сети.
IEEE transactions on power delivery ISSN 0885-8977, 2005, vol. 20, no1, pp. 366-374 [9 page(s) (article)] (22 ref.)
Там нет однако оценок спектра
http://rapidshare.com/files/315859486/Powe...qEstimation.pdf

Большое спасибо за статью.
Я хотел бы попросить у TigerSHARK какую-нибудь реализацию с отклонением по частоте не менее 0.1 Гц. Я попробую определить нужные параметры и скажу как.
Сэмплы менее 24бит не интересны.
Эффект "растекания спектра" не относится ни к Вашим данным (согласно Котельникову Вы имеете все что нужно) ни к тому, что Вы желаете получить. Это эффект БПФ и окна,
Следовательно применив для реальных сигналов более прямой способ выделения известных гармоник - умножение на синус, на косинус, сглаживание... - Мы получим идеализированный невозможный для БПФ результат. (для ослабления эффекта Гиббса потребуются данные о предыстории и послеистории сигнала, всего на длину сглаживающего фильтра)
В случае реальных сигналов результат БПФ менее точный из-за неизбежной амплитудной и еще хуже незначительной частотной модуляции и из-за несинхронной части шумов. И чем больше число выборок тем хуже, уменьшается только шум.
Наверное существует оптимальный объем выборки для каждого сигнала, при котором шум достаточно уменьшен, а модуляции не сильно затронули соседние частоты.
Ресэмплинг мне не понравился. Самый лучший это подстановка sinx/x для некоторых значений выборки около требуемого значения времени - приводило к слишком длительным вычислениям, а применение Фэрроу слишком обогащало спектр Моего сигнала. А сигнал был принципиально несинхронным. В целом знание пред и послеистории значительно повышает точность ресэмплинга, делайте запас сразу.
Самое простое - попробовать. А вдруг у Вас ни шумов ни модуляций ни посторонних сигналов, а процессор Атом, тогда все должно быть шоколадно.
А еще мне встретилась задача где точность 10% = очень хорошо, а 6% = идеально. Удачи.

Готов предоставить всё, что в моих силах.
Немного не понял вопроса. Что вы имеете ввиду под "реализацией"...

Честно признаюсь, что Я - дилетант, и вокруг никто курс ЦОС не прослушал. Жаргон такой образовался: реализация = набор сэмплов полученных одновременно с нескольких датчиков при одном испытании, например "вчерашняя реапизация" или "реализация где датчик отвалился". Если у Вас есть правильный термин постараюсь выучить.
Постараюсь помочь, хотя времени мало - сегодня начальник выбрал корпус для прибора, с понедельника работа под давлением до макетного экземпляра прибора. Математика и фреймы написаны.
Мне нужно:
1. Какой АЦП.
2. Реализация(файл содержащий результаты АЦП без промежутков и заголовков)

... у меня никакой аппаратной реализации((
АЦП по расчётам 12 бит минимум.
Моделирую всё в MATLAB. Подаю на вход алгоритма полигармонический синус содержащий 40 гармоник кратных 50 Гц (50, 100, ... до 2000) с одинаковыми амплитудами и нулевыми фазами (это наихудший теоретический случай, позволяющий оценить погрешность сверху) Ставлю кратную частоту дискретизации (минимум 6400Гц для БПФ). Получается эталонный спектр.
Затем меняю частоту основной гармоники (скажем 47.35 ставлю), не меняя частоты дискретизации. Смотрю на спектр - искажён.... вот тут, получается, и надо как-то частоту дискретизации поменять, и, соответсвенно, временное окно анализа.

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

Так вот главный вопрос всё ещё актуален: фильтр фарроу меняет длину временного окна? как я понял он только выборки сдвигает во времени и всё((

Так он может сдвигать на разные величины сдвига. Поэтому можно получить любое количество промежуточных отсчетов в любых промежуточных точках в каких пожелаете. В частности, можно провести ресамплинг на любую частоту дискретизации c*Fd.

К-во отсчетов меняет. А продолжительность в секундах - нет. Это уже сами добавляете или отбрасываете.

И полифазные фильтры позволяют перейти на любую "рациональную" частоту (m/n)*Fd. Полифазные фильтры, как говорилось, практически не вносят искажений в АЧХ, в отличие от. Но Вам нужно было бы поднимать частоту дискретизации на очень высокий уровень, чтобы иметь возможность отслеживать основной тон с достаточной точностью.

Фильтр фарроу количество отсчетов не меняет. Он только сдвигает сигнал, а частота дискретизации остается та же. Меняет количество отсчетов ресэмплер, если он сделан на основе фильтра фарроу.

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

Вот есть отсчеты x(0) x(1) x(2) x(3) ...
Мы хотим увеличить частоту дискретизации в 1.1 раза
Получим отсчеты x(0) x(1.1) x(2.2) x(3.3) ...

Отсчет x(1.1) находится между x(1) и x(2)
Интерполятор 3-й степени вычисляет его так:
h(0)*x(0)+h(1)*x(1)+h(2)*x(2)+h(3)*x(3)

Для отсчета x(2.2)
h(0)*x(1)+h(1)*x(2)+h(2)*x(3)+h(3)*x(4)
итд

К-ты h(0...3) вычисляются для каждого нового отсчета исходя из дробной части его индекса.

Сам по себе фильтр фильтр фарроу это просто элемент дробной задержки, у которого есть параметр mu, который определяет, сколько эта задержка составляет относительно интвервала дискретизации входного сигнала. А на основе фильтра фарроу можно построить переменную задержку или ресэмплер. Чтобы сделатьресэмплер недостаточно просто фильтра, нужна еще хитрая схема синхронизации.

Это уже философия. Так Вы человека совсем запутаете.
Фильтр - не меняет. Постоянное mu.
А интерполятор - меняет. Хоть на каждом такте mu можно менять.
Самый обычный интерполятор B-сплайном.
Хоть для сдвига (фильтра), хоть ресамплинга.
Короче, если человек хочет поменять, то он может поменять с помощью интерполятора фарроу.
Трудно даже задачу представить, где mu не меняется. Скорее фильтр - это когда mu меняется медленно, а
интерполятор, когда mu меняется быстро

Хотя может если в FPGA то и хитрая синхронизация...

Прошу прощения Я в матлабе не работаю. Всегда на живом процессоре. Сейчас LPC2478.
Я думаю разработка или публикация по практической тематике должны производиться в следующем порядке:
1. При помощи своего или чужого прибора получить запись интересующего процесса.
2. На живом процессоре или в матлабе(по желанию) получить желаемое.
3. Проверить на других объектах.
4. Приступить к созданию прибора или публикации.
В Вашем случае п. 1 не выполнен.
А кто Вам сказал, что эти синусы там имеются и какой величины и что 12бит хватит?
У WU Jiekang 16бит. У меня 24 и меньше не получится.
А что касается Фарроу, то это способ пересчитать сэмплы в нужные моменты времени.(Ресэмплинг)
Просто берете выборку с запасом и пересчитываете.
Еще не ясно как реальный сигнал на это издевательство прореагирует.
Если выполните п.1 помогу с п.2.

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

Я же хочу оценить погрешность сверху... где я возму сеть с амплитудой 40-й гармоники равной 1-й... она по-идее вообще копеечная))
Предлагаете в розетку врубиться и оцифровать сигнал???
Так это получится частный случай... тем более из домашней электросети...

Для постоянной дробной задержки я предпочту использовать более качественный (длинный) фильтр НЧ.
Полифазный фильтр с одной фазой.
Для этого случае фарроу не нужен, поскольку это плохой фильтр НЧ с плохой АЧХ.
Интерполяторы придумали как раз из-за того, что меняется дробная mu.
Разве, что у Вас предельно высокие частоты.


В данной конкретной задаче автору хочется провести ресамплинг входного сигнала в пропорции основной частоты сети к 50 гц.
Может это не лучший, не единственный метод, но так ему хочется. Чтобы гармоники легли на бины DFT
Это небольшое изменение частоты, но не сдвиг. Для 45 гц отсчёты нужно пересчитать на сетку с шагом 1.1 (mu=[x+0.1*k]).
Количество отсчётов изменится в этой пропорции. Если было 8192, станет 7372. 820 отсчёта нужно будет добавить нулями или лучше с самого начала захватывать отсчеты с запасом, 9011. В противоположном случае 55 гц, шаг 0.9, 8192 отсчёта превращается в 9011. 820 отсчетов нужно отбросить если использовать FFT.
Так что число отсчетов меняется.



Я тоже так думаю. Этот спектр оценивается при произвольном сэмплинге, без привязки к бинам DFT.
Если только оценивать нужно только параметры кратных гармоник, а всё что между ними неважно.
Как это делают в Стэнфорде. Про идею я уже говорил раньше: квадратичная интерполяция спектра.
Гауссово окно - оно гауссово и в спектре и после логарифмирования точки лягут вблизи максимумов на параболы очень точно. Ширину гауссова окна подобрать так, чтобы гармоники друг-на-друга не влияли. Нулей добавить побольше, раз в 8 больше чем данных, чтобы минимизировать систематическую ошибку определения максимумов.
Вот эта работа, она общедоступна. Там ещё отличная библиография
https://ccrma.stanford.edu/STANM/stanms/sta...14/stanm114.pdf
https://ccrma.stanford.edu/STANM/stanms/sta...15/stanm115.pdf
https://ccrma.stanford.edu/STANM/stanms/sta...16/stanm116.pdf

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



А с этим я даже и не пытался спорить Ресэмплер на основе фарроу естественно изменит объем выборки в количестве отсчетов.

Можно применить сдвиг частоты так, как это делается в радиотехнике для OFDM. Корреляторы бинов ДФТ сдвигают свою опорную частоту. Так проще , чем интерполировать входной сигнал.