add synchro fucntions for NPSS and NSSS signal for NB-IoT

openair1/PHY/LTE_REFSIG/first_synchro_NB_IoT.m

+function [ theta_estim, estim_CFO ] = Fc_first_synchro( observation, L_frame, L_sub_frame, FFT_size, L_symbol, N_subframe_observation, L_CP, SNR, type_first_estim )
+% This function performs the estimation of the beginning of symbols as well
+% as the estimation of CFO. It allows for a coarse synchronization 
+
+gamma = zeros(1,L_frame);
+epsilon = zeros(1,L_frame); 
+
+for n = 1 : 1 : length(gamma)
+     
+    gamma(n) = sum(observation(n:n+L_CP-1).*conj(observation(n+FFT_size:n+FFT_size+L_CP-1)));
+    epsilon(n) = sum(abs(observation(n:n+L_CP-1)).^2 + abs(observation(n+FFT_size:n+FFT_size+L_CP-1)).^2);
+end
+
+rho = 10^(SNR/20)/(10^(SNR/20)+1);
+theta = 2*abs(gamma)-rho*(epsilon); 
+
+% Estimation of the symbol start and the corresponding CFO
+
+theta_reshape = reshape(theta,L_symbol,N_subframe_observation*L_sub_frame); 
+% gamma_reshape = reshape(gamma,L_symbol,N_subframe_observation*L_sub_frame); % useful for estimation of CFO
+[~,index_max] = max(theta_reshape); % where theta is max symbol by symbol
+
+switch type_first_estim
+    case 1
+        %estimation by mean
+        theta_estim = sum(index_max)/length(index_max); % estimation by mean
+        estim_CFO = -1/(2*pi)*atan(imag(gamma(round(theta_estim)))/real(gamma(round(theta_estim)))); 
+    case 2
+        %estimation by majority
+        counter_index = zeros(1,L_symbol);  
+        for k = 1 : 1 : length(index_max)
+    
+        counter_index(index_max(k)) = counter_index(index_max(k)) + 1; % add the number of index_max
+    
+        end
+
+        [~,theta_estim] = max(counter_index); % get the max of index max -> theta estim
+%         index_index_max = find(index_max == theta_estim); 
+%         estim_CFO = -1/(2*pi)*atan(imag(gamma(round(theta_estim)))/real(gamma(round(theta_estim))));
+        estim_CFO_vec = -1/(2*pi)*atan(imag(gamma(round(theta_estim:L_symbol:end)))./real(gamma(round(theta_estim:L_symbol:end))));
+%         estim_CFO_vec2 = -1/(2*pi)*atan(imag(gamma_reshape(theta_estim,index_max(index_index_max)))./real(gamma_reshape(theta_estim,index_max(index_index_max))));
+        estim_CFO = sum(estim_CFO_vec)/length(estim_CFO_vec); 
+%         estim_CFO_2 = sum(estim_CFO_vec2)/length(estim_CFO_vec2); 
+    otherwise
+        print('error: type of estimation not defined')
+end
+

openair1/PHY/LTE_REFSIG/main_synchro_NPSS_NB_IoT.m

