минимизация погрешности, при плавающей частоте Опции

[TigerSHARC]

Есть сигнал. Частота основной гармоники 50 гц.
Используется алгоритм ДПФ. При смещении частоты на 1-5 Гц. Появляется дополнительная погрешность по модулю ДПФ. Оконное сглаживание не годится.
частота дискретизации 1200гц задана жёстко и не меняется, выборок строго 24 за период.
Слышал про алгоритм интерполяции, но так ничего конкретного не нашёл.
Может кто сталкивался или известен источник.????
Необходимо минимизировать погрешность БПФ при смещении частоты сигнала при постоянной частоте дискретизации.

Сообщение отредактировал TigerSHARC - Sep 4 2009, 10:49

[Den]

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