ChannelEstimator

 Estimates the channel impulse response relying on pilots distributed
 among the subcarriers.

 function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para)
 function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para, options)
 function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para, options, causality)

 The function works for SISO and MIMO systems.

 Input arguments:

   z: demodulated symbols. Size: matrix [2*Para.M, ~] if Para.N_R == 1,
   multidimensional array [para.N_R, 2*Para.M,~] if Para.N_R > 1.

   L: channel impulse response length.

   d_TR: vector [2*Para.M, 1] of pilot symbols.

   Para: structure containing the modulation parameters.

   options: optional, 'ScatteredPilots' or 'Preamble' (by default).

   causality: optional (default value set to 1), estimated a causal or non
   causal channel impulse response. If set to 0, non causal channel
   estimated meaning that the channel impulse response is centered in zero
   and has non zero coefficients for taps before 0. Then, C_est has a size
   para.nSubcarriers and end coefficients contains negative taps.

 Outputs arguments:

   C_est: channel impulse response. Multidimensional array [Para.N_R, Para.N_T, ~].

   N0_dB_est: estimated noise variance.

 Reference: [1] F. Rottenberg, F. Horlin, E. Kofidis and J. Louveaux,
 "Generalized optimal pilot allocation for channel estimation in
 multicarrier systems." In IEEE 17th International Workshop on Signal
 Processing Advances in Wireless Communications (SPAWC), 2016, Edinburgh.

