... не совсем понимаю, что даёт интерполяция сама по себе ...
С помощью интерполяции (добавление нулей и фильтрация) мы получаем дополнительные точки для решения задачи во временной области. Очень важна фильтрация, которая позволяет убрать весь шум выше искомой частоты и точно расставить интерполяционные точки. Как вариант, фильтр должен делать свертку фрагмента выборки с добавленными нулями с импульсной характеристикой фильтра на каждый сдвиг, то есть с увеличенной Fs. Размер фильтра, как предлагал bahurin, может быть больше выборки. Допустим, коэффициент интерполяции равен 10, тогда точек в выборке будет 24х10=240. Размер фильтра может быть и больше, что желательно для лучшего подавления частот выше первой гармоники (то есть, можно подавлять со 100 Гц). Кстати, фильтр лучше сделать и до интерполяции, так как 24 точек хватит, тем более, что импульсную характеристику лучше взять длиннее. Еще раз, интерполяция нужна для решения задачи во временной области.
Насчет FFT на каждые 24 точки - это ведь не означает, что поток точек прерывистый. Применение FFT на 24 или на большее число точек нужно делать скользящим, т.е каждый раз со сдвигом на каждую точку (в 24 раза больше). После этого можно провести дополнительную обработку с целью получить пользу из этой избыточности.
TigerSHARC, Вы бы дали выборку для анализа, было бы легче. Только размер нужен, естественно, хотя бы в несколько периодов подряд.

... зачем нулями-то? Лучше же приближённо реальные точки найти... линейно, или сплайнами...
смоделирую в Матлаба завтра всё то, что вы сказали...
а насчёт выборки... ну вот примерно так в матлаб и задавать нужно

f = 50+(rand(1,1)-0.5*10); %случайная частота 50+-5 Гц
s0 = 1.4*sin(2*pi*f*t).*(1+0.5*cos(0.5*pi*t+rand(1,1)*2*pi));% сигнал с небольшим гулянием амплитуды

Как раз только НУЛЯМИ - это и есть настоящая интерполяция. Так мы не портим спектр, а только клонироем его на высоких частотах, а потом все эти компоненты ликвидируем фильтрацией (естественно, с приемлемой точностью, определяющейся затуханием интерполяционного фильтра). В результате мы увидим, что число точек во временной области увеличилось в N раз и они расставились после фильтрации именно наилучшим образом. В этом вся суть интерполяции. Ничего не надо делать приближенно, по крайней мере, вначале для получения максимального результата. Только добавление нулей и фильтрация. В Матлабе при просмотре структуры интерполяционного фильтра можно увидеть эти нули.

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

Вы только что сами все сказали интерполяция просто клонирует спектр. И если по исходному спектру нельзя точно померить, то и по "клонированному" тоже не получится, так как все повторяющиеся куски полностью совпадают с исходным. После фильтрации в идеале опять останется исходный спектр, по которому нельзя точно померить. А поскольку идеального фильтра нет, то вы только уменьшите точность оценки амплитуды. Еще раз повторяю: интерполяция не вносит никакой дополнительной информации, поэтому не может улучшить точности оценки.

Что касается фильтрации до интерполяции, то опять таки отсутствие идеального фильтра не полностью подавит высшие гармоники и использовать геометрические "точные" методы также не получиться. Этот вопрос уже обсуждался на форуме экспоненты.

Так я ж не спорю, что интерполяция не генерирует дополнительную информацию - ей неоткуда взяться на том же числе точек. Интерполяция позволяет с максимальной точностью (используя предоставленную выборку) решить задачу на временной области (не на частотной, а на временной. И только). А предварительная фильтрация все равно улучшит сигнал (насколько сможет). Допустим, гармоника имеют -40 дБ (что хорошо), а фильтр -60. Наверняка хватит и отличаться от максимально возможного почти не будет.
На частотной области решения должны быть с той же точностью при условии использования того же временного окна. Для увеличения разрешения по частоте нужно выборку дополнить нулями и взять больший размер FFT (это тоже обсуждалось).
Если сравнивать FFT или свертку, то разницы конечного результата не будет при всех равных условиях.

