Commit 847a4203 authored by Florian Kaltenberger's avatar Florian Kaltenberger

moving project CORRIDOR to extras

git-svn-id: http://svn.eurecom.fr/openair4G/trunk@7399 818b1a75-f10b-46b9-bf7c-635c3b92a50f
parent fd45609d
function La = LanczosKernel(a,x0,samples)
min = floor(x0) - a;
max = ceil(x0) + a;
min(min<1)=1;
max(max>size(samples,2))=size(samples,2);
La=zeros(size(x0,2),1);
for i=1:size(x0,2)
La(i)=sinc((-x0(i)+(min(i):max(i)))).*sinc(((-x0(i)+(min(i):max(i)))/a))*samples(min(i):max(i))';
end
\ No newline at end of file
function received_f = OFDM_RX(received,num_carriers,useful_carriers,prefix_length,num_symbols_frame)
nant = size(received,2);
ofdm_symbol_length = num_carriers + prefix_length;
received_f = zeros(num_symbols_frame,useful_carriers,nant);
for j=0:num_symbols_frame-1;
ifblock=received(j*ofdm_symbol_length+(1:ofdm_symbol_length),:);
ifblock(1:prefix_length,:)=[];
fblock=fft(ifblock,[],1)/sqrt(num_carriers);
received_f(j+1,:,:) = [fblock(2:useful_carriers/2+1,:); fblock(end-useful_carriers/2+1:end,:)];
end
\ No newline at end of file
function [sig,sig_length] = OFDM_TX(num_carriers,num_zeros,prefix_length,input)
% OFDM Transmitter - DC removed
% sig is the output signal
% length is the length of the output signal
% num_carriers - number of sub-carriers (power of 2)
% num_zeros - number of zeros minus 1 (DC) in output spectrum (odd)
% prefix_length - length of cyclic prefix
% input - input dimensions (length = number_carriers - num_zeros - 1)
if (length(input) + num_zeros + 1 ~= num_carriers)
fprintf('error in lengths\n');
return;
end
ext_input = [0 input(1:length(input)/2) zeros(1,num_zeros) input((1+length(input)/2) : length(input))];
output_1 = ifft(ext_input);
sig = [output_1((num_carriers - prefix_length + 1) : num_carriers) output_1];
sig_length = length(sig);
function [sig, sig_f] = OFDM_TX_FRAME(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length)
% sig - output signal
% sig_length - output signal length
% num_carriers - number of sub-carriers
% num_zeros - number of zero carriers minus 1 (DC)
% prefix_length - length of cyclic prefix
% num_symbols_frame - number of symbols per OFDM frame
% preamble_length - length of 4-QAM preamble
num_useful_carriers = num_carriers - num_zeros -1;
sig = zeros(1,(num_carriers+prefix_length)*num_symbols_frame);
sig_f = zeros(num_symbols_frame,num_useful_carriers);
for k=1:preamble_length
QAM4_preamble = QAM_MOD(4,floor(256*abs(rand(1,num_useful_carriers/4))));
sig((k-1)*(num_carriers+prefix_length)+1:k*(num_carriers+prefix_length)) = OFDM_TX(num_carriers,num_zeros,prefix_length,QAM4_preamble);
sig_f(k,:) = QAM4_preamble;
end
for k=preamble_length+1:num_symbols_frame
QAM_data = QAM_MOD(256,floor(256*abs(rand(1,num_useful_carriers))));
sig((k-1)*(num_carriers+prefix_length)+1:k*(num_carriers+prefix_length)) = OFDM_TX(num_carriers,num_zeros,prefix_length,QAM_data);
sig_f(k,:) = QAM_data;
end
function [sig, sig_f] = OFDM_TX_FRAME_MIMO(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length,num_ant)
% sig - output signal
% sig_length - output signal length
% num_carriers - number of sub-carriers
% num_zeros - number of zero carriers minus 1 (DC)
% prefix_length - length of cyclic prefix
% num_symbols_frame - number of symbols per OFDM frame
% preamble_length - length of 4-QAM preamble
% numant - number of antennas
num_useful_carriers = num_carriers - num_zeros -1;
if (num_ant ==1)
t_dec = 1;
f_dec = 1;
elseif (num_ant ==2)
t_dec = 1;
f_dec = 2;
elseif (num_ant == 4)
t_dec = 2;
f_dec = 4;
else
error('Only 1, 2 or 4 antennas supported');
end
sig = zeros(num_ant,(num_carriers+prefix_length)*num_symbols_frame);
sig_f = zeros(num_ant,num_symbols_frame,num_useful_carriers);
for a=1:num_ant
for k=(floor((a-1)/2)+1):t_dec:preamble_length
QAM4_preamble = zeros(1,num_useful_carriers);
QAM4_preamble(2*mod(a-1,2)+1:f_dec:num_useful_carriers) = QAM_MOD(4,floor(4*abs(rand(1,num_useful_carriers/f_dec))));
sig(a,(k-1)*(num_carriers+prefix_length)+1:k*(num_carriers+prefix_length)) = OFDM_TX(num_carriers,num_zeros,prefix_length,QAM4_preamble);
sig_f(a,k,:) = QAM4_preamble;
end
for k=preamble_length+1:num_symbols_frame
QAM_data = QAM_MOD(256,floor(4*abs(rand(1,num_useful_carriers))));
sig(a,(k-1)*(num_carriers+prefix_length)+1:k*(num_carriers+prefix_length)) = OFDM_TX(num_carriers,num_zeros,prefix_length,QAM_data);
sig_f(a,k,:) = QAM_data;
end
end
function [sig,sig_length] = QAM_MOD(size,input)
% sig - output symbols
% size - modulation size (4,16,256)
% input - vector of bytes to be modulated
AM2 = [-1 1];
AM4 = [-3 -1 1 3]; AM4 = 2*AM4/sqrt(AM4*AM4');
AM16 = [-15 -13 -11 -9 -7 -5 -3 -1 1 3 5 7 9 11 13 15]; AM16 = 4*AM16/sqrt(AM16*AM16');
sig = zeros(1,length(input));
sig_length = length(input);
for l=1:length(input)
if (size == 256)
sig(l) = (AM16(1+ floor((input(l)/16))) + sqrt(-1)*AM16(1+rem(input(l),16)))/sqrt(2);
elseif (size == 16)
sig(l) = (AM4(1+rem(floor(input(l)/4),4)) + sqrt(-1)*AM4(1+rem(input(l),4)))/sqrt(2);
elseif (size == 4)
sig(l) = (AM2(1+rem(floor(input(l)/2),2)) + sqrt(-1)*AM2(1+rem(input(l),2)))/sqrt(2);
end
end
close all;
clear all;
%n_carriers=1;% 1 for UHF files, 2 for 2.6GHz files
%file='E:\EMOS\corridor\postprocessed data\eNB_data_UHF_20140519_run4.mat'; % mat file
%file='E:\EMOS\corridor\postprocessed data\eNB_data_20140331_UHF_run2.mat'; % mat file
n_carriers=2;
n_trials=2;
n_runs=1;
%file='E:\EMOS\corridor\postprocessed data\eNB_data_20140331_2.6GHz_run2.mat'; % mat file
%file='eNB_data_UHF_20140519_run4.mat'; % mat file
file='eNB_data_20140519_2.6GHz_run1.mat'; % mat file
post_processed_data=load(file, 'PDD_totala','PDP_totala','delay_doppler_profile_beforea','delay_doppler_profile_duringa','delay_doppler_profile_aftera');
%post_processed_data=load(file, 'PDD_totala','PDP_totala');
if(n_carriers==2)
post_processed_data=load(file, 'PDD_totala','PDP_totala','delay_doppler_profile_beforea','delay_doppler_profile_duringa','delay_doppler_profile_aftera','PDD_totalb','PDP_totalb','delay_doppler_profile_beforeb','delay_doppler_profile_duringb','delay_doppler_profile_afterb');
end
PDDta=post_processed_data(1,1).PDD_totala;
PDPta=post_processed_data(1,1).PDP_totala;
delay_doppler_profile_beforea=post_processed_data(1,1).delay_doppler_profile_beforea;
delay_doppler_profile_duringa=post_processed_data(1,1).delay_doppler_profile_duringa;
delay_doppler_profile_aftera=post_processed_data(1,1).delay_doppler_profile_aftera;
if(n_carriers==2)
PDDtb=post_processed_data(1,1).PDD_totalb;
PDPtb=post_processed_data(1,1).PDP_totalb;
delay_doppler_profile_beforeb=post_processed_data(1,1).delay_doppler_profile_beforeb;
delay_doppler_profile_duringb=post_processed_data(1,1).delay_doppler_profile_duringb;
delay_doppler_profile_afterb=post_processed_data(1,1).delay_doppler_profile_afterb;
end
if n_trials==1
if n_runs==1
block_before=50;
block_during=90;
block_after=130;
end
if n_runs==2
if n_carriers==1% we have changed the orientation of the antennas for the UHF channel in Trial 1 Run 2
block_before=60;
block_during=155;
block_after=190;
end
if n_carriers==2
block_before=60;
block_during=107;
block_after=140;
end
end
end
if n_trials==2
if n_runs==1
block_before=50;
block_during=91;
block_after=140;
end
if n_runs==2
block_before=45;
block_during=77;
block_after=120;
end
if n_runs==3
block_before=45;
block_during=83;
block_after=120;
end
if n_runs==4
block_before=34;
block_during=43;
block_after=90;
end
end
%% Doppler spectrum (choose the block you want to read)
block = 20;
figure(1)
for i=1:size(PDDta,3)
for j=1:size(PDDta,4)
subplot(size(PDDta,3),size(PDDta,4),(i-1)*size(PDDta,4)+j)
F=-(100*120/2-1)*7.68E6/(2*100*120/2)/1280:7.68E6/(100*120/2)/1280:(100*120/2-1)*7.68E6/(2*100*120/2)/1280;
if (n_carriers==2)
F=-(50*120/2-1)*30.72E6/(2*50*120/2)/5120:30.72E6/(50*120/2)/5120:(50*120/2-1)*30.72E6/(2*50*120/2)/5120;
end
plot(F,10*log(PDDta(:,block,i,j)));
ylabel('power [dB]')
xlabel('Doppler shift [Hz]')
end
end
if(n_carriers==2)
figure(2)
for i=1:size(PDDtb,3)
for j=1:size(PDDtb,4)
subplot(size(PDDtb,3),size(PDDtb,4),(i-1)*size(PDDtb,4)+j)
F=-(50*120/2-1)*15.36E6/(2*50*120/2)/2560:15.36E6/(50*120/2)/2560:(50*120/2-1)*15.36E6/(2*50*120/2)/2560;
plot(F,10*log(PDDtb(:,block,i,j)));
end
end
end
%% Power Delay Profile (choose the frame you want to read)
frame = 3000;
figure(3)
for i=1:size(PDDta,3)
for j=1:size(PDDta,4)
T=1:1:(size(PDPta,1));
tau=linspace(0,300/4/4.5E6,300);
T=1:1:(size(PDPta,1));
if n_carriers==2
tau=linspace(0,1200/4/18E6,1200/4);
end
subplot(size(PDDta,3),size(PDDta,4),(i-1)*size(PDDta,4)+j)
plot(10*log10(PDPta(frame,:,i,j)));
xlabel('delay [s]')
ylabel('time [*10 ms]')
end
end
if(n_carriers==2)
figure(4)
for i=1:size(PDDtb,3)
for j=1:size(PDDtb,4)
T=1:1:(size(PDPtb,1));
tau=linspace(0,600/4/18E6,600/4);
subplot(size(PDDtb,3),size(PDDtb,4),(i-1)*size(PDDtb,4)+j)
plot(10*log10(PDPtb(frame,:,i,j)));
xlabel('delay [s]')
ylabel('time [*10 ms]')
end
end
end
%% Total doppler spectrum in pseudocolor plot
doppler_profile_figures_dir = 'E:\byiringi\Matlab Plots\Doppler Shift pcolor plots new sync\';
h=figure(5);
hold off
for i=1:size(PDDta,3)
for j=1:size(PDDta,4)
T=1:1:size(PDDta,2);
F=-(100*120/2-1)*7.68E6/(2*100*120/2)/1280:7.68E6/(100*120/2)/1280:(100*120/2-1)*7.68E6/(2*100*120/2)/1280;
filename=sprintf('Trial %d Run %d UHF.fig',n_trials,n_runs);
if(n_carriers==2)
F=-(50*120/2-1)*30.72E6/(2*50*120/2)/5120:30.72E6/(50*120/2)/5120:(50*120/2-1)*30.72E6/(2*50*120/2)/5120;
filename=sprintf('Trial %d Run %d 2.6 GHz Carrier 2a.fig',n_trials,n_runs);
end
subplot(size(PDDta,3),size(PDDta,4),(i-1)*size(PDDta,4) + j);
pcolor(T,F,10*log10( PDDta(:,:,i,j)));
shading flat
%colormap hot
bara=colorbar;
%ylim([])
%xlim([])
xlabel('time [s]')
ylabel('Doppler shift [Hz]')
ylabel(bara,'Power [dB]')
end
end
saveas(h,strcat(doppler_profile_figures_dir, filename));
% for i=1:size(PDDta,1)
% for j=1:size(PDDta,2)
% if 10*log10(PDDta(i,j,1,1))<115
% PDDta(i,j,1,1)=10^11.5;
% end
% end
%
% end
% figure (15)
%
% T=1:1:size(PDDta,2);
% F=-(100*120/2-1)*7.68E6/(2*100*120/2)/1280:7.68E6/(100*120/2)/1280:(100*120/2-1)*7.68E6/(2*100*120/2)/1280;
% if(n_carriers==2)
% F=-(50*120/2-1)*30.72E6/(2*50*120/2)/5120:30.72E6/(50*120/2)/5120:(50*120/2-1)*30.72E6/(2*50*120/2)/5120;
% end
%
%
%
%
% pcolor(T,F,10*log10( PDDta(:,:,1,1)));
% shading flat
% colormap hot
% bara=colorbar;
% %ylim([])
% %xlim([])
% xlabel('time [s]')
% ylabel('Doppler shift [Hz]')
if(n_carriers==2)
h=figure(6);
for i=1:size(PDDtb,3)
for j=1:size(PDDtb,4)
T=1:1:size(PDDtb,2);
F=-(50*120/2-1)*15.36E6/(2*50*120/2)/2560:15.36E6/(50*120/2)/2560:(50*120/2-1)*15.36E6/(2*50*120/2)/2560;
filename=sprintf('Trial %d Run %d 2.6 GHz Carrier 2b.fig',n_trials,n_runs);
subplot(size(PDDtb,3),size(PDDtb,4),(i-1)*size(PDDtb,4) + j);
pcolor(T,F,10*log10( PDDtb(:,:,i,j)));
shading flat
barb=colorbar;
%colormap hot
%ylim([])
%xlim([])
xlabel('time [s]')
ylabel('Doppler shift [Hz]')
ylabel(barb,'Power [dB]')
end
end
saveas(h,strcat(doppler_profile_figures_dir, filename));
end
%% Total Power Delay Profile in pseudocolor
power_delay_profile_figures_dir = 'E:\byiringi\Matlab Plots\PDP pcolor plots new sync\';
h=figure(7);
for i=1:size(PDDta,3)
for j=1:size(PDDta,4)
tau=linspace(0,300/4/4.5E6,300/4);
T=1:1:(size(PDPta,1));
filename=sprintf('Trial %d Run %d UHF.fig',n_trials,n_runs);
if n_carriers==2
tau=linspace(0,1200/4/18E6,1200/4);
filename=sprintf('Trial %d Run %d 2.6 GHz Carrier 2a.fig',n_trials,n_runs);
end
subplot(size(PDDta,3),size(PDDta,4),(i-1)*size(PDDta,4)+j)
pcolor(tau,T,10*log10(PDPta(:,:,i,j)));
bara=colorbar;
shading flat
%colormap hot
xlabel('delay [s]')
ylabel('time [*10 ms]')
ylabel(bara,'Power [dB]')
end
end
saveas(h,strcat(power_delay_profile_figures_dir, filename));
% for i=1:size(PDPta,1)
% for j=1:size(PDPta,2)
% if 10*log10(PDPta(i,j,1,1))<57
% PDPta(i,j,1,1)=10^5.7;
% end
% end
%
% end
% figure (17)
% tau=linspace(0,300/4/4.5E6,300/4);
% T=1:1:(size(PDPta,1));
% if n_carriers==2
% tau=linspace(0,1200/4/18E6,1200/4);
% end
%
% pcolor(tau,T,10*log10(PDPta(:,:,1,1)));
% bara=colorbar;
% shading flat
% colormap hot
% xlabel('delay [s]')
% ylabel('time [*10 ms]')
if(n_carriers==2)
h=figure(8);
for i=1:size(PDDtb,3)
for j=1:size(PDDtb,4)
tau=linspace(0,600/4/9E6,600/4);
T=1:1:(size(PDPtb,1));
filename=sprintf('Trial %d Run %d 2.6 GHz Carrier 2b.fig',n_trials,n_runs);
subplot(size(PDDtb,3),size(PDDtb,4),(i-1)*size(PDDtb,4)+j)
pcolor(tau,T,10*log10(PDPtb(:,:,i,j)));
barb=colorbar;
shading flat
%colormap hot
xlabel('delay [s]')
ylabel('time [*10 ms]')
ylabel(barb,'Power [dB]')
end
end
saveas(h,strcat(power_delay_profile_figures_dir, filename));
end
%% Delay Doppler Spectrum before, during and after the passing of the train
figure(9)
tau=linspace(0,300/4/4.5E6,300/4);
F=-(100*120/2-1)*7.68E6/(2*100*120/2)/1280:7.68E6/(100*120/2)/1280:(100*120/2-1)*7.68E6/(2*100*120/2)/1280;
if(n_carriers==2)
tau=linspace(0,1200/4/18E6,1200/4);
F=-(50*120/2-1)*30.72E6/(2*50*120/2)/5120:30.72E6/(50*120/2)/5120:(50*120/2-1)*30.72E6/(2*50*120/2)/5120;
end
subplot(1,3,1)
pcolor(tau,F,10*log10(delay_doppler_profile_beforea(:,:)))
shading flat
colormap hot
bar1=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar1,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for UHF-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_before));
if(n_carriers==2)
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 1-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_before));
end
subplot(1,3,2)
pcolor(tau,F,10*log10(delay_doppler_profile_duringa(:,:)))
shading flat
colormap hot
bar2=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar2,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for UHF-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_during));
if(n_carriers==2)
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 1-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_during));
end
subplot(1,3,3)
pcolor(tau,F,10*log10(delay_doppler_profile_aftera(:,:)))
shading flat
colormap hot
bar3=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar3,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for UHF-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_after));
if(n_carriers==2)
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 1-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_after));
end
if(n_carriers==2)
figure(10)
tau=linspace(0,600/4/9E6,600/4);
F=-(50*120/2-1)*15.36E6/(2*50*120/2)/2560:15.36E6/(50*120/2)/2560:(50*120/2-1)*15.36E6/(2*50*120/2)/2560;
subplot(1,3,1)
pcolor(tau,F,10*log10(delay_doppler_profile_beforeb(:,:)))
shading flat
colormap hot
bar4=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar4,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 2-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_before));
subplot(1,3,2)
pcolor(tau,F,10*log10(delay_doppler_profile_duringb(:,:)))
shading flat
bar5=colorbar;
colormap hot
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar5,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 2-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_during));
subplot(1,3,3)
pcolor(tau,F,10*log10(delay_doppler_profile_afterb(:,:)))
shading flat
bar6=colorbar;
colormap hot
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
ylabel(bar6,'Power [dB]')
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 2-Trial %d-Run %d-Block %d ',n_trials,n_runs,block_after));
end
\ No newline at end of file
function H_dpss = dpss_smooth(H,V,Dopt,f_start)
alpha = V(f_start:4:end,1:Dopt)'*H.';
H_dpss = V(:,1:Dopt)*alpha;
\ No newline at end of file
close all
clear all
enable_plots=2; %eanbles figures
%% preload and init data
addpath('../../../openair1/PHY/LTE_REFSIG');
primary_synch; %loads the primary sync signal
pss_t = upsample(primary_synch0_time,4);
%load('E:\EMOS\corridor\ofdm_pilots_sync_2048_v7.mat');
load('ofdm_pilots_sync_30MHz.mat');
ofdm_symbol_length = num_carriers + prefix_length;
frame_length = ofdm_symbol_length*num_symbols_frame;
useful_carriers = num_carriers-num_zeros-1;
%filename = 'E:\EMOS\corridor2\eNB_data_20140321_184441.EMOS';
filename = 'D:\711MHz\eNB_data_20140324_113931.EMOS';
samples_slot = 7680/2;
slots_per_frame = 20;
nframes = 100;
nant=3;
d = dir(filename);
nblocks = floor(d.bytes/(samples_slot*slots_per_frame*nframes*nant*4));
PDP_total = zeros(nblocks*nframes,useful_carriers/4);
%% main loop
fid = fopen(filename,'r');
block = 1;
while ~feof(fid)
fprintf(1,'Processing block %d of %d',block,nblocks);
[v,c]=fread(fid, samples_slot*slots_per_frame*nframes*nant*2, 'int16',0,'ieee-le');
if (c==0)
break
end
v1 = double(v(1:2:end))+1j*double(v(2:2:end));
v2 = zeros(samples_slot*slots_per_frame*nframes,nant);
for slot=1:slots_per_frame*nframes
for a=1:nant
v2((slot-1)*samples_slot+1:slot*samples_slot,a) = ...
v1((slot-1)*samples_slot*nant+(a-1)*samples_slot+1:...
(slot-1)*samples_slot*nant+ a *samples_slot,1);
end
end
if enable_plots>=2
figure(1)
plot(20*log10(abs(fftshift(fft(v2)))))
end
%% frame start detection
[corr,lag] = xcorr(v2(:,1),pss_t);
%[m,idx]=max(abs(corr));
[m,idx]=peaksfinder(corr,frame_length);
if enable_plots>=2
figure(2);
hold off
plot(lag,abs(corr));
hold on
plot(lag(idx),m,'ro')
end
%%
for i=1:size(idx,2)-1; % the last frame is not complite
fprintf(1,'.');
frame_start = lag(idx(i))-prefix_length;
% frame_start = lag(i) - prefix_length;
%% ofdm receiver
received_f = OFDM_RX(v2(frame_start:frame_start+frame_length,:),num_carriers,useful_carriers,prefix_length,num_symbols_frame);
%% channel estimation
H=conj(squeeze(f3(1,3:2:end,1:4:end))).*received_f(3:2:end,1:4:end,1);
Ht = ifft(H,[],2);
PDP = mean(abs(Ht).^2,1);
PDP_total((block-1)*nframes+i+1,:) = PDP;
if enable_plots>=1
figure(3)
surf(20*log10(abs(Ht)))
xlabel('time [OFDM symbol]')
ylabel('delay time [samples]')
zlabel('power [dB]')
shading interp
figure(4)
plot(10*log10(PDP))
xlabel('delay time [samples]')
ylabel('power [dB]')
end
end
fprintf(1,'\n');
block = block+1;
end
fclose(fid);
close all
clear all
enable_plots=2; %enables figures
%% preload and init data
addpath('../../../openair1/PHY/LTE_REFSIG');
primary_synch; %loads the primary sync signal
pss1_t = upsample(primary_synch0_time,4*4);
pss2_t = upsample(primary_synch0_time,4*2);
%load('E:\EMOS\corridor\ofdm_pilots_sync_2048_v7.mat');
load('ofdm_pilots_sync_30MHz.mat');
%filename = 'E:\EMOS\corridor\eNB_data_20140319_133327.EMOS';
filename = 'D:\2.6GHz\eNB_data_20140324_171904.EMOS';
nb_rb1 = 100; %this can be 25, 50, or 100
num_carriers1 = 2048/100*nb_rb1;
num_zeros1 = num_carriers1-(12*nb_rb1+1);
prefix_length1 = num_carriers1/4; %this is extended CP
ofdm_symbol_length1 = num_carriers1 + prefix_length1;
frame_length1 = ofdm_symbol_length1*num_symbols_frame;
useful_carriers1 = num_carriers1-num_zeros1-1;
nb_rb2 = 50; %this can be 25, 50, or 100
num_carriers2 = 2048/100*nb_rb2;
num_zeros2 = num_carriers2-(12*nb_rb2+1);
prefix_length2 = num_carriers2/4; %this is extended CP
ofdm_symbol_length2 = num_carriers2 + prefix_length2;
frame_length2 = ofdm_symbol_length2*num_symbols_frame;
useful_carriers2 = num_carriers2-num_zeros2-1;
nant1 = 2;
nant2 = 2;
samples_slot1 = 7680*2;
samples_slot2 = 7680;
samples_slot_agg = nant1*samples_slot1 + nant2*samples_slot2;
nframes = 50;
slots_per_frame = 20;
d = dir(filename);
nblocks = floor(d.bytes/(samples_slot_agg*slots_per_frame*nframes*4));
PDP1_total = zeros(nblocks*nframes,useful_carriers1/4);
PDP2_total = zeros(nblocks*nframes,useful_carriers2/4);
%% main loop
fid = fopen(filename,'r');
block = 1;
while ~feof(fid)
fprintf(1,'Processing block %d of %d',block,nblocks);
%%
[v,c]=fread(fid, 2*samples_slot_agg*slots_per_frame*nframes, 'int16',0,'ieee-le');
block = block+1;
if (c==0)
break
end
v0 = double(v(1:2:end,:))+1j*double(v(2:2:end,:));
v1 = zeros(samples_slot1*slots_per_frame*nframes,nant1);
v2 = zeros(samples_slot2*slots_per_frame*nframes,nant2);
for slot=1:slots_per_frame*nframes
for a1=1:nant1
v1((slot-1)*samples_slot1+1:slot*samples_slot1,a1) = ...
v0((slot-1)*samples_slot_agg+(a1-1)*samples_slot1+1:...
(slot-1)*samples_slot_agg+ a1 *samples_slot1,1);
end
for a2=1:nant2
v2((slot-1)*samples_slot2+1:slot*samples_slot2,a2) = ...
v0((slot-1)*samples_slot_agg+nant1*samples_slot1+(a2-1)*samples_slot2+1:...
(slot-1)*samples_slot_agg+nant1*samples_slot1+ a2 *samples_slot2,1);
end
end
if enable_plots>=2
figure(1)
plot(20*log10(abs(fftshift(fft(v1)))))
figure(2)
plot(20*log10(abs(fftshift(fft(v2)))))
end
%% frame start detection
[corr1,lag1] = xcorr(v1(:,1),pss1_t);
[corr2,lag2] = xcorr(v2(:,1),pss2_t);
%[m,idx]=max(abs(corr));
[m1,idx1]=peaksfinder(corr1,frame_length1);
[m2,idx2]=peaksfinder(corr2,frame_length2);
if (enable_plots>=2)
figure(20);
hold off
plot(lag1,abs(corr1));
hold on
plot(lag1(idx1),m1,'ro')
figure(21);
hold off
plot(lag2,abs(corr2));
hold on
plot(lag2(idx2),m2,'ro')
end
%%
for i=1:size(idx1,2)-1; % the last frame is not complite
fprintf(1,'.');
%frame_start2 = lag(i) - prefix_length2;
frame_start1 = lag1(idx1(i))-prefix_length1;
frame_start2 = lag2(idx2(i))-prefix_length2;
%frame_start1 = frame_start2*2;
% ofdm receiver
received_f1 = OFDM_RX(v1(frame_start1:frame_start1+frame_length1,:),num_carriers1,useful_carriers1,prefix_length1,num_symbols_frame);
received_f2 = OFDM_RX(v2(frame_start2:frame_start2+frame_length2,:),num_carriers2,useful_carriers2,prefix_length2,num_symbols_frame);
% channel estimation (SISO)
H1=conj(squeeze(f1(1,3:2:end,1:4:end))).*received_f1(3:2:end,1:4:end,1);
H2=conj(squeeze(f2(1,3:2:end,1:4:end))).*received_f2(3:2:end,1:4:end,1);
H1t = ifft(H1,[],2);
H2t = ifft(H2,[],2);
PDP1 = mean(abs(H1t).^2,1);
PDP2 = mean(abs(H2t).^2,1);
PDP1_total((block-1)*nframes+i+1,:) = PDP1;
PDP2_total((block-1)*nframes+i+1,:) = PDP2;
if enable_plots>=1
figure(3)
surf((abs(H1t)))
xlabel('time [OFDM symbol]')
ylabel('delay time [samples]')
zlabel('power [dB]')
title('H1t')
shading interp
figure(4)
plot(10*log10(PDP1))
xlabel('delay time [samples]')
ylabel('power [dB]')
title('PDP1')
figure(30)
surf((abs(H2t)))
xlabel('time [OFDM symbol]')
ylabel('delay time [samples]')
zlabel('power [dB]')
title('H2t')
shading interp
figure(40)
plot(10*log10(PDP2))
xlabel('delay time [samples]')
ylabel('power [dB]')
title('PDP2')
end
end
fprintf(1,'\n');
end
fclose(fid);
close all
clear all
global symbols_per_slot slots_per_frame;
addpath('E:\Synchro\kaltenbe\My Documents\Matlab\dpss_chanest')
addpath('../../../openair1/PHY/LTE_REFSIG');
enable_plots=0; %enables figures
record=0; %put 1 to enable the video record of the delay doppler profile
%% preload and init data
primary_synch; %loads the primary sync signal
%load('E:\EMOS\corridor\ofdm_pilots_sync_2048_v7.mat');
load('ofdm_pilots_sync_30MHz.mat');
n_carriers = 2; % use 1 for UHF and 2 for 2.6GHz
n_trials=2;%use 1 for trial1 and 2 for trial2
n_run=1;
symbols_per_slot = 6;
slots_per_frame = 20;
sourcedir = 'F:\trials2 train extracted\';
destdir = 'F:\trials2 train processed\';
switch n_carriers
case 1,
p = init_params(25,3,4); %this can be 25, 50, or 100
pss_t = upsample(primary_synch0_time,4);
if n_trials==1
filename = sprintf('eNB_data_20140331_UHF_run%d',n_run);
else
filename = sprintf('eNB_data_UHF_20140519_run%d',n_run);
end
nframes = 100; % frames in one block
threshold = 3e+4 ; % maybe should change that !!!!
case 2,
p(1) = init_params(100,2,4);
p(2) = init_params(50,2,4);
pss_t = upsample(primary_synch0_time,4*4); % this assumes we are doing the sync on the second carrier, which is 10MHz
if (n_trials==1)
filename = sprintf('eNB_data_20140331_2.6GHz_run%d',n_run);
else
filename = sprintf('eNB_data_20140519_2.6GHz_run%d',n_run);
end
nframes = 50; % frames in one block
threshold = 3e+4 ; % maybe should change that !!!!
end
% derived parameters
samples_slot_agg = sum([p.nant_rx].*[p.samples_slot]);
num_symbols_frame = symbols_per_slot*slots_per_frame;
d = dir(fullfile(sourcedir,[filename '.EMOS']));
nblocks = floor(d.bytes/(samples_slot_agg*slots_per_frame*nframes*4));
%frequency offset
if(n_carriers==1)
if(n_trials==1)
f_offset=840;
end
if(n_trials==2)
%f_offset=;
end
end
if(n_carriers==2)
if(n_trials==1)
%f_offset=;
end
if(n_trials==2)
%f_offset=;
end
end
doppler_freq_of_max_a=zeros(1,nblocks);
doppler_freq_of_max_b=zeros(1,nblocks);
if(n_carriers==1)
fm_total=zeros(1,nblocks);%vector containing the mean Doppler Shift for each block
freqOffset_total=zeros(1,nblocks);%vector containing the mean frequency offset for each block
TGVspeed_total=zeros(1,nblocks);%vector containing the TGV speed for each block
end
PDP_totala = zeros(nblocks*nframes,p(1).useful_carriers/4,p(1).nant_tx,p(1).nant_rx);
PDD_totala = zeros(nframes*num_symbols_frame/2,nblocks,p(1).nant_tx,p(1).nant_rx);
% delay doppler spectrum
delay_doppler_profile_videoa=VideoWriter(sprintf('Trial%d_Run%d_UHF_delayDopplerProfile.avi',n_trials,n_run));%variable used to make a video of the evolution of the delay doppler profile
if n_carriers==2
delay_doppler_profile_videoa=VideoWriter(sprintf('Trial%d_Run%d_2.6GHzCarrier2a_delayDopplerProfile.avi',n_trials,n_run));
end
delay_doppler_profile_beforea=zeros(nframes*num_symbols_frame/2,p(1).useful_carriers/4);%contains the delay doppler spectrum for a block before the passing of the train
delay_doppler_profile_duringa=zeros(nframes*num_symbols_frame/2,p(1).useful_carriers/4);%contains the delay doppler spectrum for a block during the passing of the train
delay_doppler_profile_aftera=zeros(nframes*num_symbols_frame/2,p(1).useful_carriers/4);%contains the delay doppler spectrum for a block after the passing of the train
if n_carriers==2
delay_doppler_profile_videob=VideoWriter(sprintf('Trial%d_Run%d_2.6GHzCarrier2b_delayDopplerProfile.avi',n_trials,n_run));
delay_doppler_profile_beforeb=zeros(nframes*num_symbols_frame/2,p(2).useful_carriers/4);%contains the delay doppler spectrum for a block before the passing of the train
delay_doppler_profile_duringb=zeros(nframes*num_symbols_frame/2,p(2).useful_carriers/4);%contains the delay doppler spectrum for a block during the passing of the train
delay_doppler_profile_afterb=zeros(nframes*num_symbols_frame/2,p(2).useful_carriers/4);%contains the delay doppler spectrum for a block after the passing of the train
end
if n_trials==1
if n_run==1
block_before=50;
block_during=90;
block_after=130;
end
if n_run==2
if n_carriers==1% we have changed the orientation of the antennas for the UHF channel in Trial 1 Run 2
block_before=60;
block_during=155;
block_after=190;
end
if n_carriers==2
block_before=60;
block_during=107;
block_after=140;
end
end
end
if n_trials==2
if n_run==1
block_before=50;
block_during=91;
block_after=140;
end
if n_run==2
block_before=45;
block_during=77;
block_after=120;
end
if n_run==3
block_before=45;
block_during=83;
block_after=120;
end
if n_run==4
block_before=34;
block_during=43;
block_after=90;
end
end
if(n_carriers==2)
PDP_totalb = zeros(nblocks*nframes,p(2).useful_carriers/4,p(2).nant_tx,p(2).nant_rx);
PDD_totalb=zeros(nframes*num_symbols_frame/2,nblocks,p(2).nant_tx,p(2).nant_rx);
interesting_delay_doppler_profileb=zeros(nframes*num_symbols_frame/2,p(2).useful_carriers/4);%contains the delay doppler spectrum for the wanted block
if (n_trials==2)
interesting_block=60;%contains the value of one interesting block for the delay_doppler_spectrum
end
end
syncblock=0;%variable containing the number of the synchronization block
%% init DPSS parameters
max_tau = 1e-6;
for carrier=1:n_carriers
p_dpss(carrier) = init_dpss(p(carrier),max_tau);
end
%% main loop
fid = fopen(fullfile(sourcedir,[filename '.EMOS']),'r');
vStorage1 = [];
vStorage2 = [];
block = 1;
flag1 = 1;
start = 1; % Maybe 2; if it works with 1, then the variable is useless
%fseek(fid,samples_slot_agg*slots_per_frame*nframes*60*2,'bof'); %advance 30 sec
NFRAMES = 100;
if(n_carriers==2)
NFRAMES=50;
end
nframes = NFRAMES;
open(delay_doppler_profile_videoa);
if n_carriers==2
open(delay_doppler_profile_videob);
end
noise1 = zeros(nblocks*NFRAMES,p(1).nant_rx);
if n_carriers==2
noise2 = zeros(nblocks*NFRAMES,p(2).nant_rx);
end
H_power1 = zeros(nblocks*NFRAMES,p(1).nant_tx,p(1).nant_rx);
if n_carriers==2
H_power2 = zeros(nblocks*NFRAMES,p(2).nant_tx,p(2).nant_rx);
end
while ~feof(fid)
fprintf(1,'Processing block %d of %d',block,nblocks);
[v,c]=fread(fid, 2*samples_slot_agg*slots_per_frame*nframes, 'int16',0,'ieee-le');
if (c<2*samples_slot_agg*slots_per_frame*nframes)
break
end
v0 = double(v(1:2:end))+1j*double(v(2:2:end));
v1 = zeros(p(1).samples_slot*slots_per_frame*nframes,p(1).nant_rx);
if n_carriers==2
v2 = zeros(p(2).samples_slot*slots_per_frame*nframes,p(2).nant_rx);
end
for slot=1:slots_per_frame*nframes
for a=1:p(1).nant_rx
v1((slot-1)*p(1).samples_slot+1:slot*p(1).samples_slot,a) = ...
v0((slot-1)*samples_slot_agg+(a-1)*p(1).samples_slot+1:...
(slot-1)*samples_slot_agg+ a *p(1).samples_slot,1);
end
if n_carriers==2
for a=1:p(2).nant_rx
v2((slot-1)*p(2).samples_slot+1:slot*p(2).samples_slot,a) = ...
v0((slot-1)*samples_slot_agg+p(1).nant_rx*p(1).samples_slot+(a-1)*p(2).samples_slot+1:...
(slot-1)*samples_slot_agg+p(1).nant_rx*p(1).samples_slot+ a *p(2).samples_slot,1);
end
end
end
v1 = [vStorage1; v1] ;
if size(v1,1) > p(1).frame_length*nframes ;
nframes = floor(size(v1,1) / p(1).frame_length) ;
vStorage1 = v1(p(1).frame_length*nframes+1:end,:) ;
v1(p(1).frame_length*nframes + 1 : end,:) = [] ;
start = 1 ;
end
if n_carriers==2
v2 = [vStorage2; v2] ;
if size(v2,1) > p(2).frame_length*nframes ;
nframes = floor(size(v2,1) / p(2).frame_length) ;
vStorage2 = v2(p(2).frame_length*nframes+1:end,:) ;
v2(p(2).frame_length*nframes + 1 : end,:) = [] ;
start = 1 ;
end
end
if enable_plots>=2
if(n_carriers==1)
figure(1)
title('');
plot(20*log10(abs(fftshift(fft(v1)))))
end
if(n_carriers==2)
figure(1)
subplot(1,2,1);
plot(20*log10(abs(fftshift(fft(v1)))))
subplot(1,2,2);
plot(20*log10(abs(fftshift(fft(v2)))))
end
end
%% frame start detection
if flag1==1
[corr,lag] = xcorr(v1(:,1),pss_t);
% if(n_carriers==2)
% [corrb,lagb] = xcorr(v2(:,1),pss_t);
% end
%[m,idx]=max(abs(corr));
%[m,idx]=peaksfinder(corr,frame_length);
tmp = corr(nframes*slots_per_frame*p(1).samples_slot:end);
tmp2 = reshape(tmp,slots_per_frame*p(1).samples_slot,nframes);
[m,idx] = max(abs(tmp2),[],1);
% if(n_carriers==2)
% tmp = corrb(nframes*slots_per_frame*p(2).samples_slot:end);
% tmp2 = reshape(tmp,slots_per_frame*p(2).samples_slot,nframes);
% [m,idx] = max(abs(tmp2),[],1);
% end
idx(m < threshold) = [];
if size(idx,2) <= 3
flag1 = 1 ;
flag2 = 0 ;
vStorage1 = [];
vStorage2 = [];
% elseif size(idx,2) == nframes
%
% flag1 = 0;
% flag2 = 1;
else
flag1 = 0 ;
flag2 = 1 ;
syncblock=block;
end
frame_offset = round(median(idx)) - p(1).prefix_length;
% if(n_carriers==2)
% frame_offset = round(median(idx)) - p(2).prefix_length;
% end
if enable_plots>=2
figure(2)
hold off
plot(lag,abs(corr));
hold on
plot(frame_offset,m(1),'ro')
end
else
frame_offset = 0;
end
%%
if flag2 == 1
H1a=[];
if(n_carriers==2)
H1b=[];
end
fma=0;%maximum of the doppler spectrum
sa=0;
if(n_carriers==2)
fmb=0;%maximum of the doppler spectrum
sb=0;
end
max2=0;
if(n_carriers==1)
fm1=0;%First maximum of the doppler spectrum
s1=0;
fm2=0;%Second maximum of the doppler spectrum
s2=0;
max1=0;%variable containing a maximum
end
H_dpss1 = zeros(nframes*num_symbols_frame/2,p(1).useful_carriers,p(1).nant_tx,p(1).nant_rx);
if (n_carriers==2)
H_dpss2 = zeros(nframes*num_symbols_frame/2,p(2).useful_carriers,p(2).nant_tx,p(2).nant_rx);
end
for i=start:nframes;
fprintf(1,'.');
frame_start1 = (slots_per_frame*p(1).samples_slot)*(i-1)+frame_offset+1;
if n_carriers==2
frame_start2 = (slots_per_frame*p(2).samples_slot)*(i-1)+round(frame_offset/2)+1;
end
if i<nframes
%% ofdm receiver
received_f1 = OFDM_RX(v1(frame_start1:frame_start1+p(1).frame_length,:),p(1).num_carriers,p(1).useful_carriers,p(1).prefix_length,num_symbols_frame);
if n_carriers==2
received_f2 = OFDM_RX(v2(frame_start2:frame_start2+p(2).frame_length,:),p(2).num_carriers,p(2).useful_carriers,p(2).prefix_length,num_symbols_frame);
end
else
vStorage1 = [v1(frame_start1:end,:) ; vStorage1]; %%
if n_carriers==2
vStorage2 = [v2(frame_start2:end,:) ; vStorage2]; %%
end
%break
end
%disp(i);
%% noise estimation
noise1((block-1)*NFRAMES+i,:) = squeeze(10*log10(mean(abs(received_f1(1,61:end-60,:)).^2)));
if (n_carriers==2)
noise2((block-1)*NFRAMES+i,:) = squeeze(10*log10(mean(abs(received_f2(1,61:end-60,:)).^2)));
end
%% MIMO channel estimation
if (n_carriers==1)
transmit_f1 = f3;
else
transmit_f1 = f1;
transmit_f2 = f2;
end
for carrier=1:n_carriers
if (carrier==1)
transmit_f = transmit_f1;
received_f = received_f1;
else
transmit_f = transmit_f2;
received_f = received_f2;
end
% standard channel estimation
H = zeros(num_symbols_frame/2,p(carrier).useful_carriers/4,p(carrier).nant_tx,p(carrier).nant_rx);
H_dpss = zeros(num_symbols_frame/2,p(carrier).useful_carriers/4,p(carrier).nant_tx,p(carrier).nant_rx);
for itx=1:p(carrier).nant_tx
% f_start and t_start indicate the start of the pilots in time
% and frequency according to the specifications (see .doc file).
% t_start has to be >=2, since the first symbol is the PSS.
f_start = mod(itx-1,2)*2+1;
t_start = floor((itx-1)/2)+1;
for irx=1:p(carrier).nant_rx
H(:,:,itx,irx)=conj(squeeze(transmit_f(itx,t_start:2:end,f_start:4:end))).*received_f(t_start:2:end,f_start:4:end,irx);
end
end
if (carrier==1)
H_power1((block-1)*NFRAMES+i,:,:) = 10*log10(mean(mean(abs(H).^2,1),2));
else
H_power2((block-1)*NFRAMES+i,:,:) = 10*log10(mean(mean(abs(H).^2,1),2));
end
%% dpss channel estimation
snr = zeros(size(H_power1));
for itx=1:p(carrier).nant_tx
f_start = mod(itx-1,2)*2+1;
for irx=1:p(carrier).nant_rx
snr((block-1)*NFRAMES+i,itx,irx)=squeeze(H_power1((block-1)*NFRAMES+i,itx,irx))-noise1((block-1)*NFRAMES+i,irx);
snr_rounded = min(max(round(snr((block-1)*NFRAMES+i,itx,irx)),0),30);
for it=1:num_symbols_frame/2
if (carrier==1)
H_dpss1((i-1)*num_symbols_frame/2+it,:,itx,irx) = dpss_smooth(H(it,:,itx,irx),p_dpss(carrier).V,p_dpss(carrier).Dopt(snr_rounded+1),f_start);
else
H_dpss2((i-1)*num_symbols_frame/2+it,:,itx,irx) = dpss_smooth(H(it,:,itx,irx),p_dpss(carrier).V,p_dpss(carrier).Dopt(snr_rounded+1),f_start);
end
end
end
end
%% compute delay Doppler power profiles
Ht = ifft(H,[],2)*sqrt(size(H,2)); %impulse response
PDP = mean(abs(Ht).^2,1);
PDP_all = squeeze(mean(mean(PDP,3),4));
PDD=sum(abs(fftshift(fft(Ht,[],1)/sqrt(size(Ht,1)))).^2,2);
if(carrier==1)
PDP_totala((block-1)*NFRAMES+i,:,:,:) = PDP;
Ha=H;
end
if(carrier==2)
PDP_totalb((block-1)*NFRAMES+i,:,:,:) = PDP;
Hb=H;
end
% % frequency offset correction
% Hprime=H*exp(2*i*pi*167E-6*f_offset);
% Htprime = ifft(Hprime,[],2); %impulse response
% PDPprime = mean(abs(Htprime).^2,1);
%
% PDDprime=sum(abs(fftshift(fft(Htprime,[],1))).^2,2);
% if(carrier==1)
% PDP_totala((block-1)*NFRAMES+i,:,:,:) = PDP;
% Ha=H;
% end
%
% if(carrier==2)
% PDP_totalb((block-1)*NFRAMES+i,:,:,:) = PDP;
% Hb=H;
% end
if enable_plots>=1
figure(3+3*(carrier-1))
for itx=1:p(carrier).nant_tx
for irx=1:p(1).nant_rx
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
surf(20*log10(abs(Ht(:,:,itx,irx))))
ylabel('time [OFDM symbol]')
xlabel('delay time [samples]')
zlabel('power [dB]')
shading interp
end
end
figure(4+3*(carrier-1))
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
plot(10*log10(PDP(:,:,itx,irx)))
ylim([50 80])
xlim([0 75])
xlabel('delay time [samples]')
ylabel('power [dB]')
end
end
figure(5+3*(carrier-1))
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
if(n_carriers==1)
F=-(num_symbols_frame/2-1)*7.68E6/(2*num_symbols_frame/2)/1280:7.68E6/(num_symbols_frame/2)/1280:(num_symbols_frame/2)*7.68E6/(2*num_symbols_frame/2)/1280;
end
if(n_carriers==2)
if(carrier==1)
F=-(num_symbols_frame/2-1)*30.72E6/(2*num_symbols_frame/2)/5120:30.72E6/(num_symbols_frame/2)/5120:(num_symbols_frame/2)*30.72E6/(2*num_symbols_frame/2)/5120;
end
if(carrier==2)
F=-(num_symbols_frame/2-1)*15.36E6/(2*num_symbols_frame/2)/2560:15.36E6/(num_symbols_frame/2)/2560:(num_symbols_frame/2)*15.36E6/(2*num_symbols_frame/2)/2560;
end
end
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
plot(F,10*log10(PDD(:,:,itx,irx)))
%ylim([])
%xlim([])
xlabel('F=f-ftx [Hz]')
ylabel('power [dB]')
end
end
drawnow
end
if(n_carriers==1)
if carrier==1
% adjust frame offset base on channel estimate to compensate for
% timing drift. We try to keep the peak of the impulse response at
% sample prefix_length/8.
[m,idx] = max(fft(ifft(PDP_all),p(carrier).num_carriers));
offset = idx - p(carrier).prefix_length/8;
if offset > p(carrier).prefix_length
offset = offset - p(carrier).num_carriers;
end
if abs(offset) > 5
frame_offset = frame_offset + round(offset/4);
end
end
end
if(n_carriers==2)
if carrier==2
% adjust frame offset base on channel estimate to compensate for
% timing drift. We try to keep the peak of the impulse response at
% sample prefix_length/8.
[m,idx] = max(fft(ifft(PDP_all),p(carrier).num_carriers));
offset = idx - p(carrier).prefix_length/8;
if offset > p(carrier).prefix_length
offset = offset - p(carrier).num_carriers;
end
if abs(offset) > 5
frame_offset = frame_offset + round(offset/4);
end
end
end
end
H1a=cat(1,H1a,Ha);
if(n_carriers==2)
H1b=cat(1,H1b,Hb);
end
end
%save postprocessed channels
for carrier=1:n_carriers
varname1 = sprintf('H_dpss%d',carrier);
varname2 = sprintf('H_dpss_block%d',block);
eval([varname2 '=' varname1 ';']);
filename_H = sprintf('_H_dpss_carrier%d.mat',carrier);
if exist(fullfile(destdir,[filename filename_H]),'file')
save(fullfile(destdir,[filename filename_H]),varname2,'-append');
else
save(fullfile(destdir,[filename filename_H]),varname2);
end
clear(varname2);
end
Ht1a=ifft(H1a,[],2);
PDD1a=sum(abs(fftshift(fft(Ht1a,[],1))).^2,2);
delayPDD1a=mean(mean(abs(fftshift(fft(Ht1a,[],1))).^2,3),4);
if(n_carriers==2)
Ht1b=ifft(H1b,[],2);
PDD1b=sum(abs(fftshift(fft(Ht1b,[],1))).^2,2);
delayPDD1b=mean(mean(abs(fftshift(fft(Ht1b,[],1))).^2,3),4);
end
if(enable_plots>=2)
figure(9)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
F=-(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280:7.68E6/(NFRAMES*num_symbols_frame/2)/1280:(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280;
if(n_carriers==2)
F=-(NFRAMES*num_symbols_frame/2-1)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120:30.72E6/(NFRAMES*num_symbols_frame/2)/5120:(NFRAMES*num_symbols_frame/2)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120;
end
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
title(sprintf('Doppler Spectrum for UHF-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
if n_carriers==2
title(sprintf('Doppler Spectrum for 2.6GHz Carrier 1-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
end
plot(F,10*log10(PDD1a(:,:,itx,irx)))
xlabel('F=f-ftx [Hz]')
ylabel('power [dB]')
end
end
if(n_carriers==2)
figure(10)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
F=-(NFRAMES*num_symbols_frame/2-1)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560:15.36E6/(NFRAMES*num_symbols_frame/2)/2560:(NFRAMES*num_symbols_frame/2)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560;
subplot(p(2).nant_tx,p(2).nant_rx,(itx-1)*p(2).nant_rx + irx);
title(sprintf('Doppler Spectrum for 2.6GHz Carrier 2-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
plot(F,10*log10(PDD1b(:,:,itx,irx)))
xlabel('F=f-ftx [Hz]')
ylabel('power [dB]')
end
end
end
end
if record==1
ha=figure(20);
set(gca,'nextplot','replacechildren');
set(gcf,'Renderer','zbuffer');
tau=linspace(0,p(1).useful_carriers/4/4.5E6,p(1).useful_carriers/4);
F=-(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280:7.68E6/(NFRAMES*num_symbols_frame/2)/1280:(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280;
if(n_carriers==2)
tau=linspace(0,p(1).useful_carriers/4/18E6,p(1).useful_carriers/4);
F=-(NFRAMES*num_symbols_frame/2-1)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120:30.72E6/(NFRAMES*num_symbols_frame/2)/5120:(NFRAMES*num_symbols_frame/2)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120;
end
title(sprintf('Delay Doppler Spectrum for UHF-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
if(n_carriers==2)
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 1-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
end
pcolor(tau,F,10*log10(delayPDD1a(:,:)))
shading flat
bara=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
framea = getframe(ha);
writeVideo(delay_doppler_profile_videoa,framea);
if(n_carriers==2)
hb=figure(21);
set(gca,'nextplot','replacechildren');
set(gcf,'Renderer','zbuffer');
tau=linspace(0,p(2).useful_carriers/4/9E6,p(2).useful_carriers/4);
F=-(NFRAMES*num_symbols_frame/2-1)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560:15.36E6/(NFRAMES*num_symbols_frame/2)/2560:(NFRAMES*num_symbols_frame/2)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560;
title(sprintf('Delay Doppler Spectrum for 2.6GHz Carrier 2-Trial %d-Run %d-Block %d ',n_trials,n_run,block));
pcolor(tau,F,10*log10(delayPDD1b(:,:)))
shading flat
barb=colorbar;
xlabel('delay [s]')
ylabel('Doppler shift [Hz]')
frameb = getframe(hb);
writeVideo(delay_doppler_profile_videob,frameb);
end
end
PDD_totala(:,block,:,:)=PDD1a;
if(block==block_before)
delay_doppler_profile_beforea=delayPDD1a;
end
if(block==block_during)
delay_doppler_profile_duringa=delayPDD1a;
end
if(block==block_after)
delay_doppler_profile_aftera=delayPDD1a;
end
if(n_carriers==2)
PDD_totalb(:,block,:,:)=PDD1b;
if(block==block_before)
delay_doppler_profile_beforeb=delayPDD1b;
end
if(block==block_during)
delay_doppler_profile_duringb=delayPDD1b;
end
if(block==block_after)
delay_doppler_profile_afterb=delayPDD1b;
end
end
%%
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
for i=1:NFRAMES*num_symbols_frame/2
if(10*log10(PDD1a(i,:,itx,irx))>max2)
max2=10*log10(PDD1a(i,:,itx,irx));
fma=i;
end
end
sa=sa+fma;
end
end
sa=sa/(p(1).nant_tx*p(1).nant_rx)-2999.5;
doppler_freq_of_max_a(block)=sa;
max2=0;
if(n_carriers==2)
for itx=1:p(2).nant_tx
for irx=1:p(2).nant_rx
for i=1:NFRAMES*num_symbols_frame/2
if(10*log10(PDD1a(i,:,itx,irx))>max2)
max2=10*log10(PDD1a(i,:,itx,irx));
fmb=i;
end
end
sb=sb+fmb;
end
end
sb=sb/(p(2).nant_tx*p(2).nant_rx)-2999.5;
doppler_freq_of_max_b(block)=sb;
end
%% Doppler shift, tgv speed and frequency offset for trial1 UHF run1
if(n_carriers==1)
if(n_trials==1)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
for i=1940:1960
if(10*log10(PDD1a(i,:,itx,irx))>max1)
max1=10*log10(PDD1a(i,:,itx,irx));
fm1=i;
end
end
s1=s1+fm1;
end
end
s1=s1/(p(1).nant_tx*p(1).nant_rx);
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
for i=2340:2370
if(10*log10(PDD1a(i,:,itx,irx))>max1)
max1=10*log10(PDD1a(i,:,itx,irx));
fm2=i;
end
end
s2=s2+fm2;
end
end
s2=s2/(p(1).nant_tx*p(1).nant_rx);
fm=(s2-s1)/2;
if(abs(300-fm*3/7.7715*3.6)<50)
fm_total(block)=fm;
TGVspeed_total(block)=fm*3/7.7715*3.6;
freqOffset_total(block)=abs((s1+fm)-3000.5);
end
end
end
end
fprintf(1,'\n');
block = block+1;
if (size(vStorage1,1)>=p(1).frame_length)
nframes=NFRAMES-floor((size(vStorage1,1))/(p(1).frame_length));
else
nframes=NFRAMES;
end
end
close(delay_doppler_profile_videoa);
if n_carriers==2
close(delay_doppler_profile_videob);
end
%%
if(enable_plots>=2)
figure(11)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
T=1:1:block-1;
F=-(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280:7.68E6/(NFRAMES*num_symbols_frame/2)/1280:(NFRAMES*num_symbols_frame/2-1)*7.68E6/(2*NFRAMES*num_symbols_frame/2)/1280;
if(n_carriers==2)
F=-(NFRAMES*num_symbols_frame/2-1)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120:30.72E6/(NFRAMES*num_symbols_frame/2)/5120:(NFRAMES*num_symbols_frame/2)*30.72E6/(2*NFRAMES*num_symbols_frame/2)/5120;
end
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
title(sprintf('Doppler spectrum UHF Trial %d-Run %d',n_trials,n_run));
if n_carriers==2
title(sprintf('Doppler spectrum 2.6GHz Carrier 1 Trial %d-Run %d',n_trials,n_run));
end
pcolor(T,F,10*log10( PDD_totala(:,:,itx,irx)));
shading flat
bara=colorbar;
ylabel(bara, 'dBm')
ylabel('F=f-ftx [Hz]')
xlabel('time [s]')
end
end
if(n_carriers==2)
figure(12)
for itx=1:p(2).nant_tx
for irx=1:p(2).nant_rx
T=1:1:block-1;
F=-(NFRAMES*num_symbols_frame/2-1)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560:15.36E6/(NFRAMES*num_symbols_frame/2)/2560:(NFRAMES*num_symbols_frame/2)*15.36E6/(2*NFRAMES*num_symbols_frame/2)/2560;
subplot(p(2).nant_tx,p(2).nant_rx,(itx-1)*p(2).nant_rx + irx);
title(sprintf('Doppler spectrum 2.6GHz Carrier 2 Trial %d-Run %d',n_trials,n_run));
pcolor(T,F,10*log10( PDD_totalb(:,:,itx,irx)));
shading flat
barb=colorbar;
ylabel(barb, 'dBm')
ylabel('F=f-ftx [Hz]')
xlabel('time [s]')
end
end
end
end
%%
%% Mean Delay
Pma=zeros((block-1)*NFRAMES,1,p(1).nant_tx,p(1).nant_rx);% zeroth-order moment
Pma1=zeros((block-1)*NFRAMES,1,p(1).nant_tx,p(1).nant_rx);
atau=linspace(0,p(1).useful_carriers/4/4.5E6,p(1).useful_carriers/4);
if(n_carriers==2)
atau=linspace(0,p(1).useful_carriers/4/18E6,p(1).useful_carriers/4);
end
for i=1:p(1).useful_carriers/4
Pma(:,1,:,:)=Pma(:,1,:,:)+PDP_totala(:,i,:,:);
Pma1(:,1,:,:)=Pma1(:,1,:,:)+atau(i)*PDP_totala(:,i,:,:);
end
mean_delay_a=Pma1./Pma;% mean delay: first-order moment
if(n_carriers==2)
Pmb=zeros((block-1)*NFRAMES,1,p(2).nant_tx,p(2).nant_rx);
Pmb1=zeros((block-1)*NFRAMES,1,p(2).nant_tx,p(2).nant_rx);
btau=linspace(0,p(2).useful_carriers/4/9E6,p(2).useful_carriers/4);
for i=1:p(2).useful_carriers/4
Pmb(:,1,:,:)=Pmb(:,1,:,:)+PDP_totalb(:,i,:,:);
Pmb1(:,1,:,:)=Pmb1(:,1,:,:)+btau(i)*PDP_totalb(:,i,:,:);
end
mean_delay_b=Pmb1./Pmb;
end
figure(13)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
title(sprintf('Mean Delay UHF Trial %d-Run %d',n_trials,n_run));
if n_carriers==2
title(sprintf('Mean Delay 2.6GHz Carrier 1 Trial %d-Run %d',n_trials,n_run));
end
plot(mean_delay_a(:,:,itx,irx));
ylabel('delay [s]')
xlabel('time [s]')
end
end
if (n_carriers==2)
figure(14)
for itx=1:p(2).nant_tx
for irx=1:p(2).nant_rx
subplot(p(2).nant_tx,p(2).nant_rx,(itx-1)*p(2).nant_rx + irx);
title(sprintf('Mean Delay 2.6GHz Carrier 2 Trial %d-Run %d',n_trials,n_run));
plot(mean_delay_b(:,:,itx,irx));
ylabel('delay [s]')
xlabel('time [s]')
end
end
end
%% Mean Doppler Shift
PDma=zeros(block-1,1,p(1).nant_tx,p(1).nant_rx);
PDma1=zeros(block-1,1,p(1).nant_tx,p(1).nant_rx);
theta=linspace(-NFRAMES*num_symbols_frame/2/2,NFRAMES*num_symbols_frame/2/2,NFRAMES*num_symbols_frame/2);
for j=1:block-1
for i=1:NFRAMES*num_symbols_frame/2
PDma(j,1,:,:)=PDma(j,1,:,:)+PDD_totala(i,j,:,:);
PDma1(j,1,:,:)=PDma1(j,1,:,:)+theta(i)*PDD_totala(i,j,:,:);
end
end
mean_doppler_shift_a=PDma1./PDma; % mean doppler shift
if(n_carriers==2)
PDmb=zeros(block-1,1,p(2).nant_tx,p(2).nant_rx);
PDmb1=zeros(block-1,1,p(2).nant_tx,p(2).nant_rx);
for j=1:block-1
for i=1:NFRAMES*num_symbols_frame/2
PDmb(j,1,:,:)=PDmb(j,1,:,:)+PDD_totalb(i,j,:,:);
PDmb1(j,1,:,:)=PDmb1(j,1,:,:)+theta(i)*PDD_totalb(i,j,:,:);
end
end
mean_doppler_shift_b=PDmb1./PDmb;
end
figure(15)
for itx=1:p(1).nant_tx
for irx=1:p(1).nant_rx
subplot(p(1).nant_tx,p(1).nant_rx,(itx-1)*p(1).nant_rx + irx);
title(sprintf('Mean Doppler shift UHF Trial %d-Run %d',n_trials,n_run));
if n_carriers==2
title(sprintf('Mean Doppler shift 2.6GHz Carrier 1 Trial %d-Run %d',n_trials,n_run));
end
plot(mean_doppler_shift_a(:,:,itx,irx));
ylabel('f-ftx [Hz]')
xlabel('time [s]')
end
end
if (n_carriers==2)
figure(16)
for itx=1:p(2).nant_tx
for irx=1:p(2).nant_rx
subplot(p(2).nant_tx,p(2).nant_rx,(itx-1)*p(2).nant_rx + irx);
title(sprintf('Mean Doppler shift 2.6GHz Carrier 2 Trial %d-Run %d',n_trials,n_run));
plot(mean_doppler_shift_b(:,:,itx,irx));
ylabel('f-ftx [Hz]')
xlabel('time [s]')
end
end
end
%%
figure(17)
plot(doppler_freq_of_max_a);
title(sprintf('Main Doppler peak for UHF Trial %d-Run%d',n_trials,n_run));
if n_carriers==2
title(sprintf('Main Doppler peak for 2.6GHz Carrier 1 Trial %d-Run%d',n_trials,n_run));
end
xlabel('time [s]');
ylabel('f-ftx [Hz]');
if(n_carriers==2)
figure(18)
plot(doppler_freq_of_max_b);
title(sprintf('Main Doppler peak for 2.6GHz Carrier 2 Trial %d-Run%d',n_trials,n_run));
xlabel('time [s]');
ylabel('f-ftx [Hz]');
end
%%
if(n_carriers==1)
if(n_trials==1)
subplot(2,2,1);
title('variation of the mean fdop');
plot(fm_total);
xlabel('time [s]');
ylabel('fdop [Hz]');
subplot(2,2,2);
title('variation of the TGV speed');
plot(TGVspeed_total);
xlabel('time [s]');
ylabel('TGV speed [km/h]');
subplot(2,2,3);
title('variation of the mean frequency offset');
plot(freqOffset_total);
xlabel('time [s]');
ylabel('frequency Offset [Hz]');
end
end
fclose(fid);
%% save processed data
if(n_carriers==1)
save(fullfile(destdir,[filename '.mat']),'PDP_totala','PDD_totala','mean_delay_a','mean_doppler_shift_a','doppler_freq_of_max_a','delay_doppler_profile_beforea','delay_doppler_profile_duringa','delay_doppler_profile_aftera','noise1','H_power1');
end
if(n_carriers==2)
save(fullfile(destdir,[filename '.mat']),'PDP_totala','PDD_totala','mean_delay_a','mean_doppler_shift_a','doppler_freq_of_max_a','delay_doppler_profile_beforea','delay_doppler_profile_duringa','delay_doppler_profile_aftera','PDP_totalb','PDD_totalb','mean_delay_b','mean_doppler_shift_b','doppler_freq_of_max_b','delay_doppler_profile_beforeb','delay_doppler_profile_duringb','delay_doppler_profile_afterb','noise1','H_power1','noise2','H_power2');
end
\ No newline at end of file
function [MSEsubsp,minind]=findminimumd(sigma2, nM, nD, Sb)
% Changes:
% Author: Thomas Zemen
% Copyright (c) by Forschungszentrum Telekommunikation Wien (ftw.)
for i=1:nD
MSEsubsp(i)=i/nM*sigma2+1/nM*sum(Sb(i+1:nD));
end
[MSEmin,minind]=min(MSEsubsp);
\ No newline at end of file
%addpath('../../../openair1/SIMULATION/LTE_PHY/')
%addpath('../../../openair1/PHY/LTE_ESTIMATION/')
%addpath('../../../openair1/PHY/LTE_REFSIG/')
%addpath('../../../targets/ARCH/EXMIMO/USERSPACE/OCTAVE')
nb_rb = 100; %this can be 25, 50, or 100
num_carriers = 2048/100*nb_rb;
num_zeros = num_carriers-(12*nb_rb+1);
prefix_length = num_carriers/4; %this is extended CP
num_symbols_frame = 120;
preamble_length = 120;
% this generates one LTE frame (10ms) full of OFDM modulated random QPSK
% symbols
[s,f] = OFDM_TX_FRAME(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length);
% scale to conserve energy (Matlabs IFFT does not scale)
s=s*sqrt(num_carriers);
% load the LTE sync sequence, upsample it to the right frequency and insert
% it in the first symbol of the frame
primary_synch;
pss0_up = interp(primary_synch0_time,num_carriers/128);
pss0_up_cp = [pss0_up(num_carriers-prefix_length+1:end) pss0_up];
s(1:num_carriers+prefix_length) = pss0_up_cp;
plot(abs(s))
% save for later use (channel estimation and transmission with the SMBV)
save(sprintf('ofdm_pilots_sync_%d.mat',num_carriers),'-v7','s','f','num_carriers','num_zeros','prefix_length','num_symbols_frame','preamble_length');
mat2wv(s, sprintf('ofdm_pilots_sync_%d.wv',num_carriers), 30.72e6/2048*num_carriers, 1)
%% this script generates the signals for the CORRIDOR channel sounding campaing
addpath('../../../openair1/SIMULATION/LTE_PHY/')
%addpath('../../../openair1/PHY/LTE_ESTIMATION/')
addpath('../../../openair1/PHY/LTE_REFSIG/')
%addpath('../../../targets/ARCH/EXMIMO/USERSPACE/OCTAVE')
rand('seed',42); %make sure seed random numbers are alwyas the same
% load the LTE sync sequence
primary_synch;
nant = 4;
%% this generates one LTE frame (10ms) full of OFDM modulated random QPSK symbols
%% 20MHz carrier
nb_rb = 100; %this can be 25, 50, or 100
num_carriers = 2048/100*nb_rb;
num_zeros = num_carriers-(12*nb_rb+1);
prefix_length = num_carriers/4; %this is extended CP
num_symbols_frame = 120;
preamble_length = 120;
[s1,f1] = OFDM_TX_FRAME_MIMO(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length,nant);
% scale to conserve energy (Matlabs IFFT does not scale)
s1=s1*sqrt(num_carriers);
% upsample PSS to the right frequency and insert it in the first symbol of the frame
pss0_up = interp(primary_synch0_time,num_carriers/128);
pss0_up_cp = [pss0_up(num_carriers-prefix_length+1:end) pss0_up];
s1(:,1:num_carriers+prefix_length) = repmat(pss0_up_cp,nant,1);
%% 10MHz carrier
nb_rb = 50; %this can be 25, 50, or 100
num_carriers = 2048/100*nb_rb;
num_zeros = num_carriers-(12*nb_rb+1);
prefix_length = num_carriers/4; %this is extended CP
[s2,f2] = OFDM_TX_FRAME_MIMO(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length,nant);
% scale to conserve energy (Matlabs IFFT does not scale)
s2=s2*sqrt(num_carriers);
% upsample PSS to the right frequency and insert it in the first symbol of the frame
pss0_up = interp(primary_synch0_time,num_carriers/128);
pss0_up_cp = [pss0_up(num_carriers-prefix_length+1:end) pss0_up];
s2(:,1:num_carriers+prefix_length) = repmat(pss0_up_cp,nant,1);
%% 5MHz carrier
nb_rb = 25; %this can be 25, 50, or 100
num_carriers = 2048/100*nb_rb;
num_zeros = num_carriers-(12*nb_rb+1);
prefix_length = num_carriers/4; %this is extended CP
[s3,f3] = OFDM_TX_FRAME_MIMO(num_carriers,num_zeros,prefix_length,num_symbols_frame,preamble_length,nant);
% scale to conserve energy (Matlabs IFFT does not scale)
s3=s3*sqrt(num_carriers);
% upsample PSS to the right frequency and insert it in the first symbol of the frame
pss0_up = interp(primary_synch0_time,num_carriers/128);
pss0_up_cp = [pss0_up(num_carriers-prefix_length+1:end) pss0_up];
s3(:,1:num_carriers+prefix_length) = repmat(pss0_up_cp,nant,1);
%% combine the 10 and 20 MHz carriers
f1_shift = -5e6;
f2_shift = 10e6;
sample_rate = 30.72e6*2;
s = zeros(nant,sample_rate/100);
for a=1:nant
s1_up = interp(s1(a,:),2);
s1_shift = s1_up .* exp(2*1i*pi*f1_shift*(0:length(s1_up)-1)/sample_rate);
s2_up = interp(s2(a,:),4);
s2_shift = s2_up .* exp(2*1i*pi*f2_shift*(0:length(s2_up)-1)/sample_rate);
s(a,:) = s1_shift + s2_shift/sqrt(2);
end
%%
figure(1)
hold off
plot(linspace(-sample_rate/2,sample_rate/2,length(s)),20*log10(abs(fftshift(fft(s,[],2)))))
%% save for later use (channel estimation and transmission with the SMBV)
save('ofdm_pilots_sync_30MHz.mat','-v7','s1','s2','s3','f1','f2','f3','num_carriers','num_zeros','prefix_length','num_symbols_frame','preamble_length');
s_all = sum(s,1);
s_all(1:5120) = s(1,1:5120);
mat2wv(s_all, 'ofdm_pilots_sync_30MHz.wv', sample_rate, 1);
s_all = sum(s3,1);
s_all(1:640) = s3(1,1:640);
mat2wv(s_all, 'ofdm_pilots_sync_5MHz.wv', 7.68e6, 1);
function p_dpss = init_dpss(p,tau_max)
%% init DPSS
% channel parameters
theta_max = tau_max * p.frame_length*100 / p.useful_carriers; % TS = 1/(p.frame_length*100)
N_sub = floor(theta_max*p.useful_carriers)+1; % number of approximate subspace dimensions
%N_sub = 1;
% DPSS generation
nDmax=min(ceil(N_sub*1.25),p.useful_carriers);
[Udps,Sdps]=dpss(p.useful_carriers,theta_max*p.useful_carriers/2,nDmax);
% the fftshift is needed dependig in which format the frequency response is used
%V=fftshift(diag(exp(-2*pi*1i*theta_max/2*(0:(p.nN-1))))*Udps);
p_dpss.V=diag(exp(-2*pi*1i*theta_max/2*(0:( p.useful_carriers-1))))*Udps;
% compute the optimal subspace dimension based on the expected SNR
p_dpss.SNR=0:30;
for SNR=p_dpss.SNR
% calculate optimum subspace dimension as bias variance tradeoff for a
% given noise variance (assuming pilot symbols with energy one)
sigma2 = 10.^(-SNR/10);
[~,p_dpss.Dopt(SNR+1)]=findminimumd(sigma2, p.useful_carriers,nDmax,1/theta_max*Sdps);
end
% % precompute the filters for DPSS estimation for the required number of
% % basis function
% Dvec = min(Dopt):max(Dopt);
% for ind=1:length(Dvec)
% Dind = Dvec(ind);
%
% for itx=1:p.nant_tx
% % f_start and t_start indicate the start of the pilots in time
% % and frequency according to the specifications (see .doc file).
% % t_start has to be >=2, since the first symbol is the PSS.
% f_start = mod(itx-1,2)*2+1;
% t_start = floor((itx-1)/2)+1;
%
% % filter for DPSS channel estimation
% M = conj(diag(squeeze(transmit_f(itx,t_start,f_start:4:end))) * V(f_start:4:end,1:Dind));
% Mpinv = (M'*M)\M.';
% % to compute the basis coefficients do
% % psi = Mpinv*(symb0(2:6:end));
%
% end
%
% end
function p = init_params(nb_rb,nant_rx,nant_tx)
global symbols_per_slot slots_per_frame;
num_symbols_frame = symbols_per_slot*slots_per_frame;
p.nb_rb = nb_rb;
p.num_carriers = 2048/(100/p.nb_rb);
p.num_zeros = p.num_carriers-(12*p.nb_rb+1);
p.useful_carriers = p.num_carriers-p.num_zeros-1;
p.prefix_length = p.num_carriers/4; %this is extended CP
p.ofdm_symbol_length = p.num_carriers + p.prefix_length;
p.samples_slot = p.ofdm_symbol_length*symbols_per_slot;
p.frame_length = p.ofdm_symbol_length*num_symbols_frame;
p.nant_rx=nant_rx;
p.nant_tx=nant_tx;
clear all
close all
top_dir = 'E:\EMOS\CORRIDOR\trials2 train'; % needs to be updated according to your computer
d1 = dir(fullfile(top_dir,'UHF','*.log'));
d2 = dir(fullfile(top_dir,'2.6GHz','*.log'));
load exmimo2_39_0.8G.mat
G0 = permute(G0,[2 1 3]);
G0_800_interp = interp1(ALL_gain,G0,0:30);
G0_800_interp = permute(G0_800_interp,[2 1 3]);
NF0 = permute(NF0,[2 1 3]);
NF0_800_interp = interp1(ALL_gain,NF0,0:30);
NF0_800_interp = permute(NF0_800_interp,[2 1 3]);
load exmimo2_39_2.6G_v2.mat
G0 = permute(G0,[2 1 3]);
G0_2600_interp = interp1(ALL_gain,G0,0:30);
G0_2600_interp = permute(G0_2600_interp,[2 1 3]);
NF0 = permute(NF0,[2 1 3]);
NF0_2600_interp = interp1(ALL_gain,NF0,0:30);
NF0_2600_interp = permute(NF0_2600_interp,[2 1 3]);
%%
start_time = [1.400489088000000e+09 1.400493112000000e+09 1.400499696000000e+09 1.400506864000000e+09];
for idx=1:length(d1)
data1{idx}=csvread(fullfile(top_dir,'UHF',d1(idx).name),1,0);
data2{idx}=csvread(fullfile(top_dir,'2.6GHz',d2(idx).name),1,0);
frame_start1(idx) = ceil(data1{idx}(find(data1{idx}(:,1)>start_time(idx),1,'first'),3)/92160000)*100;
frame_start2(idx) = ceil(data2{idx}(find(data2{idx}(:,1)>start_time(idx),1,'first'),3)/368640000)*100;
%% find the first dataset with valid GPS signal and throw away everything before
idx1_start = find(data1{idx}(:,6)>0,1,'first');
data1{idx}(1:idx1_start-1,:)=[];
idx2_start = find(data2{idx}(:,6)>0,1,'first');
data2{idx}(1:idx2_start-1,:)=[];
rtime1 = data1{idx}(:,1) - data1{idx}(1,1);
rtime2 = data2{idx}(:,1) - data2{idx}(1,1);
%% compute the noise level based on the AGC values and the calibrated noise figures
data1{idx}(data1{idx}(:,16)==0,16) = 1;
data2{idx}(data2{idx}(:,16)==0,16) = 1;
data2{idx}(data2{idx}(:,22)==0,22) = 1;
data2{idx}(data2{idx}(:,28)==0,28) = 1;
NF1=[NF0_800_interp(sub2ind(size(NF0_800_interp),data1{idx}(:,16),data1{idx}(:,15)+1,ones(length(data1{idx}),1)))];
NF2=[NF0_2600_interp(sub2ind(size(NF0_2600_interp),data2{idx}(:,16),data2{idx}(:,15)+1,ones(length(data2{idx}),1))) ...
NF0_2600_interp(sub2ind(size(NF0_2600_interp),data2{idx}(:,22),data2{idx}(:,21)+1,ones(length(data2{idx}),1))) ...
NF0_2600_interp(sub2ind(size(NF0_2600_interp),data2{idx}(:,28),data2{idx}(:,27)+1,ones(length(data2{idx}),1)))];
G1=[G0_800_interp(sub2ind(size(G0_800_interp),data1{idx}(:,16),data1{idx}(:,15)+1,ones(length(data1{idx}),1)))];
G2=[G0_2600_interp(sub2ind(size(G0_2600_interp),data2{idx}(:,16),data2{idx}(:,15)+1,ones(length(data2{idx}),1))) ...
G0_2600_interp(sub2ind(size(G0_2600_interp),data2{idx}(:,22),data2{idx}(:,21)+1,ones(length(data2{idx}),1))) ...
G0_2600_interp(sub2ind(size(G0_2600_interp),data2{idx}(:,28),data2{idx}(:,27)+1,ones(length(data2{idx}),1)))];
% digital noise power in dB
N1 = -174 + 10*log10(7.68e6) + NF1 + G1;
N2 = repmat(-174+10*log10([30.72e6 30.72e6 15.36e6]),length(NF2),1) + NF2 + G2;
%% plot gps coordinates
h=figure(idx*10+1);
hold off
plot(data1{idx}(1:100:end,7),data1{idx}(1:100:end,8),'rx')
hold on
plot(data2{idx}(1:100:end,7),data2{idx}(1:100:end,8),'bx')
xlabel('lat [deg]');
ylabel('lon [deg]')
legend('UHF','2.6GHz')
title(sprintf('Run %d',idx));
saveas(h,sprintf('figures/gps_trace_run%d.eps',idx));
%% plot RSSI
% TODO: convert time (in unix epoch) into something more meaninful
h=figure(idx*10+2);
hold off
plot(rtime1,smooth(data1{idx}(:,13),100),'r')
hold on
plot(rtime2,smooth(data2{idx}(:,13),100),'b')
plot(rtime2,smooth(data2{idx}(:,19),100),'c')
plot(rtime2,smooth(data2{idx}(:,25),100),'m')
legend('UHF','2.6GHz card 1','2.6GHz card 2','2.6GHz card 3');
xlabel('time [seconds]')
ylabel('RSSI [dBm]')
title(sprintf('Run %d',idx));
saveas(h,sprintf('figures/rssi_vs_time_run%d.eps',idx),'epsc2');
%% plot NF
% TODO: convert time (in unix epoch) into something more meaninful
h=figure(idx*10+5);
hold off
plot(rtime1,smooth(N1,100),'r')
hold on
plot(rtime2,smooth(N2(:,1),100),'b')
plot(rtime2,smooth(N2(:,2),100),'c')
plot(rtime2,smooth(N2(:,3),100),'m')
legend('UHF','2.6GHz card 1','2.6GHz card 2','2.6GHz card 3');
xlabel('time [seconds]')
ylabel('Noise [dB]')
title(sprintf('Run %d',idx));
saveas(h,sprintf('figures/nf_vs_time_run%d.eps',idx),'epsc2');
%% measured distance (km)
% We can get the distance between the base station and the RX antenna from the GPS coordinates
distances1=zeros(size(data1{idx},1),1);
distances2=zeros(size(data2{idx},1),1);
for i=1:size(data1{idx},1)
distances1(i)=Dist_Calc_from_GPS(data1{idx}(i,7),data1{idx}(i,8),48.25073056,1.55481944);
end
for i=1:size(data2{idx},1)
distances2(i)=Dist_Calc_from_GPS(data2{idx}(i,7),data2{idx}(i,8),48.25073056,1.55481944);
end
figure (10*idx+6)
subplot(1,2,1)
plot(rtime1,distances1,'r',rtime1,smooth(data1{idx}(:,13),100),'b')
title(sprintf('Run %d with the measured distance : UHF',idx));
xlabel('time [s]');
legend('distance [km]','RSSI [dBm]');
subplot(1,2,2)
plot(rtime2,distances2,'r',rtime2,smooth(data2{idx}(:,13),100),'b')
title(sprintf('Run %d with the measured distance : 2.6GHz ',idx));
xlabel('time [s]');
legend('distance [km]','RSSI [dBm]');
%% estimated distance under the assumption of a constant speed
% We assume that the TGV speed is constant and then we find the distances to the base station with the time vector
TGV_speed=82.5;%constant TGV speed in m/s
[RSSI_max1,I_RSSI_max1]=max(data1{idx}(:,13));%we find the index corresponding to the maximum of RSSI
time01=rtime1(I_RSSI_max1)*ones(length(rtime1),1);%time corresponding to the maximim of RSSI
new_distances1=(TGV_speed*abs(rtime1-time01))/1000+min(distances1)*ones(length(rtime1),1);% new distance in km. is is minimum when the RSSI is maximum
[RSSI_max2,I_RSSI_max2]=max(data2{idx}(:,13));
time02=rtime2(I_RSSI_max2)*ones(length(rtime2),1);
new_distances2=(TGV_speed*abs(rtime2-time02))/1000+min(distances2)*ones(length(rtime2),1);% distance in km
if (idx==2)%For Run 2, there is an anomalous peak for the RSSI at the end. Here we ignore it
[RSSI_max2,I_RSSI_max2]=max(data2{idx}(1:32900,13));
time02=rtime2(I_RSSI_max2)*ones(length(rtime2),1);
new_distances2=(TGV_speed*abs(rtime2-time02))/1000+min(distances2)*ones(length(rtime2),1);% distance in km
end
figure (10*idx+7)
subplot(1,2,1)
plot(rtime1,new_distances1,'r',rtime1,smooth(data1{idx}(:,13),100),'b')
title(sprintf('Run %d with the estimated distance : UHF',idx));
xlabel('time [s]');
legend('distance [km]','RSSI [dBm]');
subplot(1,2,2)
plot(rtime2,new_distances2,'r',rtime2,smooth(data2{idx}(:,13),100),'b')
title(sprintf('Run %d with the estimated distance : 2.6GHz ',idx));
xlabel('time [s]');
legend('distance [km]','RSSI [dBm]');
%% rssi(dBm) versus distance (log scale)
% We will plot the rssi versus the distance. We separate the data
% before and after the passing of the train for run 3 and run 4 because
% all the antennas are poiting at the same direction
% we heuristically determine a starting point and a ending point for the linear fitting
if idx==1
distance_break1_start=0.4;%in km
distance_break1_end=7.5;
distance_break2_start=0.4;
distance_break2_end=7.5;
end
if idx==2
distance_break1_start=0.8;%in km
distance_break1_end=4.5;
distance_break2_start=0.8;
distance_break2_end=4.5;
end
if idx==3
distance_before_break1_start=5;%in km
distance_before_break1_end=23;
distance_before_break2_start=5;
distance_before_break2_end=23;
distance_after_break1_start=2.274;%in km
distance_after_break1_end=5.376;
distance_after_break2_start=5.996;
distance_after_break2_end=7.613;
end
if idx==4
distance_before_break1_start=0.1176;%in km
distance_before_break1_end=3.489;
distance_before_break2_start=0.1344;
distance_before_break2_end=3.78;
distance_after_break1_start=0.5;%in km
distance_after_break1_end=8;
distance_after_break2_start=0.5;
distance_after_break2_end=8;
end
if idx==1 || idx==2
% indexes of the starting and ending points
index_break1_start=1;
index_break2_start=1;
index_break1_end=1;
index_break2_end=1;
end
if idx==3 || idx==4
% indexes of the starting and ending points with the data before the passing of the
% train
index_break1_before_start=1;
index_break2_before_start=1;
index_break1_before_end=1;
index_break2_before_end=1;
% indexes of the starting and ending points with the data after the passing of the
% train
index_break1_after_start=I_RSSI_max1;
index_break2_after_start=I_RSSI_max2;
index_break1_after_end=I_RSSI_max1;
index_break2_after_end=I_RSSI_max2;
end
if idx==1 || idx==2
%starting points
while (index_break1_start<length(new_distances1)) && (new_distances1(index_break1_start)>distance_break1_start)
index_break1_start=index_break1_start+1;
end
while (index_break2_start<length(new_distances2)) && (new_distances2(index_break2_start)>distance_break2_start)
index_break2_start=index_break2_start+1;
end
%ending points
while (index_break1_end<length(new_distances1)) && (new_distances1(index_break1_end)>distance_break1_end)
index_break1_end=index_break1_end+1;
end
while (index_break2_end<length(new_distances2)) && (new_distances2(index_break2_end)>distance_break2_end)
index_break2_end=index_break2_end+1;
end
end
if idx==3 || idx==4
%starting points
while (index_break1_before_start<length(new_distances1)) && (new_distances1(index_break1_before_start)>distance_before_break1_start)
index_break1_before_start=index_break1_before_start+1;
end
while (index_break2_before_start<length(new_distances2)) && (new_distances2(index_break2_before_start)>distance_before_break2_start)
index_break2_before_start=index_break2_before_start+1;
end
%ending points
while (index_break1_before_end<length(new_distances1)) && (new_distances1(index_break1_before_end)>distance_before_break1_end)
index_break1_before_end=index_break1_before_end+1;
end
while (index_break2_before_end<length(new_distances2)) && (new_distances2(index_break2_before_end)>distance_before_break2_end)
index_break2_before_end=index_break2_before_end+1;
end
%starting points
while (index_break1_after_start<length(new_distances1)) && (new_distances1(index_break1_after_start)<distance_after_break1_start)
index_break1_after_start=index_break1_after_start+1;
end
while (index_break2_after_start<length(new_distances2)) && (new_distances2(index_break2_after_start)<distance_after_break2_start)
index_break2_after_start=index_break2_after_start+1;
end
%ending points
while (index_break1_after_end<length(new_distances1)) && (new_distances1(index_break1_after_end)<distance_after_break1_end)
index_break1_after_end=index_break1_after_end+1;
end
while (index_break2_after_end<length(new_distances2)) && (new_distances2(index_break2_after_end)<distance_after_break2_end)
index_break2_after_end=index_break2_after_end+1;
end
end
if idx==1 || idx ==2
h=figure(idx*10+3);
hold off
linearCoef1 = polyfit(10*log10(new_distances1(index_break1_end:index_break1_start)),data1{idx}(index_break1_end:index_break1_start,13),1);
linearFit1 = polyval(linearCoef1,10*log10(new_distances1(index_break1_end:index_break1_start)));
semilogx(new_distances1(1:I_RSSI_max1),data1{idx}(1:I_RSSI_max1,13),'rx',new_distances1(index_break1_end:index_break1_start),linearFit1,'r-')
display(sprintf('Run %d :slope UHF : %f',idx,linearCoef1(1)))
hold on
linearCoef2 = polyfit(10*log10(new_distances2(index_break2_end:index_break2_start)),data2{idx}(index_break2_end:index_break2_start,13),1);
linearFit2 = polyval(linearCoef2,10*log10(new_distances2(index_break2_end:index_break2_start)));
semilogx(new_distances2(1:I_RSSI_max2),data2{idx}(1:I_RSSI_max2,13),'bx',new_distances2(index_break2_end:index_break2_start),linearFit2,'b-')
display(sprintf('Run %d :slope 2.6GHz : %f',idx,linearCoef2(1)))
title(sprintf('Run %d',idx))
legend('UHF','UHF linear fit','2.6GHz','2.6GHz linear fit');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
saveas(h,sprintf('figures/rssi_vs_dist_run%d.eps',idx),'epsc2');
% Zoom on the linear fitting
h=figure(idx*10+4);
hold off
linearCoef1 = polyfit(10*log10(new_distances1(index_break1_end:index_break1_start)),data1{idx}(index_break1_end:index_break1_start,13),1);
linearFit1 = polyval(linearCoef1,10*log10(new_distances1(index_break1_end:index_break1_start)));
semilogx(new_distances1(index_break1_end:index_break1_start),data1{idx}(index_break1_end:index_break1_start,13),'rx',new_distances1(index_break1_end:index_break1_start),linearFit1,'r-')
%display(sprintf('Run %d :slope UHF : %f',idx,linearCoef1(1)))
hold on
linearCoef2 = polyfit(10*log10(new_distances2(index_break2_end:index_break2_start)),data2{idx}(index_break2_end:index_break2_start,13),1);
linearFit2 = polyval(linearCoef2,10*log10(new_distances2(index_break2_end:index_break2_start)));
semilogx(new_distances2(index_break2_end:index_break2_start),data2{idx}(index_break2_end:index_break2_start,13),'bx',new_distances2(index_break2_end:index_break2_start),linearFit2,'b-')
%display(sprintf('Run %d :slope 2.6GHz : %f',idx,linearCoef2(1)))
title(sprintf('Run %d',idx))
legend('UHF','UHF linear fit','2.6GHz','2.6GHz linear fit');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
saveas(h,sprintf('figures/rssi_vs_dist_zoom_run%d.eps',idx),'epsc2');
end
if idx==3 || idx==4
h=figure(idx*10+3);
subplot(2,1,1)
hold off
linearCoef1_before = polyfit(10*log10(new_distances1(index_break1_before_end:index_break1_before_start)),data1{idx}(index_break1_before_end:index_break1_before_start,13),1);
linearFit1_before = polyval(linearCoef1_before,10*log10(new_distances1(index_break1_before_end:index_break1_before_start)));
semilogx(new_distances1(1:I_RSSI_max1),data1{idx}(1:I_RSSI_max1,13),'rx',new_distances1(index_break1_before_end:index_break1_before_start),linearFit1_before,'r-')
display(sprintf('Run %d :slope UHF before: %f',idx,linearCoef1_before(1)))
hold on
linearCoef2_before = polyfit(10*log10(new_distances2(index_break2_before_end:index_break2_before_start)),data2{idx}(index_break2_before_end:index_break2_before_start,13),1);
linearFit2_before = polyval(linearCoef2_before,10*log10(new_distances2(index_break2_before_end:index_break2_before_start)));
semilogx(new_distances2(1:I_RSSI_max2),data2{idx}(1:I_RSSI_max2,13),'bx',new_distances2(index_break2_before_end:index_break2_before_start),linearFit2_before,'b-')
display(sprintf('Run %d :slope 2.6GHz before: %f',idx,linearCoef2_before(1)))
title(sprintf('Run %d: With the data before the passing of the train',idx))
legend('UHF','UHF:linear fitting','2.6GHz card 1','2.6GHz card 1:linear fitting');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
subplot(2,1,2)
hold off
linearCoef1_after = polyfit(10*log10(new_distances1(index_break1_after_start:index_break1_after_end)),data1{idx}(index_break1_after_start:index_break1_after_end,13),1);
linearFit1_after = polyval(linearCoef1_after,10*log10(new_distances1(index_break1_after_start:index_break1_after_end)));
semilogx(new_distances1(I_RSSI_max1:end),data1{idx}(I_RSSI_max1:end,13),'rx',new_distances1(index_break1_after_start:index_break1_after_end),linearFit1_after,'r-')
display(sprintf('Run %d :slope UHF after: %f',idx,linearCoef1_after(1)))
hold on
linearCoef2_after = polyfit(10*log10(new_distances2(index_break2_after_start:index_break2_after_end)),data2{idx}(index_break2_after_start:index_break2_after_end,13),1);
linearFit2_after = polyval(linearCoef2_after,10*log10(new_distances2(index_break2_after_start:index_break2_after_end)));
semilogx(new_distances2(I_RSSI_max2:end),data2{idx}(I_RSSI_max2:end,13),'bx',new_distances2(index_break2_after_start:index_break2_after_end),linearFit2_after,'b-')
display(sprintf('Run %d :slope 2.6GHz after: %f',idx,linearCoef2_after(1)))
title(sprintf('Run %d: With the data after the passing of the train',idx))
legend('UHF','UHF:linear fitting','2.6GHz card 1','2.6GHz card 1:linear fitting');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
saveas(h,sprintf('figures/rssi_vs_dist_run%d.eps',idx),'epsc2');
% Zoom on the linear fitting
h=figure(idx*10+4);
subplot(2,1,1)
hold off
linearCoef1_before = polyfit(10*log10(new_distances1(index_break1_before_end:index_break1_before_start)),data1{idx}(index_break1_before_end:index_break1_before_start,13),1);
linearFit1_before = polyval(linearCoef1_before,10*log10(new_distances1(index_break1_before_end:index_break1_before_start)));
semilogx(new_distances1(index_break1_before_end:index_break1_before_start),data1{idx}(index_break1_before_end:index_break1_before_start,13),'rx',new_distances1(index_break1_before_end:index_break1_before_start),linearFit1_before,'r-')
%display(sprintf('Run %d :slope UHF before: %f',idx,linearCoef1_before(1)))
hold on
linearCoef2_before = polyfit(10*log10(new_distances2(index_break2_before_end:index_break2_before_start)),data2{idx}(index_break2_before_end:index_break2_before_start,13),1);
linearFit2_before = polyval(linearCoef2_before,10*log10(new_distances2(index_break2_before_end:index_break2_before_start)));
semilogx(new_distances2(index_break2_before_end:index_break2_before_start),data2{idx}(index_break2_before_end:index_break2_before_start,13),'bx',new_distances2(index_break2_before_end:index_break2_before_start),linearFit2_before,'b-')
%display(sprintf('Run %d :slope 2.6GHz before: %f',idx,linearCoef2_before(1)))
title(sprintf('Run %d: With the data before the passing of the train',idx))
legend('UHF','UHF:linear fitting','2.6GHz card 1','2.6GHz card 1:linear fitting');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
subplot(2,1,2)
hold off
linearCoef1_after = polyfit(10*log10(new_distances1(index_break1_after_start:index_break1_after_end)),data1{idx}(index_break1_after_start:index_break1_after_end,13),1);
linearFit1_after = polyval(linearCoef1_after,10*log10(new_distances1(index_break1_after_start:index_break1_after_end)));
semilogx(new_distances1(index_break1_after_start:index_break1_after_end),data1{idx}(index_break1_after_start:index_break1_after_end,13),'rx',new_distances1(index_break1_after_start:index_break1_after_end),linearFit1_after,'r-')
%display(sprintf('Run %d :slope UHF after: %f',idx,linearCoef1_after(1)))
hold on
linearCoef2_after = polyfit(10*log10(new_distances2(index_break2_after_start:index_break2_after_end)),data2{idx}(index_break2_after_start:index_break2_after_end,13),1);
linearFit2_after = polyval(linearCoef2_after,10*log10(new_distances2(index_break2_after_start:index_break2_after_end)));
semilogx(new_distances2(index_break2_after_start:index_break2_after_end),data2{idx}(index_break2_after_start:index_break2_after_end,13),'bx',new_distances2(index_break2_after_start:index_break2_after_end),linearFit2_after,'b-')
%display(sprintf('Run %d :slope 2.6GHz after: %f',idx,linearCoef2_after(1)))
title(sprintf('Run %d: With the data after the passing of the train',idx))
legend('UHF','UHF:linear fitting','2.6GHz card 1','2.6GHz card 1:linear fitting');
xlabel('distance [km]')
ylabel('RSSI [dBm]')
saveas(h,sprintf('figures/rssi_vs_dist_zoom_run%d.eps',idx),'epsc2');
end
end
\ No newline at end of file
function [m,ind]=peaksfinder(corr,frame_length)
threshold=max(abs(corr))*0.75;
consecutivePos=[];
highCorrVal=find(abs(corr)>threshold);
num_peak=0;
k=1;
consecutivePos(1)=highCorrVal(1);
for i_high=1:size(highCorrVal,1)-1
if highCorrVal(i_high+1)-highCorrVal(i_high)==1
consecutivePos(k+1)=highCorrVal(i_high+1);
k=k+1;
else
num_peak=num_peak+1;
[m(num_peak),temp_idx]=max(abs(corr(min(consecutivePos):max(consecutivePos))));
ind(num_peak)=min(consecutivePos)-1+temp_idx;
consecutivePos=[];
consecutivePos(1)=highCorrVal(i_high+1);
k=1;
end
end
num_peak=num_peak+1;
[m(num_peak),temp_idx]=max(abs(corr(min(consecutivePos):max(consecutivePos))));
ind(num_peak)=min(consecutivePos)-1+temp_idx;
% a bigining for make code which clculate the best peak in respect to whole frlames
% corrMatrix=vec2mat(corr, frame_length);
% corrMatrix(end,:)=[];
% sumCorrMatrix=sum(abs(corrMatrix));
% [Value Pos]=max(sumCorrMatrix);
% this script sends the signals previously generated with generation_ca to
% the cards. We need in total 7 cards configured as follows
% card0 - card3: 20MHz, 1 channel each, s1, freq 2590 MHz
% card4 - card5: 10MHz, 2 channels each, s2, freq 2605 MHz
% card6: 5MHz, 4 channels, s3, freq 771.5 MHz
load('ofdm_pilots_sync_30MHz.mat');
addpath('../../../targets/ARCH/EXMIMO/USERSPACE/OCTAVE')
limeparms;
num_cards = oarf_get_num_detected_cards;
card_select = [1 1 1 1 1 1 1];
% common parameters
rf_local= rf_local*[1 1 1 1];
rf_rxdc = rf_rxdc*[1 1 1 1];
rf_vcocal=rf_vcocal_19G*[1 1 1 1];
tx_gain_all = [10 0 0 0;
10 0 0 0;
15 0 0 0;
12 0 0 0;
1 0 0 0;
0 10 10 0;
9 8 0 0];
%tx_gain_all = [11 0 0 0;
% 8 0 0 0;
% 14 0 0 0;
% 13 0 0 0;
% 5 3 0 0;
% 13 13 0 0;
% 0 0 0 0];
for card=0:min(3,num_cards-1)
disp(card)
% card 0-3: 20MHz
active_rf = [1 0 0 0];
autocal = [1 1 1 1];
resampling_factor = [0 0 0 0];
rf_mode = (RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF10+RXLPFNORM+RXLPFEN+RXLPF10+LNA1ON+LNAMax+RFBBNORM)*active_rf;
rf_mode = rf_mode+((DMAMODE_RX+DMAMODE_TX)*active_rf);
freq_rx = 2590e6*active_rf;
freq_tx = freq_rx;
tx_gain = tx_gain_all(card+1,:);
rx_gain = 0*active_rf; %1 1 1 1];
eNBflag = 0;
tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_TESTTX;
%tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_LSB;
if (card==0)
syncmode = SYNCMODE_MASTER;
else
syncmode = SYNCMODE_SLAVE;
end
%syncmode = SYNCMODE_FREE;
rffe_rxg_low = 31*active_rf; %[1 1 1 1];
rffe_rxg_final = 63*active_rf; %[1 1 1 1];
rffe_band = B19G_TDD*active_rf; %[1 1 1 1];
oarf_config_exmimo(card,freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNBflag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal,resampling_factor);
end
for card=5:min(5,num_cards-1)
disp(card)
% card 5: 10MHz
active_rf = [0 1 1 0];
autocal = [1 1 1 1];
resampling_factor = [1 1 1 1];
rf_mode = (RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF5+RXLPFNORM+RXLPFEN+RXLPF5+LNA1ON+LNAMax+RFBBNORM)*active_rf;
rf_mode = rf_mode+((DMAMODE_RX+DMAMODE_TX)*active_rf);
freq_rx = 2605e6*active_rf;
freq_tx = freq_rx;
tx_gain = tx_gain_all(card+1,:);
rx_gain = 0*active_rf; %1 1 1 1];
eNBflag = 0;
tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_TESTTX;
%tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_LSB;
syncmode = SYNCMODE_SLAVE;
%syncmode = SYNCMODE_FREE;
rffe_rxg_low = 31*active_rf; %[1 1 1 1];
rffe_rxg_final = 63*active_rf; %[1 1 1 1];
rffe_band = B19G_TDD*active_rf; %[1 1 1 1];
oarf_config_exmimo(card,freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNBflag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal,resampling_factor);
end
for card=6:min(6,num_cards-1)
disp(card)
% card 6: 10MHz
active_rf = [1 1 0 0];
autocal = [1 1 1 1];
resampling_factor = [1 1 1 1];
rf_mode = (RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF5+RXLPFNORM+RXLPFEN+RXLPF5+LNA1ON+LNAMax+RFBBNORM)*active_rf;
rf_mode = rf_mode+((DMAMODE_RX+DMAMODE_TX)*active_rf);
freq_rx = 2605e6*active_rf;
freq_tx = freq_rx;
tx_gain = tx_gain_all(card+1,:);
rx_gain = 0*active_rf; %1 1 1 1];
eNBflag = 0;
tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_TESTTX;
%tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_LSB;
syncmode = SYNCMODE_SLAVE;
%syncmode = SYNCMODE_FREE;
rffe_rxg_low = 31*active_rf; %[1 1 1 1];
rffe_rxg_final = 63*active_rf; %[1 1 1 1];
rffe_band = B19G_TDD*active_rf; %[1 1 1 1];
oarf_config_exmimo(card,freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNBflag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal,resampling_factor);
end
for card=4:min(4,num_cards-1)
disp(card)
% card 4: 5MHz
active_rf = [1 1 1 1];
autocal = [1 1 1 1];
resampling_factor = [2 2 2 2];
rf_mode = (RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM)*active_rf;
rf_mode = rf_mode+((DMAMODE_RX+DMAMODE_TX)*active_rf);
freq_rx = 771.5e6*active_rf;
freq_tx = freq_rx;
tx_gain = tx_gain_all(card+1,:);
rx_gain = 0*active_rf; %1 1 1 1];
eNBflag = 0;
tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_TESTTX;
%tdd_config = DUPLEXMODE_FDD + TXRXSWITCH_LSB;
syncmode = SYNCMODE_SLAVE;
%syncmode = SYNCMODE_FREE;
rffe_rxg_low = 31*active_rf; %[1 1 1 1];
rffe_rxg_final = 63*active_rf; %[1 1 1 1];
rffe_band = B19G_TDD*active_rf; %[1 1 1 1];
oarf_config_exmimo(card,freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNBflag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal,resampling_factor);
end
amp = pow2(14)-1;
s1p = 2*floor(amp*(s1./max([real(s1(:)); imag(s1(:))])));
s2p = 2*floor(amp*(s2./max([real(s2(:)); imag(s2(:))])));
s3p = 2*floor(amp*(s3./max([real(s3(:)); imag(s3(:))])));
%for card=min(6,num_cards-1):-1:6
% oarf_send_frame(card,s3p.',16);
%end
%for card=min(5,num_cards-1):-1:4
% oarf_send_frame(card,s2p.',16);
%end
%for card=min(3,num_cards-1):-1:0
% oarf_send_frame(card,s1p.',16);
%end
if card_select(7)
oarf_send_frame(6,s2p(3:4,:).',16);
end
if card_select(6)
oarf_send_frame(5,[zeros(153600,1) (s2p(1:2,:).')],16);
end
if card_select(5)
oarf_send_frame(4,s3p.',16);
end
if card_select(4)
oarf_send_frame(3,s1p(4,:).',16);
end
if card_select(3)
oarf_send_frame(2,s1p(3,:).',16);
end
if card_select(2)
oarf_send_frame(1,s1p(2,:).',16);
end
if card_select(1)
oarf_send_frame(0,s1p(1,:).',16);
end
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment