Точно измерить частоту. Опции

[b-volkov]

Имеем на выходе датчика синусоиду примерно 2кГц в течении 0.5 сек. , т.е. цуг из примерно 1000 периодов. В зависимости от физ. величины частота меняется на +-10% Считается, что во время цуга частота не меняется. Надо определить частоту с точностью хотя бы до 0.01Гц. На данный момент метод используетсяпрямое измерение периода каждого колебания с дальнейшей статистической обработкой. При существующем уровне помех, наводок и т.д. точность получается не лучше 0.1Гц. Вопрос: как методами ЦОС уточнить результат измерения? Это вообще возможно теоретически? Сразу должен оговориться, что я о цифровой обработке имею весьма общие представления и никакой практики. Мне в голову приходит только свертка с синусоидами от (f0 – 0.1Гц) до ( f0 + 0.1Гц) с шагом 0.01Гц (f0-приблицительно измеренная частота) и поиском "резонанса". Может быть можно как по другому, попроще?

[leksa]

Всем здравствуйте!
Тема интересная, позвольте выложить матлабовский скрипт, в котором я набросал алгоритм оценки частоты комплексной экспоненты на основе планарной фильтрации.
Ничего сверхестественного, просто один из adhoc способов. Идея в том, чтобы умножая отсчеты комплексной экспоненты на свою задержанную и взятую с комплексным сопряжением копию, получить последовательность фазоров, угол которых пропорционален частоте исходной экспоненты. Затем нужно сложить полученные фазоры как комплексные числа, получив на комплексной плоскости длиный вектор, угол которого опять же пропорционален искомой частоте(эта операция и есть планарная фильтрация). Таким образом уменьшается воздействие шума на оценку. Далее вычисляем угол этого длинного вектора, умножаем на коэффициент Fs/(2*pi) и получаем оценку частоты.
Правда полученная таким образом оценка достаточно сильно подвержена влиянию шума.
Чтобы несколько улучшить оценку, я пошел немного дальше и используя вычисленную оценку частоты, создаю комплексный КИХ фильтр с ИХ в виде 1 периода комплексной экспоненты на полученной оценкой частоте. Получается такой комплексный ППФ со средней частотой совпадающей с нашей оценкой частоты.
Затем пропускаю входной сигнал через этот ППФ и снова нахожу оценку по вышеописанному методу.
Засчет второй итерации после фильтрации, погрешность оценки может уменьшиться на порядок по сравнению с первой оценкой.
Вот такой вот способ.
Правда, вторая итерация имеет смысл только в определенном диапазоне отношения сигнал шум (если ОСШ выше 45 дБ, то полезный эффект от второй итерации уже незначителен)
Планарная фильтрация используется в ряде прямых(feedforward) методов синхронизации, описанных например в
в часто упоминающейся на форуме книге Synchronization Techniques for Digital Receivers.

Собственно, скрипт:

%test_carr_freq_ff_est
script
clear
clc
close all

Fs=6400;
Tau=1.00;%
Nlen=Tau*Fs;
t=[1/Fs:1/Fs:Nlen/Fs].';
SNR=6;
Fc=47.333333333333333333333;%
disp(['SNR = ' num2str(SNR)])
disp(['Fc = ' num2str(Fc)])
ini_phase = (rand(1)*2-1)*pi
y = 1.*exp(j*(2*pi*Fc*t+ini_phase));
Sinp = awgn(y,SNR,'measured');
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d)
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(11)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est = phase_pf_pi*Fs/2;
delta_F = Fc_est-Fc;
disp(['Оценка частоты: Fc_est = ' num2str(Fc_est)])
disp(['Погрешность оценки частоты: delta_F = ' num2str(delta_F)])
figure(12)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ сигнала после фильтра');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(13)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр сигнала после фильтра'),hold on;
% пофильтруем
h_len = round(Fs/Fc_est)*1;
hbpf = exp(j*2*pi*(Fc_est/Fs)*[1:h_len]);
Sinp_f = filter(hbpf,1,Sinp);

% еще разок
Sinp = Sinp_f;
Sinp_d = Sinp(1:end-1).*conj(Sinp(2:end));
Sinp_d_pf = sum(Sinp_d);
% цикл чисто для наглядности работы алгоритма
Sinp_d_pf_k = zeros(size(Sinp_d));
for k=2:length(Sinp_d)
Sinp_d_pf_k(k)=Sinp_d_pf_k(k-1)+Sinp_d(k);
% figure(101)
% plot([Sinp_d_pf_k(k-1) Sinp_d_pf_k(k)],'y- *'),grid on,hold on
% pause(1/25)
end
figure(21)
plot(Sinp_d_pf_k,'y- *'),grid on,hold on,title('планарная фильтрация');
phase_pf = -angle(Sinp_d_pf);
phase_pf_pi = phase_pf/pi;
Fc_est2 = phase_pf_pi*Fs/2;
delta_F2 = Fc_est2-Fc;
disp(['Оценка частоты №2: Fc_est2 = ' num2str(Fc_est2)])
disp(['Погрешность оценки частоты: delta_F2 = ' num2str(delta_F2)])
figure(22)
sz=prod(size(Sinp));
tv=[1/Fs:1/Fs:sz/Fs];
plot(tv,real(Sinp),'b- .'),grid on,hold on,title('IQ входного сигнала');
plot(tv,imag(Sinp),'r- .'),grid on,hold off
figure(23)
[PSDPSD,f]=pwelch(Sinp,[],[],[],Fs,'twosided');
plot(f-Fs/2,fftshift(10*log10(PSDPSD/max(PSDPSD))),'y- '),grid on,title('Спектр вх сигнала'),hold on;
% % характеристики фильтра
% figure(31)
% plot(real(hbpf),'b- .'),grid on,hold on
% plot(imag(hbpf),'r- .'),grid on
% figure(32)
% freqz(hbpf,1)

Сообщение отредактировал leksa - Mar 3 2010, 17:28