просто применил окно к сигналу, и получил погрешность чуть выше 1%

Код MATLAB:

T = 0.02;
Ts = 1/1200;
f = 55;
t = 0:Ts:T-Ts;
A = 10;
df = 1/T;
Fmax = 1/Ts;
n = length(Y);
w = hamming ( n, 'periodic');
f1 = 0:df:Fmax-df;
Y = A*sin(2*pi*f*t) ;
Xw = fft(Y.*w.')/sum(w);
X = fft(Y)/n;
subplot(2,1,1), stem(f1, abs(X)), grid
title('Спектр сигнала 24 выборки без окна')
subplot(2,1,2), stem (f1, abs(Xw)), grid
title('Спектр сигнала 24 выборки c окном')

... как ещё можно минимизировать погрешность???

Ещё можно адаптивной фильтрацией.
Вот здесь минимизируют через LMS предсказание синуса вместе с тремя гармониками
http://rapidshare.com/files/299646763/Nove...acking.pdf.html
Но следящая система будет работать не спешно



Так это без шума. Без шума интерполяция в спектральной области даст почти бесконечную точность.
Сделайте преобразование Гильберта. Добавить в сигнал в конце нулей (от n до 3*n штук ) по вкусу, сделайте fft, найдите максимум, по 3-м точкам в окрестности максимума проведите параболу, найдите максимум этой параболы, аргумент максимума - это и есть частота

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

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

На входе почти синусоида, сиг/шум~4, диапазон 10..90 Гц, погрешность оценки частоты +/-0.02 Гц по спектру за 200 точек (частота выборок 200 Гц, время измерения 1 с). Все утопталось в ATmega8, осталось под индикатор и обмен с PC. Алгоритм используется полгода в другом серийном устройстве.

Это же противоречит принципу неопределённости!!!
Щас придут умные дяди и по голове настучат

TigerSHARC, Вы определитесь - Вам вычисление частоты на каждом периоде нужно или, всё-таки, пересчёт действующего значения основной гармоники?
Подозреваю, что нужно именно действующее значение основной гармоники. Вот его и считайте по 24 точкам.
Что касается частоты. Поскольку речь идёт о промышленной частоте, то надо понимать, что энергосистема - штука ну очччень инертная и за один период частота не уйдёт не то что с 50 на 49 (или 51), но и даже в пределах требуемой по стандарту точности определения, а значит такого разрешения по времени и не требуется. Считайте частоту, как минимум, на 10 периодах для защит и быстрой автоматики и (см. выше у AlexU) - 0,5 - 1,0 секунду для индикации.
Все видели кадры с Саяно-Шушенской? Железная "юла" там немаленькая, правда? А уж какая тяжёлая... Так вот, в штатном режиме такая турбина каждые 20 мс свою скорость не менять не будет... а вот динамические искажения сигнала и броски фазы быть могут (локальная нагрузка коммутируется и меняет местный режим сети довольно шустро, но не частоту!) - их то и будете отлавливать весьма успешно, но для дела бестолково.
Кстати, какова разрядность АЦП и какая точность определения частоты требуется? 0.01 Гц?

На работе прибора это никак не отразится.
Пока разговор вертится вокруг ограничений FFT и цитатами из учебников будем спокойно его выпускать.

"Умные дяди" как раз утверждают, что это предельная точность, соответствующая критерию Крамера-Рао
Стандартное отклонение = примерно Fd/(SNR*NS**(3/2))
Можете сами проверить, что он именно её получил.

При измерении частоты во всём представимом диапазоне частот (от 0 до Fd/2) точность оценки улучшить нельзя,
поскольку критерий Крамера-Рао - это оценка максимального правдоподобия.
Если диапазон частот узкий, наверное, точность можно улучшить предварительной фильтрацией и передискретизацией. Только предварительная фильтрация вносит ещё и задержку, что противоречит вообще-то требованию, что измерение за минимальное время.

А о разрешении (принципе неопределённости) можно говорить только когда необходимо различать хотя бы 2 спектральные линии. Не йорничайте.
Точность измерения частоты единственной линии достигается, как говорят, априорной информацией, что она там одна. Но это уже метафизика.

Все верно, +1
Кстати, на практике увеличение частоты выборки не приводит к заметному увеличению точности. (Хотя диапазон соответственно расширяется)

Откуда сигнал и какой он природы?

Датчик (чего-то там ) частотный. На выходе от 7 до 3000 Гц. От 7 до 100 Гц сигнал очень слабый много шума s/n до 1,5. Бывают кратковременные импульсные помехи большой амплитуды. Форма сигнала зависит от частоты (5мВ синусоида <40Гц)-(15мВ синусоида ограниченная сверху или снизу, плавает постоянная составляющая <120Гц) -(20..500мВ синусоида<500Гц)-(3В почти меандр). Измерять по спектру - самое то. Хотя я бы для начала сам датчик переделал.

мне очень интересно узнать вот такую вещ:
Есть сигнал (50Гц + гармоники). Вычисляем с помощью БПФ значения гармоник кратных основной частоте(50Гц).
Частота дискретизации фиксирована, временное окно наблюдения тоже. Проблема в том, что как только основная частота не равно 50Гц, то значения гармоник искажаются (оно и понятно, т.к. в окне наблюдения умещается не полное число периодов).
Основная задача - это узнать максимально достоверную информацию о значении гармоник.
При этом приоритетными задачами являются: постоянная(аппаратная) частота дискретизации, и минимальное временное окно наблюдение (насколько это возможно).

Если аппаратно, то речь, как я понимаю идёт о ПЛИС...
Тогда я Вам советую реализовать цифровой узкополосый КИХ-фильтр фильтр для 50 Гц (Матлаб и fdatool в помощь), который порежет всё остальное. К очищенному сигналу применить простой, но достаточно эффективный метод определения периода по переходам через ноль. И от этого "плясать" дальше.
Естественно, учитывать надо будет не 1 период, а 10-20 и момент перехода через ноль интерполировать (для дискретизации 1200 - линейно). Будет полезно сверять однородность периодов (хотя бы примерно) внутри интервала, чтобы избежать погрешностей от динамических искажений фазы.

предполагалось использовать DSP-процессор.
А если требуется узнать информацию о большем количестве гармоник(20-40 гамоник), что бы на их основе коэффицент несинусоидальности посчитать...

Устройство контроля качества электроэнергии?
Вообще, на это дело есть пара иностранных стандартов и наш ГОСТ, где ширина окна наблюдения и интервалы усреднения оговариваются. Ищите...

Ну тогда 1200 Гц это не просто мало, а мало до нереализуемости. Приглядитесь к частотам дискретизации 6кГц-10кГц

)) понятно что мало... пусть даже для предотвращения наложения используется RC-цепь и передискретизация...
Как минимизировать погрешность для всех гармоник (пусть 20 гармоник) при том что основная частота сигнала 45-55Гц... (т.е. нужно как то предотвратить растекание спектра...

8-10-ая гармоника это максимум при 1200 Гц.

Да ладно Вам, см. ГОСТ 13109-97

Если частота вышла за эти пределы, то коэффициент несинусоидальности как-то малоинтересен, т.к. в системе начинают происходить совсем другие процессы.

Да знаю я теорему котельникова!!!

ладно объясню по порядку:
При ДПФ есть два источника методической погрешности:
1) наложение частот ввиду несоблюдения теоремы отсчётов;
2) неполное число периодов во временном окне;

