Синхронизация Park, правильно ли я все сделал? Опции

[Ivan55]

Добрый день!

Решил я поразбираться с синхрой по преамбуле(тренировочным символам) реализовал метод Парка, не знаю правильно или нет замоделил, подскажите если что не так.
Хочу померить MSE оценок по данной преамбуле, как это можно правильно сделать?

CODE

clear all; clc; close all;

SNR=-5:10;
N_FFT=1024;
NGI=1/4;
N_used=N_FFT/2;
N_GI=round(NGI*N_FFT);
Nofdm=N_FFT+N_GI;
Nsym = 500;

nSTOs=50;
CFO=0.4;

%% Формирование тренировочных символов
B=2;
reS=(-1).^round(rand(1,round(N_used/B)));
imS=1i.*(-1).^round(rand(1,round(N_used/B)));
s1=reS+imS;
s2=fliplr(s1);
tx_signal_noGI=[s1 s2 conj(s1) conj(s2)];
tx_signal = zeros(1,Nsym*N_FFT);
tx_signal_GI=[tx_signal_noGI(1,N_FFT-N_GI+1:N_FFT) tx_signal_noGI];
tx_signal=repmat(tx_signal_GI,1,Nsym);

%% Добавление смещения частоты и времени
if nSTOs>=0, y_STO=[tx_signal(1,nSTOs+1:end) zeros(1,nSTOs)];
else
y_STO=[zeros(1,-nSTOs) tx_signal(1,1:end+nSTOs)];
end
nn=0:length(y_STO)-1; y_CFO_STO = y_STO.*exp(1i*2*pi*CFO*nn/N_FFT);

%%
buff = [];
for k = 1:length(SNR)
recvd_signal=awgn(y_CFO_STO, SNR(k), 'measured', 'dB');
for n = 1:Nsym
%%
Len_all=N_FFT;
recvd_signal_zeropad=[zeros(1,Len_all) recvd_signal(1,(n-1)*Nofdm+1:n*Nofdm) zeros(1,Len_all) zeros(1,Len_all)];
for d=1:length(recvd_signal_zeropad)-Len_all-1
P_Park(d)=sum(recvd_signal_zeropad(d+(1:N_FFT/2)).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2)))); %
E_Park(d)=sum(abs(recvd_signal_zeropad(d+(1:N_FFT))).^2);
Q_Park(d)=sum(conj(recvd_signal_zeropad(d+(1:N_FFT/2))).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2))));
end
Park_metric=(abs(P_Park).^2)./(E_Park).^2;
buff = [buff Park_metric];
plot(Park_metric/max(Park_metric))

[~, maxipos1]=max(Park_metric-N_GI);
t_est_Park(1,n)=Nofdm - maxipos1;
FreqOffset(1,n) = angle(Q_Park(maxipos1-N_GI+1))/pi;
end
mse_CFO(k) = sqrt(sum(abs(FreqOffset - CFO).^2)/Nsym)
mse_STO(k) = sqrt(sum(abs(t_est_Park - nSTOs).^2)/Nsym)
end

Сообщение отредактировал Ivan55 - Jul 14 2017, 09:34

[Ivan55]

Если быть точнее то я щас делаю вот, так

CODE

clear all; clc; close all;

flTimeShift = 1; % 1 - с уходом времени; 0 - без ухода времени
flChann = 0; % 1 - с каналом; 0 - без канала

%choose SNR level in dB
SNR=-5:10;
%choose FFT length
N_FFT=1024;
%choose guard interval length as a percentage of N_FFT
NGI=1/4;
%choose number of used carriers
N_used=N_FFT/2;
N_GI=round(NGI*N_FFT);
Nofdm=N_FFT+N_GI;
Nsym = 500;
%frequency offset
nSTOs=50;
CFO=0.4;

Fs = 11.2e6; % Частота дискретизации сигнала
Tofdm = Nofdm/Fs; % Длительность OFDM сигнала в сек
ppm = 50*10^-6; % Стабильность тактового гененратора
Nppm = Fs*ppm; % Уход в отсчетах в сек

%% Канал 1
% PathDelays = [0 50e-9 110e-9 170e-9 290e-9 310e-9]; % задержки лучей
% U = [0 -3 -10 -18 -26 -32];
% Fd= [4.6 23 139]; % Доплеровское смещение в Гц
%% Канал 2
% PathDelays = [0 110e-9 190e-9 410e-9]; % задержки лучей
% U = [0 -9.7 -19.2 -22.8];
% Fd= [4.6 23 139]; % Доплеровское смещение в Гц
%% Канал 3
PathDelays = [0 310e-9 710e-9 1090e-9 1730e-9 2510e-9]; % задержки лучей
U = [0 -1.0 -9.0 -10.0 -15.0 -20.0];
Fd= [4.6 23 139]; % Доплеровское смещение в Гц

%% preamble Minn
B=2;
reS=(-1).^round(rand(1,round(N_used/B)));
imS=1i.*(-1).^round(rand(1,round(N_used/B)));
s1=reS+imS;
s2=fliplr(s1);
tx_signal_noGI=[s1 s2 conj(s1) conj(s2)];
tx_signal = zeros(1,Nsym*N_FFT);
tx_signal_GI=[tx_signal_noGI(1,N_FFT-N_GI+1:N_FFT) tx_signal_noGI];
tx_signal=repmat(tx_signal_GI,1,Nsym);

