Commit 20d52e19 authored by jiangx's avatar jiangx

ICC 2015 version

git-svn-id: http://svn.eurecom.fr/openair4G/trunk@5910 818b1a75-f10b-46b9-bf7c-635c3b92a50f
parent 9465ca98
%
% PURPOSE : channel estimation using least square method
%
% ARGUMENTS :
%
% m_sym_T : transmitted symbol, d_N_f x d_N_ofdm x d_N_ant_act x d_N_meas
% m_sym_R : received symbol, d_N_f x d_N_ofdm x d_N_ant_act x d_N_meas
% d_N_meas : number of measurements
%
% OUTPUTS :
%
% m_H_est : estimation of sub-channels, d_N_antR x d_N_antT x d_N_f x d_N_meas
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-29-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the function "runmeas_full_duplex" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
function m_H_est = f_ch_est(m_sym_T, m_sym_R)
%% ** initialisation **
[d_N_f,d_N_OFDM,d_N_antT,d_N_meas] = size(m_sym_T);
d_N_antR = size(m_sym_R,3);
m_H_est = zeros(d_N_antR,d_N_antT,d_N_f,d_N_meas);
%% ** estimate the subband channel for each measurement and antenna **
for d_n_meas = 1:d_N_meas
for d_n_f = 1:d_N_f
m_y = reshape(m_sym_R(d_n_f,:,:,d_n_meas),d_N_OFDM,d_N_antR).'; % squeeze: problem for antenna number = 1 case
m_s = reshape(m_sym_T(d_n_f,:,:,d_n_meas),d_N_OFDM,d_N_antT).';
m_H_est(:,:,d_n_f,d_n_meas) = m_y*m_s'/(m_s*m_s'); % LS channel estimation
end
end
end
function m_sig_T = f_ofdm_mod(m_sym_T, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rf, d_amp)
d_N_ant_act = sum(v_active_rf);
%** mapping useful data to favorable carriers **
m_sym_T_ext = zeros(d_N_FFT,d_N_OFDM,d_N_ant_act);
m_sym_T_ext(2:151,:,:) = m_sym_T(1:150,:,:);
m_sym_T_ext(362:512,:,:) = m_sym_T(151:301,:,:);
%** ifft **
m_sig_T_ = sqrt(d_N_FFT)*ifft(m_sym_T_ext,d_N_FFT,1);
%** add cyclic prefix **
m_sig_T_ = [m_sig_T_(end-d_N_CP+1:end,:,:); m_sig_T_];
d_L = (d_N_FFT+d_N_CP)*d_N_OFDM;
m_sig_T = floor(reshape(m_sig_T_,d_L,d_N_ant_act)*d_amp);
end
%
% PURPOSE : OFDM Receiver
%
% ARGUMENTS :
%
% m_sig_R : received signal with dimension ((d_N_FFT+d_N_CP)*d_N_ofdm) x d_N
% d_N_FFT : total carrier number
% d_N_CP : extented cyclic prefix
% d_N_OFDM : OFDM symbol number per frame
% v_active_rf : active RF antenna indicator
%
% OUTPUTS :
%
% m_sym_R : transmitted signal before IFFT with dimension d_N_f x d_N_ofdm x d_N_ant_act
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-29-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the function "genrandpskseq" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
function m_sym_R = f_ofdm_rx(m_sig_R, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rf)
d_N_ant_act = sum(v_active_rf);
m_sig_R_eff = m_sig_R(:,find(v_active_rf));
m_sig_R_f = reshape(m_sig_R_eff,(d_N_FFT+d_N_CP),d_N_OFDM,d_N_ant_act);
%** delete the CP **
m_sig_R_noCP = m_sig_R_f(d_N_CP+1:end,:,:);
%** fft **
m_sym_R_fft = 1/sqrt(d_N_FFT)*fft(m_sig_R_noCP,d_N_FFT,1);
%m_sym_R_fft = fft(m_sig_R_noCP,d_N_FFT,1);
m_sym_R = m_sym_R_fft([2:151 362:512],:,:);
end
%
% PURPOSE : OFDM Transmitter
%
% ARGUMENTS :
%
% d_M : modulation order
% d_N_f : carrier number carrying data
% d_N_FFT : total carrier number
% d_N_CP : extented cyclic prefix
% d_N_OFDM : OFDM symbol number per frame
% v_active_rf : active RF antenna indicator
% d_amp : amplitude
%
% OUTPUTS :
%
% m_sym_T : transmitted signal before IFFT with dimension d_N_f x d_N_OFDM x d_N_ant_act
% m_sig_T : OFDM signal with dimension ((d_N_FFT+d_N_CP)*d_N_OFDM) x d_N_ant
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-29-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the function "genrandpskseq" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
function [m_sym_T, m_sig_T] = f_ofdm_tx(d_M, d_N_f, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rf, d_amp)
d_N_ant_act = sum(v_active_rf);
%** constellation table **
v_MPSK = exp(sqrt(-1)*([1:d_M]*2*pi/d_M+pi/d_M));
%** transmitted symbol **
%v_state = [1;2;3;4];
%rand("state",v_state);
m_sym_T = v_MPSK(ceil(rand(d_N_f, d_N_OFDM, d_N_ant_act)*d_M));
%** mapping useful data to favorable carriers **
m_sym_T_ext = zeros(d_N_FFT,d_N_OFDM,d_N_ant_act);
m_sym_T_ext(2:151,:,:) = m_sym_T(1:150,:,:);
m_sym_T_ext(362:512,:,:) = m_sym_T(151:301,:,:);
%** ifft **
m_sig_T_ = sqrt(d_N_FFT)*ifft(m_sym_T_ext,d_N_FFT,1);
%** add cyclic prefix **
m_sig_T_ = [m_sig_T_(end-d_N_CP+1:end,:,:); m_sig_T_];
d_L = (d_N_FFT+d_N_CP)*d_N_OFDM;
m_sig_T = floor(reshape(m_sig_T_,d_L,d_N_ant_act)*d_amp);
end
%
% PURPOSE : TLS solution for AX = B based on alternative projection
%
% ARGUMENTS :
%
% A : observation of A
% B : observation of B
%
% OUTPUTS :
%
% X : TLS solution for X
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Mai-05-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - none.
%
%**********************************************************************************************
function [X_est A_est B_est] = f_tls_ap(A,B)
%** initlisation **
e_new = 0;
e_old = 1e14;
e_thr = 1e-5; %error threshold: what if the error cannot fall down under the error threshold
X_est = eye(size(A,2));
A_est = A;
%** alternative projection **
while(abs(e_new-e_old)>e_thr)
e_old = e_new;
% optimise X_est
X_est = (A_est'*A_est)\A_est'*B;
%optimise A_est
A_est = B*X_est'/(X_est*X_est');
e_new = norm(A_est*X_est-B)^2+norm(A_est-A)^2;
end
B_est = A_est*X_est;
end
%
% PURPOSE : TLS solution for AX = B based on SVD assuming X is diagonal
%
% ARGUMENTS :
%
% A : observation of A
% B : observation of B
%
% OUTPUTS :
%
% X : TLS solution for X (Diagonal)
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-30-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - I. Markovsky and S. V. Huffel, “Overview of total least-squares methods,” Signal Processing, vol. 87, pp.
% 2283–2302, 2007
%
%**********************************************************************************************
function [X_est A_est B_est] = f_tls_diag(A,B)
d_N = size(A,2);
X_est = zeros(d_N);
A_est = zeros(size(A));
B_est = zeros(size(B));
err_est = zeros(d_N);
for d_n = 1:d_N
[X_est(d_n,d_n) A_est(:,d_n) B_est(:,d_n)] = f_tls_svd(A(:,d_n),B(:,d_n));
end
%method 2: LS solution
% for d_n = 1:d_N
% X_est(d_n,d_n) = (A(:,d_n).'*A(:,d_n))\A(:,d_n).'*B(:,d_n);
% end
%
% PURPOSE : TLS solution for AX = B based on SVD
%
% ARGUMENTS :
%
% A : observation of A
% B : observation of B
%
% OUTPUTS :
%
% X : TLS solution for X
%
%**********************************************************************************************
% EURECOM - All rights reserved
%
% AUTHOR : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-30-2014 X. JIANG 0.1 creation of code
%
% REFERENCES/NOTES/COMMENTS :
%
% - I. Markovsky and S. V. Huffel, Overview of total least-squares methods, Signal Processing, vol. 87, pp.
% 22832302, 2007
%
%**********************************************************************************************
function [X_est A_est B_est]= f_tls_svd(A,B)
C = [A B];
n = size(A,2);
d = size(B,2);
[U S V] = svd(C,0);
V12 = V(1:n,n+1:end);
V22 = V(n+1:end,n+1:end);
S1 = S(1:n,1:n);
Z12 = zeros(n,d);
Z22 = zeros(d);
Z21 = zeros(d,n);
X_est = - V12/V22;
C_est = U*[S1 Z12;Z21 Z22]*V';
A_est = C_est(:,1:n);
B_est = C_est(:,n+1:end);
% delta_C = -U*diag([zeros(n,1);diag(S(n+1:end,n+1:end))])*V';
% delta_C = [A_est-A B_est-B]; %same as the previous calculation
% err = norm(delta_C,'fro')^2/norm(C,'fro')^2;:w
end
%
% SCRIPT ID : s_beamforming
%
% PROJECT NAME : TDD Recoprocity
%
% PURPOSE : perform beamforming based on TDD calibration
%
%**********************************************************************************************
% Eurecom - All rights reserved
%
% AUTHOR(s) : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-30-2014 X. JIANG 0.1 script creation v0.1
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the script "beamforming" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
%% -------- initilisation --------
d_M = 4; % modulation order, e.g. 4 means QPSK
%** frequency **
d_fc = 1907600000;
d_delta_f = 15000;
d_N_f = 301; % carrier number carrying data
d_N_FFT = 512; % total carrier number
d_N_CP = 128; % extented cyclic prefix
%** time **
d_N_OFDM = 60; % number of ofdm symbol per half frame
d_N_meas = 1; % measuement number
%** space **
d_N_antA = 4; % antenna number at site a
d_N_antB = 4; % antenna number at site b
v_active_rfA = [1 1 1 1];
v_active_rfB = [1 0 0 0];
v_indA = find(v_active_rfA); % active antenna index at site a
v_indB = find(v_active_rfB); % active antenna index at site b
d_N_antA_act = length(v_indA);
d_N_antB_act = length(v_indB);
%** amplitude **
d_amp = pow2(12.5)-1;
%% -------- load F --------
o_result = load('result/4a_l45_t10b.mat');
m_F_full = o_result.m_F0;
m_F_diag = o_result.m_F2;
%% -------- channel measurement --------
s_run_meas;
%% -------- signal precoding --------
v_MPSK = exp(sqrt(-1)*([1:d_M]*2*pi/d_M+pi/d_M));
m_sym_TA = v_MPSK(ceil(rand(d_N_f, d_N_OFDM*2)*d_M));
m_sym_TA_ideal = zeros(d_N_f,d_N_OFDM*2,d_N_antA_act);
for d_f = 1:d_N_f
%** ideal case **
v_H_A2B_ideal = squeeze(m_H_est_A2B(:,:,d_f));
v_P_ideal = v_H_A2B_ideal'/norm(v_H_A2B_ideal);
m_sym_TA_ideal(d_f,:,:) = (v_P_ideal*m_sym_TA(d_f,:)).';
%** identity matrix **
v_H_A2B_iden = squeeze(m_H_est_B2A(:,:,d_f)).';
v_P_iden = v_H_A2B_iden'/norm(v_H_A2B_iden);
m_sym_TA_iden(d_f,:,:) = (v_P_iden*m_sym_TA(d_f,:)).';
%** diagonal calibration **
v_H_A2B_diag = squeeze(m_H_est_B2A(:,:,d_f).')*squeeze(m_F_diag(:,:,2));
v_P_diag = v_H_A2B_diag'/norm(v_H_A2B_diag);
m_sym_TA_diag(d_f,:,:) = (v_P_diag*m_sym_TA(d_f,:)).';
%** full calibration **
v_H_A2B_full = squeeze(m_H_est_B2A(:,:,d_f).')*squeeze(m_F_diag(:,:,2));
v_P_full = v_H_A2B_full'/norm(v_H_A2B_full);
m_sym_TA_full(d_f,:,:) = (v_P_full*m_sym_TA(d_f,:)).';
end
%% -------- signal transmission --------
m_sig_TA_ideal = f_ofdm_mod(m_sym_TA_ideal,d_N_FFT,d_N_CP,d_N_OFDM*2,v_active_rfA,d_amp)*2;
m_sig_TA_iden = f_ofdm_mod(m_sym_TA_iden,d_N_FFT,d_N_CP,d_N_OFDM*2,v_active_rfA,d_amp)*2;
m_sig_TA_diag = f_ofdm_mod(m_sym_TA_diag,d_N_FFT,d_N_CP,d_N_OFDM*2,v_active_rfA,d_amp)*2;
m_sig_TA_full = f_ofdm_mod(m_sym_TA_full,d_N_FFT,d_N_CP,d_N_OFDM*2,v_active_rfA,d_amp)*2;
m_sig_TB = zeros((d_N_CP+d_N_FFT)*d_N_OFDM*2,d_N_antB);
oarf_send_frame(cardB,m_sig_TB,d_n_bit);
m_noise_RB_ = oarf_get_frame(cardB);
m_noise_RB = m_noise_RB_(1:d_N_OFDM*(d_N_FFT+d_N_CP),5:8);
m_n_sym_RB = f_ofdm_rx(m_noise_RB, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
m_sig_TB(:,:)=1+1i;
oarf_send_frame(cardB,m_sig_TB,d_n_bit);
oarf_send_frame(cardA,m_sig_TA_ideal,d_n_bit);
m_sig_RB_ideal_ = oarf_get_frame(cardB);
m_sig_RB_ideal = m_sig_RB_ideal_(1:d_N_OFDM*(d_N_FFT+d_N_CP),5:8);
m_sym_RB_ideal = f_ofdm_rx(m_sig_RB_ideal, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
oarf_send_frame(cardA,m_sig_TA_iden,d_n_bit);
m_sig_RB_iden_ = oarf_get_frame(cardB);
m_sig_RB_iden = m_sig_RB_iden_(1:d_N_OFDM*(d_N_FFT+d_N_CP),5:8);
m_sym_RB_iden = f_ofdm_rx(m_sig_RB_iden, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
oarf_send_frame(cardA,m_sig_TA_diag,d_n_bit);
m_sig_RB_diag_ = oarf_get_frame(cardB);
m_sig_RB_diag = m_sig_RB_diag_(1:d_N_OFDM*(d_N_FFT+d_N_CP),5:8);
m_sym_RB_diag = f_ofdm_rx(m_sig_RB_diag, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
oarf_send_frame(cardA,m_sig_TA_full,d_n_bit);
m_sig_RB_full_ = oarf_get_frame(cardB);
m_sig_RB_full = m_sig_RB_full_(1:d_N_OFDM*(d_N_FFT+d_N_CP),5:8);
m_sym_RB_full = f_ofdm_rx(m_sig_RB_full, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
%% -------- SNR measurement --------
%** noise measurment **
%v_P_n = mean(var(squeeze(m_n_sym_RB),0,2));
v_P_n = mean(var(squeeze(m_n_sym_RB),0,2));
%** SNR caculation
%v_P_s_ideal = zeros(301,1);
%for d_f=1:d_N_f
% v_H_A2B_ideal = squeeze(m_H_est_A2B(:,:,d_f));
% v_P_s_ideal(d_f) = norm(v_H_A2B_ideal)^2;
%end
%keyboard;
v_P_s_ideal = var(squeeze(m_sym_RB_ideal),0,2);
v_P_s_iden = var(squeeze(m_sym_RB_iden),0,2);
v_P_s_diag = var(squeeze(m_sym_RB_diag),0,2);
v_P_s_full = var(squeeze(m_sym_RB_full),0,2);
v_SNR_ideal_ = 10*log10((v_P_s_ideal-v_P_n)/v_P_n);
v_SNR_iden_ = 10*log10((v_P_s_iden-v_P_n)/v_P_n);
v_SNR_diag_ = 10*log10((v_P_s_diag-v_P_n)/v_P_n);
v_SNR_full_ = 10*log10((v_P_s_full-v_P_n)/v_P_n);
v_SNR_ideal = nan(d_N_f+1,1);
v_SNR_iden = nan(d_N_f+1,1);
v_SNR_diag = nan(d_N_f+1,1) ;
v_SNR_full = nan(d_N_f+1,1) ;
v_SNR_ideal([1:151 153:302]) = v_SNR_ideal_([151:301 1:150]);
v_SNR_iden([1:151 153:302]) = v_SNR_iden_([151:301 1:150]) ;
v_SNR_diag([1:151 153:302]) = v_SNR_diag_([151:301 1:150]) ;
v_SNR_full([1:151 153:302]) = v_SNR_full_([151:301 1:150]);
save('-v7','result/bf_gain_4x1_t3.mat','v_SNR_ideal','v_SNR_iden','v_SNR_diag','v_SNR_full');
%% -------- plot --------
v_f = d_fc-floor(d_N_f/2)*d_delta_f:d_delta_f:d_fc+ceil(d_N_f/2)*d_delta_f;
figure(5)
hold on
plot(v_f,v_SNR_ideal,'k-')
plot(v_f,v_SNR_iden,'m-.')
plot(v_f,v_SNR_diag,'r-')
plot(v_f,v_SNR_full,'b-')
hold off
ylim([30 40])
%
% SCRIPT ID : s_run_calib
%
% PROJECT NAME : TDD Recoprocity
%
% PURPOSE : channel calibration for MISO case
%
%**********************************************************************************************
% Eurecom - All rights reserved
%
% AUTHOR(s) : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-30-2014 X. JIANG 0.1 script creation v0.1
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the script "calibration" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
%% ** initilisation **
%---------- to change in experiement ---------
%s_init_params;
%clc
%clear all
close all
%d_N_f = 301; % carrier number carrying data
%d_N_meas = 10; % measuement number
%d_N_loc = 5; % Rx locations
%d_N_antM = 2; % max active antenna number for site a and site b
%----------------------------------------------
%% -------- System parameters --------
d_M = 4; % modulation order, e.g. 4 means QPSK
%** frequency **
d_N_f = 301; % carrier number carrying data
d_N_FFT = 512; % total carrier number
d_N_CP = 128; % extented cyclic prefix
%** time **
d_N_OFDM = 60; % number of ofdm symbol per half frame
d_N_meas = 10; % measuement number
%** space **
d_N_antA = 4; % antenna number at site a
d_N_antB = 4; % antenna number at site b
v_indA = find(v_active_rfA); % active antenna index at site a
v_indB = find(v_active_rfB); % active antenna index at site b
d_N_antA_act = length(v_indA);
d_N_antB_act = length(v_indB);
%** amplitude **
d_amp = pow2(12.5)-1; % to see how to be used??
%% -------- calibration parameters -------
d_N_loc = 45; % Rx locations
d_N_antM = max(sum(v_active_rfA),sum(v_active_rfB)); % max active antenna number for site a and site b
%% -------- initialisation ----------------
m_H_A2B = zeros(d_N_antM,d_N_meas*d_N_loc, d_N_f); % d_N_antA x (d_N_meas*d_N_loc) x d_N_f
m_H_B2A = zeros(d_N_antM,d_N_meas*d_N_loc, d_N_f); % d_N_antA x (d_N_meas*d_N_loc) x d_N_f
m_F0 = zeros(d_N_antM,d_N_antM,d_N_f);
m_F1 = zeros(d_N_antM,d_N_antM,d_N_f);
%% ** collect the measurement data from different locations **
for d_loc = 1:d_N_loc
% run measurement, note: uncomment "clear all"
s_run_meas;
% --- the following part is dedicated to A2B MISO -----
m_H_A2Bi = permute(squeeze(m_H_est_A2B),[1 3 2]);
m_H_B2Ai = permute(squeeze(m_H_est_B2A),[1 3 2]);
% -----------------------------------------------------
m_H_A2B(:,(d_loc-1)*d_N_meas+1:d_loc*d_N_meas,:) = m_H_A2Bi;
m_H_B2A(:,(d_loc-1)*d_N_meas+1:d_loc*d_N_meas,:) = m_H_B2Ai;
%keyboard;
disp('Please move the antenna to another location and press any key to continue');
% pause
pause(5);
end
save('-v7','result/4a_l45_t10d.mat','m_H_A2B','m_H_B2A');
%save('-v7','result/2c_l15_t2a.mat','m_H_A2B','m_H_B2A');
%% ** calibration **
for d_f = 1:d_N_f
[m_F0(:,:,d_f),m_A0_est,m_B0_est] = f_tls_svd(m_H_B2A(:,:,d_f).',m_H_A2B(:,:,d_f).'); % solve the TLS problem using SVD method
[m_F1(:,:,d_f),m_A1_est,m_B0_est] = f_tls_ap(m_H_B2A(:,:,d_f).',m_H_A2B(:,:,d_f).'); % solve the TLS problem using Alternative Projection method
end
%% ** plot **
figure(10)
hold on;
for d_f=1:size(m_F0,3);
m_F= m_F0(:,:,d_f);
plot(diag(m_F),'bo')
plot(diag(m_F,1),'k+')
plot(diag(m_F,2),'rx')
plot(diag(m_F,3),'g*')
plot(diag(m_F,-1),'y+')
plot(diag(m_F,-2),'mx')
plot(diag(m_F,-3),'c*')
end
hold off;
title('F0');
axis([-2 2 -2 2])
grid on
figure(11)
hold on;
for d_f=1:size(m_F1,3);
m_F= m_F1(:,:,d_f);
plot(diag(m_F),'bo')
plot(diag(m_F,1),'k+')
plot(diag(m_F,2),'rx')
plot(diag(m_F,3),'g*')
plot(diag(m_F,-1),'y+')
plot(diag(m_F,-2),'mx')
plot(diag(m_F,-3),'c*')
end
hold off;
title('F1');
axis([-2 2 -2 2])
grid on
close all
%clear all
d_N_f = 301; % carrier number carrying data
d_N_meas = 10;
d_file_name = 'result/4a_l45_t10a.mat';
o_result = load(d_file_name);
m_H_B2A = o_result.m_H_B2A;
m_H_A2B = o_result.m_H_A2B;
d_N_loc = size(m_H_A2B,2)/d_N_meas;
m_F0 = zeros(size(m_H_B2A,1));
m_F1 = zeros(size(m_H_B2A,1));
m_F2 = zeros(size(m_H_B2A,1));
%% ** post processing **
% ** normalisation **
for d_l = 1:d_N_loc
temp = (d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas;
d_norm_fac = mean(mean(mean(abs(m_H_B2A(:,temp,:)))));
m_H_B2A(:,temp,:) = m_H_B2A(:,temp,:)/d_norm_fac;
m_H_A2B(:,temp,:) = m_H_A2B(:,temp,:)/d_norm_fac;
end
% ** average **
d_var_thr = 0.03;
m_H_A2B_ = [];
m_H_B2A_ = [];
for d_l = 1:d_N_loc
% d_var_l = max(var(squeeze(m_H_A2B(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,1)).'));
d_var_l = max(var(m_H_A2B(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,1).'));
d_mean_l = mean(mean(abs(m_H_A2B(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,1).')));
if d_var_l/d_mean_l < d_var_thr
% m_H_A2B_ = [m_H_A2B_ mean(m_H_A2B(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,:),2)];
% m_H_B2A_ = [m_H_B2A_ mean(m_H_B2A(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,:),2)];
m_H_A2B_ = [m_H_A2B_ m_H_A2B(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,:)];
m_H_B2A_ = [m_H_B2A_ m_H_B2A(:,(d_l-1)*d_N_meas+1:(d_l-1)*d_N_meas+d_N_meas,:)];
end
end
%keyboard;
% % ** no processing **
% m_H_B2A_ = m_H_B2A;
% m_H_A2B_ = m_H_A2B;
%% ** calibration **
% temp = 1:15;
for d_f = 1:d_N_f
[m_F0(:,:,d_f),m_A0_est,m_B0_est] = f_tls_svd(m_H_B2A_(:,:,d_f).',m_H_A2B_(:,:,d_f).');
% [m_F0(:,:,d_f),m_A0_est,m_B0_est,d_err0] = f_tls_svd(m_H_B2A_(:,temp,d_f).',m_H_A2B_(:,temp,d_f).');
% [m_F1(:,:,d_f),m_A1_est,m_B1_est d_err1] = f_tls_ap(m_H_B2A_(:,:,d_f).',m_H_A2B_(:,:,d_f).');
[m_F2(:,:,d_f),m_A2_est,m_B2_est] = f_tls_diag(m_H_B2A_(:,:,d_f).',m_H_A2B_(:,:,d_f).');
end
save('-v7','-append',d_file_name,'m_F0','m_F2');
%% ** plot **
h10 = figure(10);
hold on;
for d_f=1:d_N_f
m_F= m_F0(:,:,d_f);
plot(diag(m_F),'bo')
plot(diag(m_F,1),'k+')
plot(diag(m_F,-1),'y+')
plot(diag(m_F,2),'rx')
plot(diag(m_F,-2),'mx')
plot(diag(m_F,3),'g*')
plot(diag(m_F,-3),'c*')
end
hold off;
title('F0');
% axis([-0.1 0.1 -0.1 0.1])
axis([-2 2 -2 2])
% axis([-0.3 0.3 -0.3 0.3])
grid on
% saveas(h10,'results\F_svd_2c_g20_0_t1','fig');
% print('-depsc','results\F_svd_2c_g20_0_t1.eps')
%
% h11=figure(11);
% hold on;
% for d_f=1:d_N_f
% m_F= m_F1(:,:,d_f);
% plot(diag(m_F),'bo')
% plot(diag(m_F,1),'k+')
% plot(diag(m_F,-1),'y+')
% plot(diag(m_F,2),'rx')
% plot(diag(m_F,-2),'mx')
% plot(diag(m_F,3),'g*')
% plot(diag(m_F,-3),'c*')
% end
% hold off;
% axis([-0.3 0.3 -0.3 0.3])
% % axis([-2 2 -2 2])
% grid on
% saveas(h11,'results\F_ap_2c_g20_0_t2','fig');
% print('-depsc','results\F_ap_2c_g20_0_t2.eps')
%
% h12=figure(12);
% hold on;
% for d_f=1:d_N_f
% m_F= m_F2(:,:,d_f);
% plot(diag(m_F),'bo')
% % plot(diag(m_F,1),'k+')
% % plot(diag(m_F,-1),'k+')
% % plot(diag(m_F,2),'rx')
% % plot(diag(m_F,-2),'rx')
% % plot(diag(m_F,3),'g*')
% % plot(diag(m_F,-3),'g*')
% end
% hold off;
% title('F2');
% axis([-2 2 -2 2])
% grid on
% saveas(h12,'results\F_diag_2c_g20_0_t2','fig');
% print('-depsc','results\F_diag_2c_g20_0_t2.eps')
clear all
close all
addpath([getenv('OPENAIR_TARGETS') '/ARCH/EXMIMO/USERSPACE/OCTAVE']);
%% -------- ExpressMIMO2 configuration --------
limeparms;
cardA = 0;
cardB = 1;
v_active_rfA = [1 1 1 1];
v_active_rfB = [1 0 0 0];
fc = 1907600000; %3500000000; %fc = 859.5e6; %fc = 2660000000;
freq_rxA = fc*v_active_rfA;
freq_txA = freq_rxA; %+1.92e6;
freq_rxB = fc*v_active_rfB;
freq_txB = freq_rxB; %+1.92e6;
tdd_config = DUPLEXMODE_FDD+TXRXSWITCH_LSB; %we need the LSB switching for the woduplex script, otherwise we don't receive anything
rx_gainA = 10*v_active_rfA;%[0 0 0 3];
tx_gainA = 10*v_active_rfA;%[20 20 20 20]
rx_gainB = 10*v_active_rfB;
tx_gainB = 10*v_active_rfB;
syncmodeA = SYNCMODE_MASTER;
syncmodeB = SYNCMODE_SLAVE;
eNB_flag = 0;
resampling_factorA = [2 2 2 2];%2*v_active_rfA;
resampling_factorB = [2 2 2 2];%2*v_active_rfB;
rf_modeA = (TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM+DMAMODE_TX+TXEN+DMAMODE_RX+RXEN) * v_active_rfA;
rf_rxdcA = rf_rxdc*v_active_rfA;
rf_vcocalA = rf_vcocal_19G*v_active_rfA;
rf_local = [8254744 8255063 8257340 8257340]; %eNB2tx 1.9GHz
rffe_rxg_lowA = 31*v_active_rfA;
rffe_rxg_finalA = 63*v_active_rfA;
rffe_bandA = B19G_TDD*v_active_rfA;
autocal_modeA = v_active_rfA;
oarf_stop(cardA);
sleep(0.1);
oarf_config_exmimo(cardA,freq_rxA,freq_txA,tdd_config,syncmodeA,rx_gainA,tx_gainA,eNB_flag,rf_modeA,rf_rxdcA,rf_local,rf_vcocalA,rffe_rxg_lowA,rffe_rxg_finalA,rffe_bandA,autocal_modeA,resampling_factorA);
rf_modeB = (TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM+DMAMODE_TX+TXEN+DMAMODE_RX+RXEN) * v_active_rfB;
rf_rxdcB = rf_rxdc*v_active_rfB;
rf_vcocalB = rf_vcocal_19G*v_active_rfB;
rffe_rxg_lowB = 31*v_active_rfB;
rffe_rxg_finalB = 63*v_active_rfB;
rffe_bandB = B19G_TDD*v_active_rfB;
autocal_modeB = v_active_rfB;
oarf_stop(cardB);
sleep(0.1);
oarf_config_exmimo(cardB,freq_rxB,freq_txB,tdd_config,syncmodeB,rx_gainB,tx_gainB,eNB_flag,rf_modeB,rf_rxdcB,rf_local,rf_vcocalB,rffe_rxg_lowB,rffe_rxg_finalB,rffe_bandB,autocal_modeB,resampling_factorB);
d_n_bit = 16;
%
% SCRIPT ID : s_run_meas
%
% PROJECT NAME : TDD Recoprocity
%
% PURPOSE : full transmission and receive train for TDD reciprocity calibration
%
%**********************************************************************************************
% Eurecom - All rights reserved
%
% AUTHOR(s) : Xiwen JIANG, Florian Kaltenberger
%
% DEVELOPMENT HISTORY :
%
% Date Name(s) Version Description
% ----------- ------------- ------- ------------------------------------------------------
% Apr-29-2014 X. JIANG 0.1 script creation v0.1
%
% REFERENCES/NOTES/COMMENTS :
%
% - Based on the script "run_full_duplex" created by Mirsad Cirkic, Florian Kaltenberger.
%
%**********************************************************************************************
%% ** initialisation **
%% ------------- to change in experiment ------------
%clc
close all
%clear all
%
%% ----------------------------------------------------
m_sym_TA = zeros(d_N_f,d_N_OFDM,length(v_indA),d_N_meas);
m_sym_TB = zeros(d_N_f,d_N_OFDM,length(v_indB),d_N_meas);
m_sym_RA = zeros(d_N_f,d_N_OFDM,length(v_indA),d_N_meas);
m_sym_RB = zeros(d_N_f,d_N_OFDM,length(v_indB),d_N_meas);
m_sym_TA_ = zeros(d_N_f,d_N_OFDM,length(v_indA));
m_sig_TA = zeros((d_N_CP+d_N_FFT)*d_N_OFDM*2,d_N_antA);
m_sig_TB = zeros((d_N_CP+d_N_FFT)*d_N_OFDM*2,d_N_antB);
for d_n_meas = 1:d_N_meas
%% -------- tx --------
%** tx of site A **
[m_sym_TA_(:,:,:,d_n_meas), m_sig_TA_] = f_ofdm_tx(d_M, d_N_f, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfA, d_amp);
for d_a = 1:d_N_antA_act
d_N_a = d_N_OFDM/d_N_antA_act;
d_N_b = d_N_OFDM*(d_N_FFT+d_N_CP)/d_N_antA_act;
m_sig_TA(d_N_b*(d_a-1)+1:d_N_b*d_a,v_indA(d_a)) = m_sig_TA_(d_N_b*(d_a-1)+1:d_N_b*d_a,d_a);
m_sym_TA(:,d_N_a*(d_a-1)+1:d_N_a*d_a,d_a,d_n_meas) = m_sym_TA_(:,d_N_a*(d_a-1)+1:d_N_a*d_a,d_a,d_n_meas);
end
m_sig_TA(1:end/2,v_indA)= m_sig_TA(1:end/2,v_indA)*2; %affect the LSB to 0 to set on Tx mode
m_sig_TA(end/2+1:end,v_indA)= 1+1i; %affect the LSB to 1 to set on Rx mode
%** tx of site B **
[m_sym_TB(:,:,:,d_n_meas), m_sig_TB_] = f_ofdm_tx(d_M, d_N_f, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB, d_amp);
m_sig_TB(1:end/2,v_indB)=1+1i;
m_sig_TB(end/2+1:end,v_indB)= m_sig_TB_*2;
%** Transmission from A to B **
oarf_send_frame(cardA,m_sig_TA,d_n_bit);
%** Transmission from B to A **
oarf_send_frame(cardB,m_sig_TB,d_n_bit);
%sleep(0.1);
m_sig_R = oarf_get_frame(-2);
m_sig_RA = m_sig_R(d_N_OFDM*(d_N_FFT+d_N_CP)+1:d_N_OFDM*(d_N_FFT+d_N_CP)*2,1:d_N_antA);
m_sig_RB = m_sig_R(1:d_N_OFDM*(d_N_FFT+d_N_CP),d_N_antA+1:d_N_antA+d_N_antB);
%% -------- rx --------
m_sym_RB(:,:,:,d_n_meas) = f_ofdm_rx(m_sig_RB, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfB);
m_sym_RA(:,:,:,d_n_meas) = f_ofdm_rx(m_sig_RA, d_N_FFT, d_N_CP, d_N_OFDM, v_active_rfA);
%keyboard;
end
%** channel estimation **
m_H_est_A2B = f_ch_est(m_sym_TA, m_sym_RB); %dimension: d_N_antR x d_N_antT x d_N_f x d_N_meas
m_H_est_B2A = f_ch_est(m_sym_TB, m_sym_RA);
%% -------- plot --------
%** channel estimation in frequency domain **
m_H_A2B_draw = squeeze(m_H_est_A2B(1,:,:,end)).';
m_H_B2A_draw = squeeze(m_H_est_B2A(:,1,:,end)).';
figure(1)
subplot(2,1,1)
plot(real(m_sig_RB(1:40,v_indB)),'b-');
hold on
plot(real(m_sig_RB(end-80+1:end-40,v_indB)),'g-')
hold on
plot(real(m_sig_RB(end-40+1:end,v_indB)),'r-')
title('m_sig_RB')
subplot(2,1,2)
plot(real(m_sig_RA(:,v_indA)),'-');
title('m_sig_RA')
figure(2)
subplot(2,2,1)
plot(20*log10(abs(m_H_A2B_draw)),'-');
title('|h| vs. freq (A2B)')
xlabel('freq')
ylabel('|h|')
ylim([0 100])
subplot(2,2,2)
plot(20*log10(abs(m_H_B2A_draw)),'-');
title('|h| vs. freq (B2A)')
xlabel('freq')
ylabel('|h|')
ylim([0 100])
subplot(2,2,3)
plot(angle(m_H_A2B_draw),'-');
title('angle(h) vs. freq (A2B)')
xlabel('freq')
ylabel('angle(h)')
subplot(2,2,4)
plot(angle(m_H_B2A_draw),'-');
title('angle(h) vs. freq (B2A)')
xlabel('freq')
ylabel('angle(h)')
%v_color = ['b*','g*','r*','c*'];
figure(3)
for d_a = 1:d_N_antA_act
subplot(2,2,d_a);
plot(m_sym_RA(1,:,d_a,1),'*');
title(strcat('m sym RA',num2str(d_a)));
end
figure(4)
for d_a = 1:d_N_antA_act
subplot(2,2,d_a);
plot(m_sym_RB(1,d_N_OFDM/d_N_antA_act*(d_a-1)+1:d_N_OFDM/d_N_antA_act*d_a,1,1),'*');
hold on;
title(strcat('m sym RB',num2str(d_a)));
end
...@@ -26,4 +26,4 @@ for i=0:(N/640-1) ...@@ -26,4 +26,4 @@ for i=0:(N/640-1)
s([1:640]+i*640)=floor(amp*block); s([1:640]+i*640)=floor(amp*block);
end end
end end
\ No newline at end of file
...@@ -15,8 +15,9 @@ eNB_flag = 0; ...@@ -15,8 +15,9 @@ eNB_flag = 0;
card = 0; card = 0;
Ntrx=4; Ntrx=4;
dual_tx=0; dual_tx=0;
active_rfA=[0 0 1 0]; active_rfA=[0 0 0 1];
active_rfB=[1 1 0 0]; % active_rfB=[1 1 1 0];
active_rfB=[0 1 1 0];
active_rf=active_rfA | active_rfB; active_rf=active_rfA | active_rfB;
if(active_rfA*active_rfB'!=0) error("The A and B transceive chains must be orthogonal./n") endif if(active_rfA*active_rfB'!=0) error("The A and B transceive chains must be orthogonal./n") endif
...@@ -26,6 +27,7 @@ fc = 1912600000; %1907600000; ...@@ -26,6 +27,7 @@ fc = 1912600000; %1907600000;
%fc = 859.5e6; %fc = 859.5e6;
autocal_mode=active_rf; autocal_mode=active_rf;
resampling_factor = [2 2 2 2];
%rf_mode=(RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAByp+RFBBLNA1) * active_rf; %rf_mode=(RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAByp+RFBBLNA1) * active_rf;
%rf_mode=(TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM) * active_rf; %rf_mode=(TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM) * active_rf;
% we have to enable both DMA transfers so that the switching signal in the LSB of the TX buffer gets set % we have to enable both DMA transfers so that the switching signal in the LSB of the TX buffer gets set
...@@ -47,11 +49,12 @@ freq_tx = freq_rx; %+1.92e6; ...@@ -47,11 +49,12 @@ freq_tx = freq_rx; %+1.92e6;
oarf_stop(card); oarf_stop(card);
sleep(0.1); sleep(0.1);
oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode); oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode,resampling_factor);
autocal_mode=0*active_rf; % Autocalibration is only needed the first time we conf. exmimo autocal_mode=0*active_rf; % Autocalibration is only needed the first time we conf. exmimo
amp = pow2(12.5)-1; amp = pow2(12.5)-1;
n_bit = 16; n_bit = 16;
%chanest_full = 1;
chanest_full = 1; chanest_full = 1;
paramsinitialized=true; paramsinitialized=true;
# % Author: Mirsad Cirkic
# % Organisation: Eurecom (and Linkoping University)
# % E-mail: mirsad.cirkic@liu.se
%clear
paramsinitialized=false;
limeparms;
rxgain=10;
txgain=25;
eNB_flag = 0;
card = 0;
Ntrx=4;
dual_tx=0;
active_rfA=[1 0 0 0];
active_rfB=[0 1 1 0];
active_rf=active_rfA+active_rfB;
if(active_rfA*active_rfB'!=0) error("The A and B transceive chains must be orthogonal./n") endif
%fc = 2660000000;
fc = 1912600000; %1907600000;
%fc = 859.5e6;
%rf_mode=(RXEN+TXEN+TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAByp+RFBBLNA1) * active_rf;
autocal_mode=active_rf;
rf_mode=(TXLPFNORM+TXLPFEN+TXLPF25+RXLPFNORM+RXLPFEN+RXLPF25+LNA1ON+LNAMax+RFBBNORM) * active_rf;
tdd_config = DUPLEXMODE_FDD+TXRXSWITCH_TESTTX; LSBSWITCH_FLAG=false;
%tdd_config = DUPLEXMODE_FDD+TXRXSWITCH_LSB; LSBSWITCH_FLAG=true;
syncmode = SYNCMODE_FREE;
rf_local = [8254744 8255063 8257340 8257340]; %eNB2tx 1.9GHz
rf_vcocal=rf_vcocal_19G*active_rf;
rffe_rxg_low = 61*active_rf;
rffe_rxg_final = 61*active_rf;
rffe_band = B19G_TDD*active_rf;
rf_rxdc = rf_rxdc*active_rf;
freq_rx = fc*active_rf;
freq_tx = freq_rx; %+1.92e6;
tx_gain = txgain*active_rf;
rx_gain = rxgain*active_rf;
%rx_gain = [10 10 30 0];
oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode);
autocal_mode=0*active_rf; % Autocalibration is only needed the first time we conf. exmimo
amp = pow2(14)-1;
n_bit = 16;
paramsinitialized=true;
old project by Misard and Florian
...@@ -21,7 +21,7 @@ if(paramsinitialized) ...@@ -21,7 +21,7 @@ if(paramsinitialized)
if(Niter~=1) if(Niter~=1)
error('We should only use one get_frame at each run.\n'); error('We should only use one get_frame at each run.\n');
end end
Nmeas = 10; Nmeas = 5;
%% ------- Prepare the signals for A2B ---------- %% %% ------- Prepare the signals for A2B ---------- %%
signalA2B=zeros(N,4); signalA2B=zeros(N,4);
...@@ -75,19 +75,18 @@ if(paramsinitialized) ...@@ -75,19 +75,18 @@ if(paramsinitialized)
for meas=1:Nmeas for meas=1:Nmeas
%% ------- Node A to B transmission ------- %% %% ------- Node A to B transmission ------- %%
oarf_send_frame(card,signalA2B,n_bit); oarf_send_frame(card,signalA2B,n_bit);
%keyboard
sleep(0.01);
receivedA2B=oarf_get_frame(card); receivedA2B=oarf_get_frame(card);
%oarf_stop(card); %not good, since it does a reset %oarf_stop(card); %not good, since it does a reset
sleep(0.01); sleep(0.01);
%%----------Node B to A transmission---------%% %%----------Node B to A transmission---------%%
oarf_send_frame(card,signalB2A,n_bit); oarf_send_frame(card,signalB2A,n_bit);
%keyboard
sleep(0.01);
receivedB2A=oarf_get_frame(card); receivedB2A=oarf_get_frame(card);
%oarf_stop(card); %not good, since it does a reset %oarf_stop(card); %not good, since it does a reset
%keyboard;
%% ------- Do the A to B channel estimation ------- %% %% ------- Do the A to B channel estimation ------- %%
for i=0:119; for i=0:119;
ifblock=receivedA2B(i*640+[1:640],indB); ifblock=receivedA2B(i*640+[1:640],indB);
...@@ -158,7 +157,7 @@ if(paramsinitialized) ...@@ -158,7 +157,7 @@ if(paramsinitialized)
fblock=fft(ifblock); fblock=fft(ifblock);
fblock(1,:)=[]; fblock(1,:)=[];
fblock(151:360,:)=[]; fblock(151:360,:)=[];
noise_f(i+1,:,:)=fblock; %noise_f(i+1,:,:)=fblock;
end end
...@@ -177,8 +176,9 @@ if(paramsinitialized) ...@@ -177,8 +176,9 @@ if(paramsinitialized)
end end
figure(2) figure(2)
t=[0:512-1]/512*1e-2; %t=[0:512-1]/512*1e-2;
plot(t,20*log10(abs(tchanests))) %plot(t,20*log10(abs(tchanests)))
plot(20*log10(abs(tchanests)))
xlabel('time') xlabel('time')
ylabel('|h|') ylabel('|h|')
if Nantb==3 if Nantb==3
......
# % Author: Mirsad Cirkic
# % Organisation: Eurecom (and Linkoping University)
# % E-mail: mirsad.cirkic@liu.se
if(paramsinitialized && LSBSWITCH_FLAG)
disp(["\n\n------------\nThis code is, so far, only written for single runs. Multiple " ...
"runs will overwrite the previous measurement data, i.e., the " ...
"data structures are not defined for multiple runs. You will need to " ...
"add code in order to save the intermediate measurements and the " ...
"corresponding timestamps.\n------------"])
N=76800;
M=4;
signalA2B=zeros(N,4);
signalB2A=zeros(N,4);
indA=find(active_rfA==1);
indB=find(active_rfB==1);
Nanta=length(indA);
Nantb=length(indB);
Niter=1;
Nofs=3840; % The offset in samples between the downlink and uplink sequences. It must be a multiple of 640, 3840 samples is 0.5 ms.
if(mod(Nofs,640)!=0) error("blabla") endif
if(Nanta!=1) error("Node A can only have one antenna active\n"); endif
if(Niter!=1) error("We should only use one get_frame at each run.\n"); endif
# %% ------- Prepare the signals for both A2B and B2A ------- %
maskA2B=kron(ones(1,N/(2*Nofs)),[ones(1,Nofs) zeros(1,Nofs)])';
maskB2A=ones(N,1)-maskA2B;
datamaskA2B=diag(kron(ones(1,N/(2*Nofs)),[ones(1,Nofs/640) zeros(1,Nofs/640)]));
datamaskB2A=eye(N/640)-datamaskA2B;
signalA2B=zeros(N,4);
signalB2A=zeros(N,4);
ia=1; ib=1;
Db2a_T=[];
for i=1:4
if(indA(ia)==i)
[tmpd, tmps]=generqpskseq(N,M,amp);
Da2b_T=datamaskA2B*tmpd;
signalA2B(:,i)=tmps*4.*maskA2B+maskB2A*sqrt(-1); %Added maskB2A to set the LSB correctly. The factor 4 there should make a 2 bit shift. Don't know if that is correct
if(length(indA)> ia) ia=ia+1; endif
endif
if(indB(ib)==i)
% This part could be improved by creating fully orthogonal sequences
[tmpd, tmps]=generqpskseq(N,M,amp);
signalB2A(:,i)=tmps*4.*maskB2A+maskA2B*sqrt(-1); %Added maskA2B to set the LSB correctly. The factor 4 there should make a 2 bit shift. Don't know if that is correct.
Db2a_T=[Db2a_T datamaskB2A*tmpd];
if(length(indB)> ib) ib=ib+1; endif
endif
endfor
signal = signalA2B+signalB2A;
# %% ------- Node B and A duplex transmission/reception ------- %%
rf_mode_current = rf_mode + (DMAMODE_TX+TXEN+DMAMODE_RX+RXEN)*active_rf;
oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode_current,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode);
oarf_send_frame(card,signal,n_bit);
received=oarf_get_frame(card);
oarf_stop(card);
receivedA2B=received.*repmat(maskA2B,1,4);
receivedB2A=received.*repmat(maskB2A,1,4);
# %% ------- Do the A to B channel estimation ------- %%
Da2b_R=zeros(Niter*120,Nantb*301);
for i=0:119;
ifblock=receivedA2B(i*640+[1:640],indB);
ifblock(1:128,:)=[];
fblock=fft(ifblock);
fblock(1,:)=[];
fblock(151:360,:)=[];
Da2b_R((Niter-1)*120+i+1,:)=vec(fblock);
endfor
HA2B=repmat(conj(Da2b_T),1,Nantb).*Da2b_R;
phasesA2B=unwrap(angle(HA2B));
if(mean(var(phasesA2B))>0.5)
disp("The phases of your estimates are a bit high (larger than 0.5 rad.), something is wrong.");
endif
chanestsA2B=reshape(diag(repmat(Da2b_T,1,Nantb)'*Da2b_R)/size(Da2b_T,1),301,Nantb);
fchanestsA2B=zeros(512,Nantb);
for i=1:Nantb
fchanestsA2B(:,i)=[0; chanestsA2B([1:150],i); zeros(210,1); chanestsA2B(151:301,i)];
endfor
tchanestsA2B=ifft(fchanestsA2B);
%% ------- Do the B to A channel estimation ------- %%
Db2a_R=zeros(Niter*120,Nanta*301);
for i=0:119;
ifblock=receivedB2A(i*640+[1:640],indA);
ifblock(1:128,:)=[];
fblock=fft(ifblock);
fblock(1,:)=[];
fblock(151:360,:)=[];
Db2a_R((Niter-1)*120+i+1,:)=fblock.';
endfor
HB2A=conj(repmat(Db2a_T,Niter,1)).*repmat(Db2a_R,1,Nantb);
phasesB2A=unwrap(angle(HB2A));
if(mean(var(phasesB2A))>0.5)
disp("The phases of your estimates are a bit high (larger than 0.5 rad.), something is wrong.");
endif
chanestsB2A=reshape(diag(repmat(Db2a_T,Niter,1)'*repmat(Db2a_R,1,Nantb)/(Niter*120)),301,Nantb);
# %% -- Some plotting code -- %% (you can uncomment what you see fit)
# clf
# figure(1)
# for i=1:4
# subplot(220+i);plot(20*log10(abs(fftshift(fft(receivedA2B(:,i))))));
# endfor
# figure(2)
# t=[0:512-1]/512*1e-2;
# plot(t,abs(tchanests))
# xlabel('time')
# ylabel('|h|')
# figure(3)
# % wndw = 50;
# % for i=1:5:Nantb*301 %# sliding window size
# % phamean = filter(ones(wndw,1)/wndw, 1, phases(:,i)); %# moving average
# % plot(phamean(wndw:end),'LineWidth',2);
# % title(['subcarrier ' num2str(i)]);
# % xlabel('time')
# % ylabel('phase')
# % ylim([-pi pi])
# % drawnow;
# % pause(0.1)
# % endfor
# phavar=var(phases);
# plotphavar=[];
# for i=0:Nantb-1
# plotphavar=[plotphavar; phavar([1:301]+i*301)];
# endfor
# plot([1:150 362:512],plotphavar,'o');
# %ylim([0 pi])
# xlabel('subcarrier')
# ylabel('phase variance')
# figure(4)
# plot(20*log10(abs(fchanests))), ylim([40 100])
# %end
# fprintf(' done\n')
# for i=0:(Nantb-1)
# fchanests(:,i+1)=[0; chanests(301*i+[1:150]); zeros(210,1); chanests(301*i+[151:301])];
# endfor
# tchanests=ifft(fchanests);
else
if(!LSBSWITCH_FLAG) error("You have to set the LSB switch flag (LSBSWITCH_FLAG) in initparams.m.\n")
else error("You have to run init.params.m first!")
endif
endif
# % Author: Mirsad Cirkic
# % Organisation: Eurecom (and Linkoping University)
# % E-mail: mirsad.cirkic@liu.se
if(paramsinitialized && ~LSBSWITCH_FLAG)
disp(["\n\n------------\nThis code is, so far, only written for single runs. Multiple " ...
"runs will overwrite the previous measurement data, i.e., the " ...
"data structures are not defined for multiple runs. You will need to " ...
"add code in order to save the intermediate measurements and the " ...
"corresponding timestamps.\n------------"])
N=76800;
M=4;
signalA2B=zeros(N,4);
signalB2A=zeros(N,4);
indA=find(active_rfA==1);
indB=find(active_rfB==1);
Nanta=length(indA);
Nantb=length(indB);
#%i f(Nanta!=1) error("Node A can only have one antenna active\n"); endif
Niter=1;
if(Niter!=1) error("We should only use one get_frame at each \
run.\n");
endif
# %% ------- Prepare the signals for both A2B and B2A ------- %%
signalA2B=zeros(N,4);
signalB2A=zeros(N,4);
ia=1; ib=1;
Db2a_T=[];
for i=1:4
if(indA(ia)==i)
[Da2b_T, tmps]=generqpskseq(N,M,amp);
%tmps(1024:2048) = 0;
signalA2B(:,i)=tmps;
if(length(indA)> ia) ia=ia+1; endif
endif
if(indB(ib)==i)
% This part could be improved by creating fully orthogonal sequences
[tmpd, tmps]=generqpskseq(N,M,amp);
%tmps(1024:2048) = 0;
signalB2A(:,i)=tmps;
Db2a_T=[Db2a_T tmpd];
if(length(indB)> ib) ib=ib+1; endif
endif
endfor
# %% ------- Node B to A transmission ------- %%
rf_mode_current = rf_mode + (DMAMODE_TX+TXEN)*active_rfB +(DMAMODE_RX+RXEN)*active_rfA;
oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode_current,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode);
oarf_send_frame(card,signalB2A,n_bit);
#% keyboard
receivedB2A=oarf_get_frame(card);
oarf_stop(card);
# %% ------- Node A to B transmission ------- %%
rf_mode_current = rf_mode + (DMAMODE_TX+TXEN)*active_rfA +(DMAMODE_RX+RXEN)*active_rfB;
oarf_config_exmimo(card, freq_rx,freq_tx,tdd_config,syncmode,rx_gain,tx_gain,eNB_flag,rf_mode_current,rf_rxdc,rf_local,rf_vcocal,rffe_rxg_low,rffe_rxg_final,rffe_band,autocal_mode);
oarf_send_frame(card,signalA2B,n_bit);
%keyboard
receivedA2B=oarf_get_frame(card);
oarf_stop(card);
# %% ------- Do the A to B channel estimation ------- %%
Da2b_R=zeros(Niter*120,Nantb*301);
for i=0:119;
ifblock=receivedA2B(i*640+[1:640],indB);
ifblock(1:128,:)=[];
fblock=fft(ifblock);
fblock(1,:)=[];
fblock(151:360,:)=[];
Da2b_R((Niter-1)*120+i+1,:)=vec(fblock);
endfor
HA2B=repmat(conj(Da2b_T),1,Nantb).*Da2b_R;
phasesA2B=unwrap(angle(HA2B));
if(mean(var(phasesA2B))>0.5)
disp("The phases of your estimates are a bit high (larger than 0.5 rad.), something is wrong.");
endif
chanestsA2B=reshape(diag(repmat(Da2b_T,1,Nantb)'*Da2b_R)/size(Da2b_T,1),301,Nantb);
fchanestsA2B=zeros(512,Nantb);
for i=1:Nantb
fchanestsA2B(:,i)=[0; chanestsA2B([1:150],i); zeros(210,1); chanestsA2B(151:301,i)];
endfor
tchanestsA2B=ifft(fchanestsA2B);
%% ------- Do the B to A channel estimation ------- %%
Db2a_R=zeros(Niter*120,Nanta*301);
for i=0:119;
ifblock=receivedB2A(i*640+[1:640],indA);
ifblock(1:128,:)=[];
fblock=fft(ifblock);
fblock(1,:)=[];
fblock(151:360,:)=[];
Db2a_R((Niter-1)*120+i+1,:)=fblock.';
Db2a_R((Niter-1)*120+i+1,:)=fblock.';
endfor
HB2A=conj(repmat(Db2a_T,Niter,1)).*repmat(Db2a_R,1,Nantb);
phasesB2A=unwrap(angle(HB2A));
if(mean(var(phasesB2A))>0.5)
disp("The phases of your estimates are a bit high (larger than 0.5 rad.), something is wrong.");
endif
chanestsB2A=reshape(diag(repmat(Db2a_T,Niter,1)'*repmat(Db2a_R,1,Nantb)/(Niter*120)),301,Nantb);
# %% -- Some plotting code -- %% (you can uncomment what you see fit)
# clf
# figure(1)
# for i=1:4
# subplot(220+i);plot(20*log10(abs(fftshift(fft(receivedA2B(:,i))))));
# endfor
# figure(2)
# t=[0:512-1]/512*1e-2;
# plot(t,abs(tchanests))
# xlabel('time')
# ylabel('|h|')
# figure(3)
# % wndw = 50;
# % for i=1:5:Nantb*301 %# sliding window size
# % phamean = filter(ones(wndw,1)/wndw, 1, phases(:,i)); %# moving average
# % plot(phamean(wndw:end),'LineWidth',2);
# % title(['subcarrier ' num2str(i)]);
# % xlabel('time')
# % ylabel('phase')
# % ylim([-pi pi])
# % drawnow;
# % pause(0.1)
# % endfor
# phavar=var(phases);
# plotphavar=[];
# for i=0:Nantb-1
# plotphavar=[plotphavar; phavar([1:301]+i*301)];
# endfor
# plot([1:150 362:512],plotphavar,'o');
# %ylim([0 pi])
# xlabel('subcarrier')
# ylabel('phase variance')
# figure(4)
# plot(20*log10(abs(fchanests))), ylim([40 100])
# %end
# fprintf(' done\n')
# for i=0:(Nantb-1)
# fchanests(:,i+1)=[0; chanests(301*i+[1:150]); zeros(210,1); chanests(301*i+[151:301])];
# endfor
# tchanests=ifft(fchanests);
else
if(LSBSWITCH_FLAG) error("You have to unset the LSB switch flag (LSBSWITCH_FLAG) in initparams.m.\n")
else error("You have to run init.params.m first!")
endif
endif
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