diff --git a/targets/PROJECTS/TDDREC/f_ch_est.m b/targets/PROJECTS/TDDREC/f_ch_est.m
new file mode 100644
index 0000000000000000000000000000000000000000..f0a127c72d2a8214961f702c034e5b9b651687bd
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_ch_est.m
@@ -0,0 +1,46 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/f_ofdm_mod.m b/targets/PROJECTS/TDDREC/f_ofdm_mod.m
new file mode 100644
index 0000000000000000000000000000000000000000..eba1952c221dbaafc9bf0441eaa66760996291f0
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_ofdm_mod.m
@@ -0,0 +1,18 @@
+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
diff --git a/targets/PROJECTS/TDDREC/f_ofdm_rx.m b/targets/PROJECTS/TDDREC/f_ofdm_rx.m
new file mode 100644
index 0000000000000000000000000000000000000000..a85fe3a8df1c7770b3b3d6c7933743d7e5599ab0
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_ofdm_rx.m
@@ -0,0 +1,46 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/f_ofdm_tx.m b/targets/PROJECTS/TDDREC/f_ofdm_tx.m
new file mode 100644
index 0000000000000000000000000000000000000000..646a03516f75d2e14fe9e514ce3d90d83d80c589
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_ofdm_tx.m
@@ -0,0 +1,59 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/f_tls_ap.m b/targets/PROJECTS/TDDREC/f_tls_ap.m
new file mode 100644
index 0000000000000000000000000000000000000000..5b4c5410be22a79ea549d7e8fbe8783ce17548b0
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_tls_ap.m
@@ -0,0 +1,53 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/f_tls_diag.m b/targets/PROJECTS/TDDREC/f_tls_diag.m
new file mode 100644
index 0000000000000000000000000000000000000000..cde8a577a5951dedd3cf082f5a10815e19dadf54
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_tls_diag.m
@@ -0,0 +1,46 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/f_tls_svd.m b/targets/PROJECTS/TDDREC/f_tls_svd.m
new file mode 100644
index 0000000000000000000000000000000000000000..c5604a6eefd4ab749d7b6c8111dd66e58840826a
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/f_tls_svd.m
@@ -0,0 +1,49 @@
+%
+% 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.
+% 2283–2302, 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
diff --git a/targets/PROJECTS/TDDREC/s_beamforming.m b/targets/PROJECTS/TDDREC/s_beamforming.m
new file mode 100644
index 0000000000000000000000000000000000000000..4005e657c079ae079ef506aceaae82658b440dca
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/s_beamforming.m
@@ -0,0 +1,161 @@
+%
+% 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])
diff --git a/targets/PROJECTS/TDDREC/s_calib.m b/targets/PROJECTS/TDDREC/s_calib.m
new file mode 100644
index 0000000000000000000000000000000000000000..bfcbb0c7b82ae88cd263bd0769061a4ee2b41efa
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/s_calib.m
@@ -0,0 +1,127 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/s_draw_F.m b/targets/PROJECTS/TDDREC/s_draw_F.m
new file mode 100644
index 0000000000000000000000000000000000000000..7071416ea2a5284b34466c227cdfa5e5d0d4ad78
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/s_draw_F.m
@@ -0,0 +1,115 @@
+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')
+
diff --git a/targets/PROJECTS/TDDREC/s_init_params.m b/targets/PROJECTS/TDDREC/s_init_params.m
new file mode 100644
index 0000000000000000000000000000000000000000..9683ed2d1747703901fba28e828c6f1b5a6cd793
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/s_init_params.m
@@ -0,0 +1,55 @@
+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;
+
diff --git a/targets/PROJECTS/TDDREC/s_run_meas.m b/targets/PROJECTS/TDDREC/s_run_meas.m
new file mode 100644
index 0000000000000000000000000000000000000000..67b3ab1f03e116d1a99be4e84b96550e80740197
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/s_run_meas.m
@@ -0,0 +1,139 @@
+%
+% 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
diff --git a/targets/PROJECTS/TDDREC/alterproj.m b/targets/PROJECTS/TDDREC/v0/alterproj.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/alterproj.m
rename to targets/PROJECTS/TDDREC/v0/alterproj.m
diff --git a/targets/PROJECTS/TDDREC/beamforming.m b/targets/PROJECTS/TDDREC/v0/beamforming.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/beamforming.m
rename to targets/PROJECTS/TDDREC/v0/beamforming.m
diff --git a/targets/PROJECTS/TDDREC/calibration.m b/targets/PROJECTS/TDDREC/v0/calibration.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/calibration.m
rename to targets/PROJECTS/TDDREC/v0/calibration.m
diff --git a/targets/PROJECTS/TDDREC/genorthqpskseq.m b/targets/PROJECTS/TDDREC/v0/genorthqpskseq.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/genorthqpskseq.m
rename to targets/PROJECTS/TDDREC/v0/genorthqpskseq.m
diff --git a/targets/PROJECTS/TDDREC/genrandpskseq.m b/targets/PROJECTS/TDDREC/v0/genrandpskseq.m
similarity index 99%
rename from targets/PROJECTS/TDDREC/genrandpskseq.m
rename to targets/PROJECTS/TDDREC/v0/genrandpskseq.m
index 2b23059fa6777ee234d7d5c2e30bc120e7079e24..141beabeae98f076031fd0dade031fd2774444b8 100644
--- a/targets/PROJECTS/TDDREC/genrandpskseq.m
+++ b/targets/PROJECTS/TDDREC/v0/genrandpskseq.m
@@ -26,4 +26,4 @@ for i=0:(N/640-1)
   s([1:640]+i*640)=floor(amp*block);
 end
 