+clear all 
+close all 
+
+% description: test synchro using CP for time-freq. synchro
+% and ZC sequence for beginning of the radio frame estimation
+% date : 09/03/2017
+% author : Vincent Savaux, b<>com, Rennes, France
+% email: vincent.savaux@b-com.com
+
+% Parameters
+
+u = 5; % root of ZC sequence 
+size_RB = 12; % number of sub-carrier per RB 
+N_ZC = size_RB-1; 
+L_sub_frame = 14; 
+j = 1i; 
+CFO = 0.1; % normalized CFO 
+N_frames = 4; % at least 3 
+N_sub_frame = 10*N_frames; % how many simulated sub_frame you want 
+FFT_size = 128; 
+N_zeros = (FFT_size-size_RB)/2; % Number of zero subcarriers in upper and lower frequencies 
+L_CP = round(4.6875/(66.7)*FFT_size); % Number of samples of the CP 
+L_symbol = (FFT_size + L_CP); 
+L_frame = (FFT_size + L_CP)*L_sub_frame*10; 
+L_signal = (FFT_size + L_CP)*L_sub_frame*N_sub_frame; 
+normalized_time = 0 : 1 : L_signal-1; 
+SNR_start = 0; % in dB 
+SNR_end = 30; % in dB 
+N_subframe_observation = 10; % length of observation for syncronization 
+N_loop = 1000; % number of runs, for good statistics 
+type_first_estim = 2; % 1 -> estimation by mean, 2-> estimation by majority 
+
+matrix_error_theta_1 = zeros(N_loop,SNR_end-SNR_start+1); 
+matrix_error_theta_2 = zeros(N_loop,SNR_end-SNR_start+1);  
+matrix_error_angle_1 = zeros(N_loop,SNR_end-SNR_start+1); 
+matrix_error_angle_2 = zeros(N_loop,SNR_end-SNR_start+1); 
+matrix_error_angle_3 = zeros(N_loop,SNR_end-SNR_start+1); 
+matrix_error_BOF = zeros(N_loop,SNR_end-SNR_start+1); 
+
+for SNR = SNR_start : 2 : SNR_end
+for loop = 1 : N_loop     
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Creation of the signal
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ZC sequence in frequency domain
+
+vec_n = 0:N_ZC-1; 
+f_ZC_sequence = exp(-j*pi*u*vec_n.*(vec_n+1)/N_ZC); 
+f_NPSS_symbol = [f_ZC_sequence.';0]; % one NPSS symbol 
+f_NPSS_frame = [zeros(size_RB,3),kron(ones(1,L_sub_frame-3),f_NPSS_symbol)]; 
+
+% OFDM sub_frame in frequency domain -> modulation : QPSK
+%random QPSK elements: 
+f_OFDM_frames = (2*randi([0,1],size_RB,L_sub_frame*N_sub_frame)-1) + j*(2*randi([0,1],size_RB,L_sub_frame*N_sub_frame)-1); 
+
+%replace the k*6th subframes by f_NPSS_frame:
+f_LTE_frames = f_OFDM_frames; 
+for k = 0 : N_frames-1
+   
+    N_index = k*10*L_sub_frame + 85; 
+    f_LTE_frames(:,N_index:N_index+13) = f_NPSS_frame; 
+    
+end
+
+% IFFT: get frames in time domain (Parralel representation)
+
+f_zero_carriers = zeros(N_zeros,L_sub_frame*N_sub_frame); 
+f_ups_LTE_frames = [f_zero_carriers;f_LTE_frames;f_zero_carriers]; % add zero carriers up and down the symbols
+t_P_LTE_frames = ifft(f_ups_LTE_frames,FFT_size); 
+
+% Add CP (Parralel representation) 
+
+t_P_LTE_frames_CP = [t_P_LTE_frames(end-L_CP+1:end,:);t_P_LTE_frames]; 
+
+% Parralel to series conversion 
+
+t_S_LTE_frames = reshape(t_P_LTE_frames_CP,1,[]); 
+
+% Add a channel frequency offset (CFO)
+
+t_S_received_frames = t_S_LTE_frames.*exp(j*2*pi*CFO*normalized_time/FFT_size); 
+
+% Add noise 
+
+P_signal = sum(abs(t_S_received_frames).^2)/length(t_S_received_frames); 
+P_noise = P_signal*10^(-SNR/20); 
+
+init_noise = randn(size(t_S_received_frames));
+normalized_noise = init_noise/sqrt(sum(abs(init_noise).^2)/length(init_noise));
+noise = sqrt(P_noise)*normalized_noise;
+
+t_S_noisy_frames = t_S_received_frames + noise;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Time and frequency synchronization
+% The principle is based on the croos-correlation between the received
+% sequence and the transmitted SC sequence. 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Get an observation, the duration of which is one frame length. The
+% beginning of the stored samples is at 1 frame +- 0.5 frame
+
+index_start = L_frame + randi([-L_frame/2,L_frame/2],1);
+observation = t_S_noisy_frames(index_start:index_start+1.5*L_frame-1);
+f_oversampl_NPSS_symbol = [f_zero_carriers(:,1);f_NPSS_symbol;f_zero_carriers(:,1)];
+t_NPPS_unit = ifft(f_oversampl_NPSS_symbol,FFT_size); 
+t_NPPS_correl = [t_NPPS_unit(end-L_CP+1:end);t_NPPS_unit]; 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Estimation of start of the symbols and the CFO 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+[ theta_estim, estim_CFO ] = first_synchro( observation, L_frame, L_sub_frame, FFT_size, L_symbol, N_subframe_observation, L_CP, SNR, type_first_estim ); 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Time and frequency synchronization
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+new_index_start = index_start + theta_estim - 1; 
+new_vec_time = new_index_start-1 : 1 : new_index_start+1.5*L_frame-2; 
+new_observation = t_S_noisy_frames(new_index_start:new_index_start+1.5*L_frame-1).*exp(-j*2*pi*estim_CFO*new_vec_time/FFT_size); 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Second synchronization: beginning of frame (BOF)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+[ BOF ] = second_synchro( new_observation, f_NPSS_symbol , L_frame, L_symbol, FFT_size, L_CP, N_zeros ); 
+
+
+exact_index = L_symbol - rem(index_start,L_symbol) + 2; 
+[error_theta_1,ind_error1] = min([abs(exact_index - theta_estim-L_symbol), abs(exact_index - theta_estim),abs(exact_index - theta_estim+L_symbol)]); 
+% [error_theta_2,ind_error2] = min([abs(exact_index - theta_estim_2-137), abs(exact_index - theta_estim_2),abs(exact_index - theta_estim_2+137)]); 
+
+error_CFO_1 = CFO - estim_CFO; 
+% error_CFO_2 = CFO - estim_CFO_1; 
+% error_CFO_3 = CFO - estim_CFO_2; 
+
+vec_exact_index = 85 : 140 : N_frames*140;  
+estim_BOF = ceil((index_start-1)/L_symbol) + BOF ; 
+error_BOF = min(abs(estim_BOF-vec_exact_index)); 
+
+matrix_error_theta_1(loop,SNR-SNR_start+1) = error_theta_1; 
+% matrix_error_theta_2(loop,SNR-SNR_start+1) = error_theta_2;  
+matrix_error_angle_1(loop,SNR-SNR_start+1) = error_CFO_1; 
+% matrix_error_angle_2(loop,SNR-SNR_start+1) = error_CFO_2; 
+% matrix_error_angle_3(loop,SNR-SNR_start+1) = error_CFO_3;
+matrix_error_BOF(loop,SNR-SNR_start+1) = error_BOF;
+
+end
+end
+plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_theta_1(:,1:2:SNR_end+1)).^2)/N_loop))
+hold
+% plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_theta_2(:,1:2:SNR_end+1)).^2)/N_loop))
+figure
+plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_angle_1(:,1:2:SNR_end+1)).^2)/N_loop))
+hold
+% plot(SNR_start:SNR_end,sum(abs(matrix_error_angle_2(:,1:31)))/N_loop)
+% plot(SNR_start:SNR_end,sum(abs(matrix_error_angle_3(:,1:31)))/N_loop)
+% plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_angle_2(:,1:2:SNR_end+1)).^2)/N_loop),'s')
+% plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_angle_3(:,1:2:SNR_end+1)).^2)/N_loop),'d')
+figure 
+plot(SNR_start:2:SNR_end,sqrt(sum(abs(matrix_error_BOF(:,1:2:SNR_end+1)).^2)/N_loop))
+hold
+save
+