%% Модель канала
OutChan = complex(zeros(Nsym*N_FFT,1),zeros(Nsym*N_FFT,1));
rand('seed',0);
if flChann == 1
tic
%T = [0 PathDelays(OptChan)];
%U = [0 0];
T = PathDelays;
chan = rayleighchan(1/Fs,Fd(3)/2,T,U);
chan.DopplerSpectrum=doppler.ajakes; %flat ajakes
chan.ResetBeforeFiltering=0;
OutChan=filter(chan,tx_signal);
t = toc;
fprintf('Сигнал искажен каналом...\n Время выполнения %4.2f секунд\n', t);
else
OutChan = tx_signal;
end

%% Добавление смещение частоты и времени
tic;
delay(1,1) = Nppm*Tofdm;
L = Nofdm;
for n = 2:Nsym
delay(1,n) = delay(1,n-1) + Nppm*Tofdm;
int_delay = fix(delay);
fract_delay = delay - int_delay;
end

MaxDelay = fix(max(int_delay));

OutChan(1,end:end+MaxDelay) = 0;
buff_delay = complex(zeros(1, 2*L), zeros(1, 2*L));
if flTimeShift == 1
for n = 1:Nsym
buff_delay = circshift(buff_delay, [0 -(L+int_delay(1,n))]);
buff_delay(1,L-int_delay(1,n)+1:end) = OutChan(1,(n-1)*L+1:n*L+int_delay(1,n));

re_buff_delay = sig_delay(real(buff_delay),1/Fs,fract_delay(1,n)/Fs);
im_buff_delay = sig_delay(imag(buff_delay),1/Fs,fract_delay(1,n)/Fs);

buff_delay = re_buff_delay + 1i*im_buff_delay;

y_STO(1,(n-1)*L+1:n*L) = buff_delay(1,L+1:2*L);
end
else
y_STO = OutChan;
end

if nSTOs>=0, y_STO=[y_STO(1,nSTOs+1:end) zeros(1,nSTOs)];
else
y_STO=[zeros(1,-nSTOs) y_STO(1,1:end+nSTOs)];
end

nn=0:length(y_STO)-1; y_CFO_STO = y_STO.*exp(1i*2*pi*CFO*nn/N_FFT);
t = toc;
fprintf('Осуществлен сдвиг времени и частоты сигнала...\n Время выполнения %4.2f секунд\n', t);

%%
buff = [];
for k = 1:length(SNR)
recvd_signal=awgn(y_CFO_STO, SNR(k), 'measured', 'dB');

for n = 1:Nsym
%%
Len_all=N_FFT; %
recvd_signal_zeropad=[zeros(1,Len_all) recvd_signal(1,(n-1)*Nofdm+1:n*Nofdm) zeros(1,Len_all) zeros(1,Len_all)];
for d=1:length(recvd_signal_zeropad)-Len_all-1
P_Park(d)=sum(recvd_signal_zeropad(d+(1:N_FFT/2)).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2)))); %fliplr
E_Park(d)=sum(abs(recvd_signal_zeropad(d+(1:N_FFT))).^2);
Q_Park(d)=sum(conj(recvd_signal_zeropad(d+(1:N_FFT/2))).*fliplr(recvd_signal_zeropad(d+N_FFT/2+(1:N_FFT/2))));
end
Park_metric=(abs(P_Park).^2)./(E_Park).^2;
buff = [buff Park_metric];
plot(Park_metric/max(Park_metric))

[~, maxipos1]=max(Park_metric-N_GI);
t_est_Park(1,n)=Nofdm - maxipos1;
FreqOffset(1,n) = angle(Q_Park(maxipos1-N_GI))/pi;
end

mse_CFO(k) = sqrt(sum(abs(FreqOffset - CFO).^2)/Nsym)
if flTimeShift == 1
mse_STO(k) = sqrt(sum(abs(diff(t_est_Park - nSTOs) - diff(int_delay(1,)).^2)/Nsym) %
else
mse_STO(k) = sqrt(sum(abs(t_est_Park - nSTOs).^2)/Nsym)
end

end

CODE

function sig_d = sig_delay(sig,T,dt)
% Задержка сигнала по времени
% sig - отсчеты сигнала
%T - шаг дискретизации
%dt - задержка по времени
%sig_d - задержанный сигнал
if dt==0
sig_d=sig;
return
end

N=numel(sig);
sig_d(1:N)=0;
a=fix(dt/T);
b=dt/T-a;

x=-b;
y0=0;y1=0;y2=0;y3=0;

for m=1:N
if m-a-2>0
y0=sig(m-a-2);
end
if m-a-1>0
y1=sig(m-a-1);
end
if m-a>0
y2=sig(m-a);
end
if (m-a+1>0)&&(m-a+1<=N)
y3=sig(m-a+1);
end
a3=(y3-y0)/6+(y1-y2)/2;
a1=(y3-y1)/2-a3;
a2=y3-y2-a1-a3;
a0=y2;

sig_d(m)=((a3*x+a2)*x+a1)*x+a0;
end
end

но чет у меня какие-то сомнения в адекватности модели