ТАК ВОТ МЕНЯ ИНТЕРЕСУЕТ 2 ПУНКТ!!!!
пусть условия теоремы котельникова соблюдены абсолютно(эта задача уже решена)

Я говорю, то что знаю точно!!!
Существуют промышленные приборы(видел лично).
В руководстве по эксплуатации сказано, что выдаются значения до 40-й гармоники с такой-то погрешностью при том, что частота сети 45-55Гц.

Все эти приборы работают с окном наблюдения намного превышающем период сети, а потом результаты от нескольких замеров ещё и усредняются. Значение частоты, кстати, усредняется за время 20 секунд(!).
Если на пальцах, то выбрал один из двух способов определения частоты:
а) по переходам через ноль (в тяжёлых случаях - сигнал чистить)
б) из самого же Фурье - по приращению фазы основной гармоники. Грубо говоря, преобразовали на тех же 20 мс, получили модуль1-фазу1, на следующих 20 мс опять модуль2-фазу2 и т.п. Если частота сети идеальна, то разница фаза2 - фаза1 будет равна нулю. Если же частота меньше 50Гц или больше, то появится эффект скольжения частоты и Вы будете наблюдать поворот фазы в "+" или в "-". При 45(55)Гц это скольжение будет видно на двух соседних интервалах 20мс, при отклонении частоты на 0.05Гц период наблюдения увеличивается.
Допустим, с частотой разобрались. Начинаем изменять период таймера-задатчика частоты дискретизации в соответствии с определённой нами частотой, натягивая сетку квантования времени на реальный сигнал. При этом продолжаем "измерять" частоту.
В итоге добиваемся необходимой кратности частоты дискретизации и сетевой и все прочие расчёты остаются нетронутыми - никаких пересчётов коэффициентов, переменной длины буфера, комплексного доворота и т.п.

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

А чем Вас так держит аппаратная частота?
Можно менять её и не аппаратно а только для расчётов. А исходный сигнал, если необходимо сохранять осциллограммы, записывать на железной частоте.

Спасибо за совет! Думаю это то что надо!

Правильное решение, но менять собственную частоту дискретизации - это лишнее. В GPS приемниках частоту не меняют, но слежение за фазой позволяет измерять задержки до спутника с точностью до наносекунд.

Правильное решение - именно интерполировать сигнал относительно расчитанного значения основной гармоники.

Добрый день!
Прошу прощение за поднятие давней темы.
Где сейчас возможно достать файлы, обозначенные ссылкой?
Буду очень признателен, если кто-то сможет выслать, или где то выложит.

С уважением, Den.

"достать" можно здесь
http://electronix.ru/forum/index.php?showt...mp;#entry726504

Спасибо.

Уважаемый fontp, прочитал много Ваших постов. Позвольте задать Вам вопрос, как человеку, разбирающемуся в этой теме.
Правильный ли ход моих рассуждений по поводу следующей задачи…

Сетевое напряжение 50 Гц. Необходимо вычислять фазу, частоту, амплитуду. Также необходимо проводить анализ гармонического состава.
Значение фазы, частоты и амплитуды необходимо получать как можно быстрее…

Планирую так:
АЦП 16 бит, частота дискретизации 6400Гц, FIFO буфер БПФ на 512 точек, т.е. получается 128 точек на период.
При поступлении каждой новой точки с АЦП в FIFO буфер, на него накладывается окно (Гауссовское?), производится БПФ. Далее нахожу амплитуду и фазу основной гармоники.
Остается частота основной гармоники, необходимая точность измерения которой равна 0,01Гц. Возможно ли её определение по скорости изменения мгновенной фазы основной гармоники? Или же необходимо построение параболы по максимумам?

С уважением, Den.

получается вы хотите считать 128 точечное БПФ за 156,25 мкс. успеете?

128 точечное БПФ у меня высчитывается за 1.1 мкс.
БПФ хочу делать по 4-м периодам, соответственно 512 точечное.