openair1/PHY/LTE_REFSIG/second_synchro_NB_IoT.m

+function [ BOF ] = Fc_second_synchro( new_observation, f_NPSS_symbol , L_frame, L_symbol, FFT_size, L_CP , N_zeros )
+%UNTITLED2 Summary of this function goes here
+%   Detailed explanation goes here
+%     new_obs_reshape = reshape(new_observation(1:L_frame), L_symbol, []); 
+    new_obs_reshape = reshape(new_observation, L_symbol, []); 
+    t_new_obs_CP_remov = new_obs_reshape(L_CP+1:end,:); 
+    f_new_symbols = fft(t_new_obs_CP_remov,FFT_size); 
+    
+    for n = 1 : length(f_new_symbols(1,:))
+        corr(:,n) = xcorr(f_new_symbols(N_zeros+1:N_zeros+12,n),f_NPSS_symbol); 
+        mean = sum(abs(corr(:,n))); 
+        mm(n) = sum((abs(corr(:,n))-mean).^2);
+    end
+    
+    for k = 1 : length(mm) - 13
+       min_var(k) = sum(mm(k:k+13));  
+    end
+    [~,BOF] = min(min_var);
+end
+

openair1/PHY/LTE_TRANSPORT/nsss_gen_NB_IoT.m