0001 function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para, options, causality)
0002 % Estimates the channel impulse response relying on pilots distributed
0003 % among the subcarriers.
0004 %
0005 % function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para)
0006 % function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para, options)
0007 % function [ C_est, N0_dB_est ] = ChannelEstimator(z, L, d_TR, Para, options, causality)
0008 %
0009 % The function works for SISO and MIMO systems.
0010 %
0011 % Input arguments:
0012 %
0013 %   z: demodulated symbols. Size: matrix [2*Para.M, ~] if Para.N_R == 1,
0014 %   multidimensional array [para.N_R, 2*Para.M,~] if Para.N_R > 1.
0015 %
0016 %   L: channel impulse response length.
0017 %
0018 %   d_TR: vector [2*Para.M, 1] of pilot symbols.
0019 %
0020 %   Para: structure containing the modulation parameters.
0021 %
0022 %   options: optional, 'ScatteredPilots' or 'Preamble' (by default).
0023 %
0024 %   causality: optional (default value set to 1), estimated a causal or non
0025 %   causal channel impulse response. If set to 0, non causal channel
0026 %   estimated meaning that the channel impulse response is centered in zero
0027 %   and has non zero coefficients for taps before 0. Then, C_est has a size
0028 %   para.nSubcarriers and end coefficients contains negative taps.
0029 %
0030 % Outputs arguments:
0031 %
0032 %   C_est: channel impulse response. Multidimensional array [Para.N_R, Para.N_T, ~].
0033 %
0034 %   N0_dB_est: estimated noise variance.
0035 %
0036 % Reference: [1] F. Rottenberg, F. Horlin, E. Kofidis and J. Louveaux,
0037 % "Generalized optimal pilot allocation for channel estimation in
0038 % multicarrier systems." In IEEE 17th International Workshop on Signal
0039 % Processing Advances in Wireless Communications (SPAWC), 2016, Edinburgh.
0040 
0041 % This file is part of WaveComBox: www.wavecombox.com and is distributed under
0042 % the terms of the MIT license. See accompanying LICENSE file.
0043 % Original author: François Rottenberg, May 4, 2018.
0044 % Contributors:
0045 % Change log:
0046 
0047 if exist('options','var') == 0
0048     options='Preamble';
0049 end
0050 switch options
0051     case 'ScatteredPilots'
0052         PilotSubcarriers=Para.ScatteredPilotSubcarriers;
0053     case 'Preamble'
0054         PilotSubcarriers=Para.PreamblePilotSubcarriers;
0055     otherwise
0056         error('Options type not defined for this function')
0057 end
0058 
0059 switch Para.Modulation
0060     case 'FBMC-OQAM'
0061         factor=2;
0062     case 'CP-OFDM'
0063         factor=1;
0064     otherwise
0065         error('Modulation not implemented')
0066 end
0067 
0068 if exist('causality','var') == 0
0069     causality=1;
0070 end
0071 switch causality
0072     case 1 
0073         Sigma=[eye(L); zeros(Para.nSubcarriers-L,L)];
0074         C_est=zeros(Para.N_R,Para.N_T,L);
0075     case 0
0076         Sigma=[eye(L), zeros(L,L);zeros(Para.nSubcarriers-2*L,2*L);zeros(L,L) eye(L)];
0077         C_est=zeros(Para.N_R,Para.N_T,2*L);
0078     otherwise
0079         error('causality argument should be set to zero or one')
0080 end
0081 
0082 if Para.N_R==1
0083     H_pilot_subcarriers=zeros(Para.nSubcarriers,1);
0084     H_pilot_subcarriers(PilotSubcarriers)=z(PilotSubcarriers,1)./d_TR(PilotSubcarriers);
0085     F=exp(-1j*2*pi/Para.nSubcarriers*(0:Para.nSubcarriers-1)'*(0:Para.nSubcarriers-1));
0086     v=zeros(Para.nSubcarriers,1);
0087     v(PilotSubcarriers)=1;
0088     S=diag(v);
0089     Gamma=S*F*Sigma;
0090     error_temp=zeros(Para.N_R,Para.N_T);
0091     for index_N_T=1:Para.N_T
0092         for index_N_R=1:Para.N_R
0093             C_est(index_N_R,index_N_T,:)=Gamma\H_pilot_subcarriers;
0094             temp_reconstruct=Gamma*squeeze(C_est(index_N_R,index_N_T,:));
0095             error_temp(index_N_R,index_N_T)=mean(abs( temp_reconstruct(PilotSubcarriers)-squeeze(H_pilot_subcarriers(PilotSubcarriers))).^2);
0096         end
0097     end
0098     N0_dB_est=10*log10(mean(mean(error_temp)));
0099 else
0100     F_pilot=1/sqrt(Para.N_T)*exp(-1j*2*pi/Para.N_T*(0:Para.N_T-1)'*(0:Para.N_T-1)  );
0101     H_pilot_subcarriers=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0102     for m=PilotSubcarriers
0103         H_pilot_subcarriers(:,:,m)= squeeze(z(:,m,1:factor:factor*Para.N_T))*F_pilot'./d_TR(m);
0104     end
0105     
0106     F=exp(-1j*2*pi/Para.nSubcarriers*(0:Para.nSubcarriers-1)'*(0:Para.nSubcarriers-1));
0107     v=zeros(Para.nSubcarriers,1);
0108     v(PilotSubcarriers)=1;
0109     S=diag(v);
0110     Gamma=S*F*Sigma;
0111     error_temp=zeros(Para.N_R,Para.N_T);
0112     for index_N_T=1:Para.N_T
0113         for index_N_R=1:Para.N_R
0114             C_est(index_N_R,index_N_T,:)=Gamma\squeeze(H_pilot_subcarriers(index_N_R,index_N_T,:) );
0115             temp_reconstruct=Gamma*squeeze(C_est(index_N_R,index_N_T,:));
0116             error_temp(index_N_R,index_N_T)=mean(abs( temp_reconstruct(PilotSubcarriers)-squeeze(H_pilot_subcarriers(index_N_R,index_N_T,PilotSubcarriers))).^2);
0117         end
0118     end
0119     N0_dB_est=10*log10(mean(mean(error_temp)));
0120 end
0121 
0122 if causality==0
0123     temp=C_est;
0124     C_est=zeros(Para.N_R,Para.N_T,Para.nSubcarriers);
0125     for index_N_T=1:Para.N_T
0126         for index_N_R=1:Para.N_R
0127             C_est(index_N_R,index_N_T,:)=Sigma*squeeze(temp(index_N_R,index_N_T,:));
0128         end
0129     end
0130 end
0131 end
0132