[Den 0]

Спасибо.

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

 

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

[Serg76 0]

Уважаемые знатоки!

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

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

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

описываются и т.д., да и вообще есть ли какие-нибудь источники по этой теме. спасибо.

[fontp 0]

Уважаемые знатоки!

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

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

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

описываются и т.д., да и вообще есть ли какие-нибудь источники по этой теме. спасибо.

 

Для одиночной комплексной синусоиды 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_sinusoids_CRLB_best_worst.pdf

[Serg76 0]

2 fontp

 

Спасибо большое за такой развернутый ответ! интуитивно я так себе все это и представлял.

действительно, наверное с практической точки зрения проще не заморачиваться и использовать

однотональные алгоритмы, но при этом всегда надо помнить об "удаленности" паразитных гармоник

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

на случай полигармонического воздействия, во всех алгоритмах точность упала где-то на порядок в

рабочем диапазоне С/Ш при использовании прямоугольного окна. при других окнах точность еще

хуже, в отдельных алгоритмах уменьшение наблюдалось еще на два порядка. Но Маклеод лучший :),

как и рекомендуют буквари.

 

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

бигармоническом воздействии.

 

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

[Serg76 0]

2 fontp

 

пока тема не "остыла", хотел еще спросить: помнится в одной из веток (не могу найти в какой :() Вы рекомендовали

при оценке частотной отстройки M-PSK сигналов использовать алгоритм Viterbi & Viterbi вместо возведения в степень.

Не могли бы Вы по-подробнее рассказать суть алгоритма и в чем его достоинство? Спасибо.

[fontp 0]

2 fontp

 

пока тема не "остыла", хотел еще спросить: помнится в одной из веток (не могу найти в какой :() Вы рекомендовали

при оценке частотной отстройки M-PSK сигналов использовать алгоритм Viterbi & Viterbi вместо возведения в степень.

Не могли бы Вы по-подробнее рассказать суть алгоритма и в чем его достоинство? Спасибо.

 

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

Я приводил эти.

 

http://electronix.ru/forum/index.php?showt...mp;#entry720790

Преимущество в том, что работает при более низком сигнал/шуме, но конечно слепая подстройка частоты конечно

очень сильно уступает подстройке по данным. Нужен очень длинный "пакет" данных, да и то ниже некоторого порога

сигнал/шум уже не накопляется даже увеличением длины пакета

 

2 fontp

 

Спасибо большое за такой развернутый ответ! интуитивно я так себе все это и представлял.

действительно, наверное с практической точки зрения проще не заморачиваться и использовать

однотональные алгоритмы, но при этом всегда надо помнить об "удаленности" паразитных гармоник

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

на случай полигармонического воздействия, во всех алгоритмах точность упала где-то на порядок в

рабочем диапазоне С/Ш при использовании прямоугольного окна. при других окнах точность еще

хуже, в отдельных алгоритмах уменьшение наблюдалось еще на два порядка. Но Маклеод лучший :),

как и рекомендуют буквари.

 

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

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

Но она есть, и она будет определяющей при высоком отношении сигнал шум > 40 дб

 

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

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

 

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

Десятка два параметрических методов описано в книге Марпла по спектральному анализу.

Кухня с интерполяцией спектра ДПФ в основном нашла применение для низкого сигнал-шума, в основном для модемов

[Serg76 0]

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

Я приводил эти.

 

http://electronix.ru/forum/index.php?showt...mp;#entry720790

Преимущество в том, что работает при более низком сигнал/шуме, но конечно слепая подстройка частоты конечно

очень сильно уступает подстройке по данным. Нужен очень длинный "пакет" данных, да и то ниже некоторого порога

сигнал/шум уже не накопляется даже увеличением длины пакета

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

 

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

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

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

 

Но она есть, и она будет определяющей при высоком отношении сигнал шум > 40 дб

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

 

Кухня с интерполяцией спектра ДПФ в основном нашла применение для низкого сигнал-шума, в основном для модемов

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

 

[Den 0]

Добрый день, уважаемый 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.

 

[fontp 0]

Добрый день, уважаемый 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

 

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

[Den 0]

 

Добрый день.

 

Пытаюсь всё разобраться. Вот ещё код.

Имею сигнал 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.