-end
\ No newline at end of file
+end
diff --git a/targets/PROJECTS/TDDREC/initparams.m b/targets/PROJECTS/TDDREC/v0/initparams.m
similarity index 90%
rename from targets/PROJECTS/TDDREC/initparams.m
rename to targets/PROJECTS/TDDREC/v0/initparams.m
index 05e5674ec2003eedf4ea023e7e601a86f8cf8c30..5bfae4fa68be99ff23c186f88775e423bdecc4c9 100644
--- a/targets/PROJECTS/TDDREC/initparams.m
+++ b/targets/PROJECTS/TDDREC/v0/initparams.m
@@ -15,8 +15,9 @@ eNB_flag = 0;
 card = 0;
 Ntrx=4;
 dual_tx=0;
-active_rfA=[0 0 1 0];
-active_rfB=[1 1 0 0];
+active_rfA=[0 0 0 1];
+% active_rfB=[1 1 1 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
@@ -26,6 +27,7 @@ fc  = 1912600000; %1907600000;
 %fc = 859.5e6;
 
 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=(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
@@ -47,11 +49,12 @@ freq_tx = freq_rx; %+1.92e6;
 
 oarf_stop(card);
 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
 amp = pow2(12.5)-1;
 n_bit = 16;
 
+%chanest_full = 1;
 chanest_full = 1;
 
 paramsinitialized=true;
diff --git a/targets/PROJECTS/TDDREC/v0/initparams2.m b/targets/PROJECTS/TDDREC/v0/initparams2.m
new file mode 100644
index 0000000000000000000000000000000000000000..1d044033321c6f8755403ecddf023a3a58d3f13f
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/v0/initparams2.m
@@ -0,0 +1,49 @@
+# % 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;
diff --git a/targets/PROJECTS/TDDREC/ofdm_mod.m b/targets/PROJECTS/TDDREC/v0/ofdm_mod.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/ofdm_mod.m
rename to targets/PROJECTS/TDDREC/v0/ofdm_mod.m
diff --git a/targets/PROJECTS/TDDREC/v0/readme.txt b/targets/PROJECTS/TDDREC/v0/readme.txt
new file mode 100644
index 0000000000000000000000000000000000000000..17000d88128e86a8645676136c675ae07e6a744d
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/v0/readme.txt
@@ -0,0 +1 @@
+old project by Misard and Florian
diff --git a/targets/PROJECTS/TDDREC/runmeas_full_duplex.m b/targets/PROJECTS/TDDREC/v0/runmeas_full_duplex.m
similarity index 97%
rename from targets/PROJECTS/TDDREC/runmeas_full_duplex.m
rename to targets/PROJECTS/TDDREC/v0/runmeas_full_duplex.m
index b1348f7238b0129cbb59efeb08ae924d3f0cfa5b..ebcd375b39e5b6b6ad02150c960f52cba394e23e 100644
--- a/targets/PROJECTS/TDDREC/runmeas_full_duplex.m
+++ b/targets/PROJECTS/TDDREC/v0/runmeas_full_duplex.m
@@ -21,7 +21,7 @@ if(paramsinitialized)
     if(Niter~=1)
         error('We should only use one get_frame at each run.\n');
     end
-    Nmeas = 10;
+    Nmeas = 5;
     
     %% ------- Prepare the signals for A2B ---------- %%
     signalA2B=zeros(N,4);
@@ -75,19 +75,18 @@ if(paramsinitialized)
     for meas=1:Nmeas
         %% ------- Node A to B transmission ------- %%
         oarf_send_frame(card,signalA2B,n_bit);
-        %keyboard
-        sleep(0.01);
         receivedA2B=oarf_get_frame(card);
         %oarf_stop(card); %not good, since it does a reset
         sleep(0.01);
         
         %%----------Node B to A transmission---------%%
         oarf_send_frame(card,signalB2A,n_bit);
-        %keyboard
-        sleep(0.01);
+
         receivedB2A=oarf_get_frame(card);
         %oarf_stop(card); %not good, since it does a reset
         
+        %keyboard;
+
         %% ------- Do the A to B channel estimation ------- %%
         for i=0:119;
             ifblock=receivedA2B(i*640+[1:640],indB);
@@ -158,7 +157,7 @@ if(paramsinitialized)
       fblock=fft(ifblock);
       fblock(1,:)=[];
       fblock(151:360,:)=[];
-      noise_f(i+1,:,:)=fblock;
+      %noise_f(i+1,:,:)=fblock;
     end
 
 
@@ -177,8 +176,9 @@ if(paramsinitialized)
     end
     
     figure(2)
-    t=[0:512-1]/512*1e-2;
-    plot(t,20*log10(abs(tchanests)))
+    %t=[0:512-1]/512*1e-2;
+    %plot(t,20*log10(abs(tchanests)))
+    plot(20*log10(abs(tchanests)))
     xlabel('time')
     ylabel('|h|')
     if Nantb==3
diff --git a/targets/PROJECTS/TDDREC/runmeas_long_chanest.m b/targets/PROJECTS/TDDREC/v0/runmeas_long_chanest.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/runmeas_long_chanest.m
rename to targets/PROJECTS/TDDREC/v0/runmeas_long_chanest.m
diff --git a/targets/PROJECTS/TDDREC/runmeas_mimo.m b/targets/PROJECTS/TDDREC/v0/runmeas_mimo.m
similarity index 100%
rename from targets/PROJECTS/TDDREC/runmeas_mimo.m
rename to targets/PROJECTS/TDDREC/v0/runmeas_mimo.m
diff --git a/targets/PROJECTS/TDDREC/v0/runmeas_wduplex.m b/targets/PROJECTS/TDDREC/v0/runmeas_wduplex.m
new file mode 100644
index 0000000000000000000000000000000000000000..25310a814fa952b5868155804970f3d68519c9a3
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/v0/runmeas_wduplex.m
@@ -0,0 +1,155 @@
+# % 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
+	   	
diff --git a/targets/PROJECTS/TDDREC/v0/runmeas_woduplex.m b/targets/PROJECTS/TDDREC/v0/runmeas_woduplex.m
new file mode 100644
index 0000000000000000000000000000000000000000..164d24241dba55c01e63b773b89639bd0b8cff01
--- /dev/null
+++ b/targets/PROJECTS/TDDREC/v0/runmeas_woduplex.m
@@ -0,0 +1,157 @@
+# % 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