С уважением, Den.

Амплитуда и фаза максимального бина будут скакать из-за размывания спектра.
Наоборот, сначала вы строите параболу по максимальному и соседним бинам. Определяете положение максимума (и возможно амплитуду).Потом уже гармонику точной частоты накладываете на входной сигнал и точно определяете фазу и амплитуду.
Известно, что такой способ определения частоты и фазы достигает предельно возможной точности (максимального правдоподобия), при заданной длительности сигнала и умереннных отношениях сигнал шум (0-30 дб). Большинство других методов, тот же ФАПЧ, не оптимальны с точки зрения измерений по коротким последовательностям. То есть точность может быть достигнута, но время захвата не соответствует критериям Крамера-Рао

Я приводил ссылки на серию исчерпывающих исследований Стэнфордского университета по поводу частного способа интерполяции спектра
с добавлением нулей
http://electronix.ru/forum/index.php?showt...amp;hl=Stanford

Там есть всё - выбор окон, влияние нестабильности частоты сигнала, и т.д. и т.п.

А разве наложение весового окна не устранит (не минимизирует) растекание спектра?

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

При интерполяции спектра окна используются совсем с другой целью - изолировать область максимума спектра от других удаленных спектральных пиков, завалить им хвосты. Причем доказано, что если априорно известно, что других пиков в спектре нет - то никаких окон не нужно, прямоугольное окно оптимально. А верхушка любого окна, обычно так или иначе в области максимума аппроксимируется параболой по трем точкам. Если точки взять поближе к спектральному максимуму. Поскольку для трех точек на параболе есть простое аналитическое решение, а парабола верхушку приближает обычно не намного хуже, чем известная форма окна по многим точкам (большинство этих многих точек всё равно тонут в шуме).

Большое спасибо, буду разбираться.

С уважением, Den.



Ещё один вопрос.

Под наложением гармоники точной частоты, Вы имеете ввиду свертку с комплексной синусоидой, частоты определенной по максимуму параболы?

ну да. Только это называется не свертка, а корреляция. Она будет определять соответствующую амплитуду и фазу. Эта не будет скакать, поскольку в ней исключено влияние дискретности ДПФ - она соответствует физике

Спасибо.
Буду пробывать.

С уважением, Den.

Уважаемые знатоки!
Для обсуждаемых здесь алгоритмов оценки частотной отстройки в случае моногармонических сигналов в принципе понятно.
В продолжение темы хотелось бы узнать: какие существуют на данный момент практические алгоритмы в случае воздействия
более сложных сигналов, например, полигармонических. какие в этом случае существуют потенциальные оценки точности, какими выражениями
описываются и т.д., да и вообще есть ли какие-нибудь источники по этой теме. спасибо.

Для одиночной комплексной синусоиды Rife и Burstin первыми дали статистическую оценку Крамера_Рао для одиночной синусоиды и
привели алгоритм, который в широких пределах параметров её достигает
Я загружал когда-то эту статью
http://electronix.ru/forum/index.php?showt...mp;#entry693650

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

Что значит "достаточно далеко"? Интуитивно можно догадаться, что это много дальше чем разрешение по критерию Релея 1/N
Насколько много не очень известно и зависит от многих деталей реализации.

Также существует теоретический результат, о том, что ситуацию нельзя в принципе улучшить, когда пики находятся достаточно близко.
Нижняя граница Крамера-Рао существует для любого многомерного случайного распределения. Это статистический критерий, который может , в принципе, быть получен в явном виде для любого многомерного нормального распределения. Такая граница была получена в том числе и для двух синусоид в белом гауссовском шуме. Результат очень громоздкий и не очень полезный для практических применений. Но ассимптотики он даёт такие как надо интуиции. Когда спектральные пики находятся далеко результат не отличается от критерия Крамера-Рао для одной синусоиды.
Когда пики приближаются к друг другу на расстояние соизмеримое с "разрешением" 1/N граница резко падает от
корень из(6/(N*(N-1)*(N-1)*(Es/No)))
к 1/N и хуже

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

