Ivan55 0 14 июля, 2017 Опубликовано 14 июля, 2017 (изменено) · Жалоба Добрый день! Решил я поразбираться с синхрой по преамбуле(тренировочным символам) реализовал метод Парка, не знаю правильно или нет замоделил, подскажите если что не так. Хочу померить MSE оценок по данной преамбуле, как это можно правильно сделать? 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 Изменено 14 июля, 2017 пользователем Ivan55 Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Ivan55 0 14 июля, 2017 Опубликовано 14 июля, 2017 · Жалоба Если быть точнее то я щас делаю вот, так 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 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 но чет у меня какие-то сомнения в адекватности модели Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
KalashKS 0 14 июля, 2017 Опубликовано 14 июля, 2017 · Жалоба Уберите защитный интервал из обучающей оследовательности. Он вам дает лишний пик в метрике. При ее использовании для обнаружения он будет мешаться. Параметр fd у объекта rayleighchan задает не сдвиг по частоте, а доплеровское рассеяние, моделируя одновременный приход множества лучей с разной доплеровской частотой, объединяя их в один кластер. По факту генерится комплексный гауссовский случайный процесс с СПМ, заданной параметром DopplerSpectrum и умножается на входной сигнал. Это я на всякий случай написал, может вам это и нужно было. В принципе дальше по коду частотный сдвиг моделируется отдельно. В остальном криминала не вижу. Качество оценок оцениваете в лоб как срадений квадрат разности между тем, что намеряли и тем, что задавали в модели. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
Ivan55 0 14 июля, 2017 Опубликовано 14 июля, 2017 (изменено) · Жалоба Уберите защитный интервал из обучающей оследовательности. Он вам дает лишний пик в метрике. При ее использовании для обнаружения он будет мешаться. да я в курсе. Но без защитного интервала там вобще два побочных пика получается Вот в статье Парка10.1109_LCOMM.2003.812181.pdf у него тоже этот пик тока с другой стороны, видимо потому что суммирование в другом порядке делал отсюда я и предположил что префикс нужен А в жизни преамбулы разве без защитного передают? Изменено 14 июля, 2017 пользователем Ivan55 Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться
KalashKS 0 14 июля, 2017 Опубликовано 14 июля, 2017 · Жалоба В вашей модели у меня не получилось. Один хороший пик. В статье нарисовали два, но дальше на эту тему не наспространялись, а жаль. Сам я этот алгоритм не использовал. Все, что сложнее автокоррелятора не очень люблю, очень сложная схема получается. Может на DSP нормально можно реализовать. ЗИ в преамбуле нужен, если она же будет использоваться для оценивания канала. В остальных случаях он скорее мешает. Цитата Поделиться сообщением Ссылка на сообщение Поделиться на другие сайты Поделиться