+clear all
+% nsss_gen / matlab
+% Copyright 2016 b<>com. All rights reserved.
+
+% description: generation of NSSS subframe
+% Reference: 3GPP TS36.211 release 13
+% author: Vincent Savaux, b<>com, Rennes, France
+% email: vincent.savaux@b-com.com
+
+% Input : \
+% Output : matrix NSSS_frame
+% Parameters
+% frame_number = 100;
+% cellID = 200;
+% % % Mapping results to estimated u-3
+
+
+SNR_start = -10; 
+SNR_end = 2; 
+vec_SNR = SNR_start : 2 : SNR_end; 
+N_loop = 40; 
+Proba_fail = zeros(1,length(vec_SNR)); 
+
+mat_bn = zeros(4,128); % mat_bn contains the 4 possible Hadamard sequences defined in the standard
+mat_bn(1,:) = ones(1,128);
+mat_bn(2,:) = [1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  1  ...
+    -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  ...
+    1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  1  -1  1  1  -1  1  -1  ...
+    -1  1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  -1  1  1  -1  -1  ...
+    1  1  -1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  ...
+    -1  1  1  -1  1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  ...
+    1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1];
+mat_bn(3,:) = [1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  1  ... 
+    -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  -1  1  1  ...
+    -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  -1  1  ...
+    1  -1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  -1  ...
+    1  1  -1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  ...
+    -1  1  1  -1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  ...
+    -1  1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  1];
+mat_bn(4,:) = [1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  1  ...
+    -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  -1  1  1  ...
+    -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  -1  1  ...
+    1  -1  -1  1  1  -1  1  -1  -1  1  -1  1  1  -1  1  -1  -1  1  1  ...
+    -1  -1  1  -1  1  1  -1  1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  ...
+    1  -1  -1  1  1  -1  -1  1  -1  1  1  -1  -1  1  1  -1  1  -1  -1  ...
+    1  -1  1  1  -1  1  -1  -1  1  1  -1  -1  1  -1  1  1  -1];
+mat_bn = [mat_bn,mat_bn(:,1:4)]; % see the definition of m in stadard
+
+mat_theta_f = zeros(4,132); % mat_bn contains the 4 possible phase sequences defined in the standard
+mat_theta_f(1,:) = ones(1,132); 
+mat_theta_f(2,:) = repmat([1,-j,-1,j],1,33);  
+mat_theta_f(3,:) = repmat([1,-1],1,66); 
+mat_theta_f(4,:) = repmat([1,j,-1,-j],1,33); 
+
+mat_16_theta = round(kron(mat_theta_f,ones(4,1))); % mat_bn contains the 4x4=16 possible pseudo-random sequences
+mat_16_bn = repmat(mat_bn,4,1); 
+mat_16 = mat_16_theta.*mat_16_bn;
+corresponding_values = zeros(16,2); % first column for q, second for theta_f
+corresponding_values(:,1) = repmat([0;1;2;3],4,1); % mapping column to q 
+corresponding_values(:,2) = kron([0;1;2;3],ones(4,1)); % mapping column to theta_f
+
+for k = 1 : length(vec_SNR) % loop on the SNR
+    N_fail = 0; 
+for loop = 1 : N_loop    
+SNR = vec_SNR(k); 
+frame_number = 2*randi([0,3],1);
+cellID = randi([0,503],1);    
+    
+% function NSSS_subframe = nsss_gen(frame_number,cellID)
+theta_f = 33/132*mod(frame_number/2,4); % as defined in stadard
+u = mod(cellID,126) + 3; % root of ZC sequence, defined in standard
+q = floor(cellID/126);
+size_RB = 12; % number of sub-carrier per RB
+N_ZC = 131; 
+L_sub_frame = 14; % number of OFDM symbols per subframe
+j = 1i; 
+vec_n = 0:N_ZC; 
+vec_n1 = mod(vec_n,131); 
+vec_bq = mat_bn(q+1,:); 
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Creation of the signal
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% ZC sequence in frequency domain
+
+ZC_sequence = exp(-j*pi*u*vec_n1.*(vec_n1+1)/N_ZC); 
+had_sequence = exp(-j*2*pi*theta_f*vec_n);
+vec_bq_had = vec_bq.*had_sequence; 
+P_noise = 10^(-SNR/10); % SNR in dB to noise power
+noise = sqrt(P_noise/2)*randn(1,132)+sqrt(P_noise/2)*j*randn(1,132); 
+vec_d = vec_bq.*had_sequence.*ZC_sequence + noise; 
+mat_NSSS = flipud(reshape(vec_d,size_RB,L_sub_frame-3));
+NSSS_subframe = [zeros(size_RB,3),mat_NSSS]; 
+
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Exhaustive cell ID research
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
+
+sequence_r = repmat(vec_d,16,1).*conj(mat_16); % this remove the phase component
+vec_u = 3 : 128; 
+mat_u = repmat(vec_u.',1,length(vec_n1));  
+mat_n1 = repmat(vec_n1,126,1); 
+sequence_ZC = exp(-j*pi*mat_u.*mat_n1.*(mat_n1+1)/N_ZC); 
+
+matrix_max_correl = zeros(126,16); % this will be filled by the maximum of correlation value 
+for s_ = 1 : 16 
+    seq_ref = sequence_r(s_,:); 
+    for u_ = 1 : 126  
+       correl = xcorr(seq_ref,sequence_ZC(u_,:)); 
+       [val_max,ind_max] = max(abs(correl)); 
+       matrix_max_correl(u_,s_) = val_max; 
+    end
+end
+
+max_correl = max(max(matrix_max_correl));  % get the max of all correlation values 
+index_max = find(matrix_max_correl==max_correl); 
+estim_u_ = mod(index_max,126)-1;
+index_column = (index_max-mod(index_max,126))/126+1; 
+estim_q_ = corresponding_values(index_column); 
+estim_cell_ID = q*126 + estim_u_; 
+
+if cellID ~= estim_cell_ID
+    N_fail = N_fail + 1; 
+end
+
+end
+Proba_fail(k) = N_fail/N_loop; 
+end
+
+plot(vec_SNR,Proba_fail)