непрактичный пруф для 2-х синусоид

2 fontp

Спасибо большое за такой развернутый ответ! интуитивно я так себе все это и представлял.
действительно, наверное с практической точки зрения проще не заморачиваться и использовать
однотональные алгоритмы, но при этом всегда надо помнить об "удаленности" паразитных гармоник
от основной. точность при этом, конечно падает, пробовал гонять разные однотональные алгоритмы
на случай полигармонического воздействия, во всех алгоритмах точность упала где-то на порядок в
рабочем диапазоне С/Ш при использовании прямоугольного окна. при других окнах точность еще
хуже, в отдельных алгоритмах уменьшение наблюдалось еще на два порядка. Но Маклеод лучший ,
как и рекомендуют буквари.

кстати, в рекомендуемой Вами работе Маклеода рассматривается работа "оценщиков" при
бигармоническом воздействии.

Еще раз большое спасибо.

2 fontp

пока тема не "остыла", хотел еще спросить: помнится в одной из веток (не могу найти в какой ) Вы рекомендовали
при оценке частотной отстройки M-PSK сигналов использовать алгоритм Viterbi & Viterbi вместо возведения в степень.
Не могли бы Вы по-подробнее рассказать суть алгоритма и в чем его достоинство? Спасибо.

Смысл в том что фаза учетверяется, но амплитуда не возводится в степень. Вы гуглом можете найти много статей V&V
Я приводил эти.

http://electronix.ru/forum/index.php?showt...mp;#entry720790
Преимущество в том, что работает при более низком сигнал/шуме, но конечно слепая подстройка частоты конечно
очень сильно уступает подстройке по данным. Нужен очень длинный "пакет" данных, да и то ниже некоторого порога
сигнал/шум уже не накопляется даже увеличением длины пакета



Лучший из простых оценщиков непосредственно по ДПФ. У любого алгоритма есть своя область применимости.
Маклеод устойчив к сильным шумам и у него из простых самая меньшая систематическая ошибка.
Но она есть, и она будет определяющей при высоком отношении сигнал шум > 40 дб

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

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

ок, спасибо. почитаю


Да, наблюдал это при моделировании.


к сожалению, о таком канале можно только мечтать


мне это и надо Еще раз спасибо.

Добрый день, уважаемый fontp!

Что-то я совсем запутался.
Накидал код в Matlab.


clear ;
%--------------------------------------------------------------------------
A = 50 ; % Amplitude 1-harmonics 50V , Идеальная амплитуда 1-й гармоники ;
Fmain = 50 ; % F 1 harmonics, Идеальная частота 1-й гармоники ;
%--------------------------------------------------------------------------
SNR = 80 ;
%--------------------------------------------------------------------------
Fs = 6400 ; % F sampling , частота дискретизации
N = 512; % количество отсчетов
%--------------------------------------------------------------------------
t = (0:N-1)/Fs; % вектор времени
freq = 0:Fs/N:200 ; % вектор частоты
%--------------------------------------------------------------------------
xx = A*sin(2*pi*Fmain*t) ; %
x = awgn(xx, SNR, 'measured');
%x = x + randn(size(x)) ;
%--------------------------------------------------------------------------
NFFT = abs(fft(x, N))./(N/2);
%--------------------------------------------------------------------------



subplot(2,1,1), plot(t,x), axis([0 t(N) -(1.4*A) (1.4*A)]) ;
hold on ;

subplot(2,1,2), stem(freq, NFFT(1:17)), grid on ;
hold on ;


subplot(2,1,1);
text(0, -20, sprintf('Идеальная частота %2.3f Гц\n',Fmain));

text(0, -50, sprintf('Идеальная амплитуда %2.3f В\n',A));



for k = 0.8:0.01:1.2
xx = A*sin(2*pi*(Fmain*k)*t) ;
x = awgn(xx, SNR, 'measured');
NFFT = abs(fft(x, N))./(N/2);

[Ampl_max N_bin_max] = max(NFFT) ; % максимальный бин и его значение

y_l = NFFT(N_bin_max - 1) ; % амплитуда левого бина
y_c = NFFT(N_bin_max) ; % амплитуда максимального центрального бина
y_r = NFFT(N_bin_max + 1) ; % амплитуда правого бина

x_l = (N_bin_max - 2)*(Fs/N) ; % частотное положение бинов
x_c = (N_bin_max - 1)*(Fs/N) ;
x_r = N_bin_max *(Fs/N) ;


inp_vector_x = [x_l x_c x_r] ;
inp_vector_y = [y_l y_c y_r] ;

out_vector_x = x_l:0.01:x_r ;

out_vector_y = lagrange(inp_vector_x, inp_vector_y, out_vector_x);


[ymax xmax] = max(out_vector_y) ;

A_calc = ymax ;

Freq_izmer = xmax ;



errorA = (abs(A - A_calc))*100/A ;


subplot(2,1,1) ;
h(1) = plot(t, x, 'r') ;
h(3) = text(0, -35, sprintf('Измереная частота %2.3f Гц\n', Freq_izmer));


h(4) = text(0, -65, sprintf('Измереная амплитуда %2.3f В\n', A_calc));
h(5) = text(0.009, -65, sprintf('( %2.3f )\n', errorA));






subplot(2,1,2);
h(2) = stem(freq, NFFT(1:17), 'r') ;
h(6) = plot(out_vector_x, out_vector_y, 'r') ;



drawnow ;
pause(0.5) ;
if k == 1.19
return
end
delete(h) ;
end




function yy=lagrange(x,y,xx)
% вычисление интерполяционного полинома в форме Лагранжа
% x - массив координат узлов
% y - массив значений интерполируемой функции
% xx - массив значений аргумента, для которых надо вычислить значения полинома
% yy - массив значений полинома в точках xx

% узнаем число узлов интерполяции (N=n+1)
N=length(x);
% создаем нулевой массив значений интерполяционного полинома
yy=zeros(size(xx));
% в цикле считаем сумму по узлам
for k=1:N
% вычисляем произведения, т.е. функции Psi_k
t=ones(size(xx));
for j=[1:k-1, k+1:N]
t=t.*(xx-x(j))/(x(k)-x(j));
end
% накапливаем сумму
yy = yy + y(k)*t;
end



Получается, что построением параболы, амплитуду точно определить не получается?
Или я совсем уже туплю…
Ткните, пожалуйста, носом, в нужную сторону, а лучше конечно в пример Matlab.

С уважением, Den.

Так Вы ничего и не считаете... Зачем Вам амплитуда центрального бина?
Но мы говорили о комплексной экспоненте. А не действительном синусе. От действительного синуса еще нужно как-то получить квадратуру. Такие статьи тоже приводились. Там
http://electronix.ru/forum/index.php?showt...66966&st=30
Для комплексной экспоненты берете любой метод из tst3.zip, лучше маклеод. Получаете интерполированое значение частоты.
Умножаете свой сигнал, на комплексную экспоненту полученой частоты комплексно сопряженную и суммируете для усреднения.
Или добавляете нулей , делаете дПФ, проводите параболу через три точки вокруг центрального бина. Берете не центральный бин, а верхушку параболы. Как здесь.
https://ccrma.stanford.edu/STANM/stanms/sta...14/stanm114.pdf

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

Добрый день.

Пытаюсь всё разобраться. Вот ещё код.
Имею сигнал cod_ADC.


clear ;
%--------------------------------------------------------------------------
A = 50 ; % Amplitude 1-harmonics 50V , Идеальная амплитуда 1-й гармоники ;
Fmain = 51; % F 1 harmonics, Идеальная частота 1-й гармоники ;
%--------------------------------------------------------------------------
SNR = 90 ;
%--------------------------------------------------------------------------
Fs = 6400 ; % F sampling , частота дискретизации
N = 256 ; % количество отсчетов
%--------------------------------------------------------------------------
t = (0:N-1)/Fs; % вектор времени
freq = 0:Fs/N:200 ; % вектор частоты
[ee, N_points_freq] = size(freq) ;
%--------------------------------------------------------------------------

x_exp = A*exp(j*(2*pi*Fmain*t));
cod_ADC = A*cos(2*pi*Fmain*t) ;


%cod_ADC_complex = hilbert(cod_ADC) ;
cod_ADC_complex = cod_ADC ;


real_cod_ADC = real(cod_ADC_complex) ;
imag_cod_ADC = imag(cod_ADC_complex) ;


real_x_exp = real(x_exp) ;
imag_x_exp = imag(x_exp) ;


%--------------------------------------------------------------------------
subplot(3,1,1), plot(t,real_cod_ADC, 'r'), axis([0 t(N) -(1.4*A) (1.4*A)]), grid on ;
hold on ;
plot(t,imag_cod_ADC) ;
%--------------------------------------------------------------------------
subplot(3,1,2), plot(t,real_x_exp, 'r'), axis([0 t(N) -(1.4*A) (1.4*A)]), grid on ;
hold on ;
plot(t,imag_x_exp) ;
%--------------------------------------------------------------------------
NFFT_cod_ADC = fft(cod_ADC_complex) ; % Спектр
Amplituds_cod_ADC = abs(NFFT_cod_ADC)/(N) ; % Амплитуды

NFFT_x_exp = fft(x_exp) ; % Спектр
Amplituds_x_exp = abs(NFFT_x_exp)/(N) ; % Амплитуды

%--------------------------------------------------------------------------
subplot(3,1,3), stem(freq, Amplituds_cod_ADC(1:N_points_freq), 'r'), grid on ;
hold on ;

stem(freq, Amplituds_x_exp(1:N_points_freq)), grid on ;
hold on ;
%--------------------------------------------------------------------------
[Ampl_max_cod_ADC N_bin_max_ampl_cod_ADC] = max(Amplituds_cod_ADC) ;
[Ampl_max_x_exp N_bin_max_ampl_x_exp] = max(Amplituds_x_exp) ;

%--------------------------------------------------------------------------
piak3vect_cod_ADC(1:3)=NFFT_cod_ADC(N_bin_max_ampl_cod_ADC-1:N_bin_max_ampl_cod_ADC+1); % Isolated 3 samples around peak.
piak3vect_x_exp(1:3)=NFFT_x_exp(N_bin_max_ampl_x_exp-1:N_bin_max_ampl_x_exp+1);

%--------------------------------------------------------------------------
macleod_error_cod_ADC = macleod(piak3vect_cod_ADC) ;
macleod_error_x_exp = macleod(piak3vect_x_exp) ;

%--------------------------------------------------------------------------
% quin_error = quin(piak3vect) ;
% quadterp_error = quadterp(piak3vect) ;
% quin2_error = quin2(piak3vect) ;

min_er_x = Fmain - 0.1 ;
max_er_x = Fmain + 0.1 ;

figure ;


subplot(2,1,1), stem(Fmain, 1, 'r'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;
hold on ;
stem(((N_bin_max_ampl_cod_ADC-1)*(Fs/N))+(macleod_error_cod_ADC*(Fs/N)), 1, 'g'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;

subplot(2,1,2), stem(Fmain, 1, 'r'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;
hold on ;
stem(((N_bin_max_ampl_x_exp-1)*(Fs/N))+(macleod_error_x_exp*(Fs/N)), 1, 'g'),axis([min_er_x max_er_x -1.1 1.1]), grid on ;


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

С уважением, Den.