% WEBDEMO .m-file, Institute of Telecommunications, University of Stuttgart, Germany
% Author: Nagarjuna Siddam

%In this file estimation of Rayleigh fading channel by 
%Linear Interpolation, Least Squares and Weiner LMMSE methods is implemented.

% EXAMPLE:  channel_estimation(5, 10, 100, 1.2e-3, 10, 100, 200e6, 1, 12345, 8, [0 1 0], 0, 'filename')

% PARAMETERS:
%   Es_No_dB             % signal to Noise ratio in dB
%   pilot_spacing        % spacing between the pillots
%   nb_pilots            % number of pilots. Total number of samples is calculated based on nb_pilots and pilot_spacing 
%   Ts,                  % duration of one discrete-time sample in sec. (set to 1 for normalization)
%   fD,                  % max. Doppler freq. shift in Hz       (note: only Ts*fD relevant)
%   speed_v_kmh,         % alternatively, Doppler shift can be defined through vehicle speed;  in km/h
%   fc,                  % carrier frequency in Hz
%   speed_en,            % 1=use speed input for computing Doppler frequency, not fD
%   LI_LMMSE             % Set 0 for LI and 1 for LMMSE
%   webdemo_en,          % if 0, then test-mode in Octave/Matlab commandline, no webserver running
%   filename             % output file for webdemo


% BEGIN, main function 
function channel_estimation(Es_No_dB, pilot_spacing, total_coeff, Ts, fD, speed_v_kmh, fc, speed_en , method_sel, snr_est_off, webdemo_en, filename)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% BEGIN, set up graphical output

	% Create an invisible figure.
	if (webdemo_en)
		fig = figure(1); set(fig, 'visible', 'off');
	else
		fig = figure(1); set(fig, 'visible', 'on');
	end
	
	landscape = '-S930,350';
	portrait  = '-S640,480';

	output_format = landscape;     % default
	
	set(0, "defaultlinelinewidth", 2);
	set(0, "defaulttextfontname", "Helvetica");
	set(0, "defaulttextfontsize", 8);
	
	if (webdemo_en)
		set(0, "defaultaxesfontsize", 8);
		set(0, "defaultaxesfontname", "Helvetica");
		set(0, "defaulttextfontsize", 8);
		set(0, "defaulttextfontname", "Helvetica");
		const_lw = 2;
		const_msize = 4;
	else
		set(0, "defaultaxesfontsize", 10);
		set(0, "defaultaxesfontname", "Helvetica");
		set(0, "defaulttextfontsize", 8);
		set(0, "defaulttextfontname", "Helvetica");
		const_lw = 2;
		const_msize = 2;
	endif
	
	% RGB color definitions, http://www.rapidtables.com/web/color/RGB_Color.htm
	gray 	   = [[0,0,0]; [32,32,32]; [64,64,64]; [96,96,96]; [128,128,128]; [160,160,160]; [192,192,192]; [224,224,224]; [255,255,255]]/255;
	violett    = [[51,0,25]; [102,0,51]; [153,0,76]; [204,0,102]; [255,0,127]; [255,51,153]; [255,102,178]; [255,153,204]; [255,204,229]]/255;
	magenta    = [[51,0,51]; [102,0,102]; [153,0,153]; [204,0,204]; [255,0,255]; [255,51,255]; [255,102,255]; [255,153,255]; [255,204,255]]/255;
	lilac      = [[25,0,51]; [51,0,102]; [76,0,153]; [102,0,204]; [127,0,255]; [153,51,255]; [178,102,255]; [204,153,255]; [229,204,255]]/255;
	blue       = [[0,0,51]; [0,0,102]; [0,0,153]; [0,0,204]; [0,0,255]; [51,51,255]; [102,102,255]; [153,153,255]; [204,204,255]]/255;
	seablue    = [[0,25,51]; [0,51,102]; [0,76,153]; [0,102,204]; [0,128,255]; [51,153,255]; [102,178,255]; [153,204,255]; [204,229,255]]/255;
	cyan       = [[0,51,51]; [0,102,102]; [0,153,153]; [0,204,204]; [0,255,255]; [51,255,255]; [102,255,255]; [153,255,255]; [204,255,255]]/255;
	seagreen   = [[0,51,25]; [0,102,51]; [0,153,76]; [0,204,102]; [0,255,128]; [51,255,153]; [102,255,178]; [153,255,204]; [204,255,229]]/255;
	green      = [[0,51,0]; [0,102,0]; [0,153,0]; [0,204,0]; [0,255,0]; [51,255,51]; [102,255,102]; [153,255,153]; [204,255,204]]/255;
	lightgreen = [[25,51,0]; [51,102,0]; [76,153,0]; [102,204,0]; [128,255,0]; [153,255,51]; [178,255,102]; [204,255,153]; [229,255,204]]/255;
	yellow     = [[51,51,0]; [102,102,0]; [153,153,0]; [204,204,0]; [255,255,0]; [255,255,51]; [255,255,102]; [255,255,153]; [255,255,204]]/255;
	amber      = [[51,25,0]; [102,51,0]; [153,76,0]; [204,102,0]; [255,128,0]; [255,153,51]; [255,178,102]; [255,204,153]; [255,229,204]]/255;
	red        = [[51,0,0]; [102,0,0]; [153,0,0]; [204,0,0]; [255,0,0]; [255,51,51]; [255,102,102]; [255,153,153]; [255,204,204]]/255;
  
  % END, set up graphical output
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  
	%calculate number of pilots
	nb_pilots = fix(total_coeff/ pilot_spacing);
  
    %calculating total number of samples
    total_coeff = nb_pilots*pilot_spacing-pilot_spacing+1;
    
    %generate the rayleigh fading coefficients
	t_offset = -12345; % time offset to improve randomness; if -1, then random offset
	N0 = 8; % nb of scatterers used to generate complex Rayleigh fading phasor
    rayleigh_coeff = rayleigh_fading_chest(total_coeff, Ts, fD, speed_v_kmh, fc, speed_en, t_offset, N0);

    %pilot symbol instances
    pilots_instances = 1:pilot_spacing:total_coeff;
	
	%generate the BPSK symbols NOTE: these are only pilot symbols
    x = rand(1,nb_pilots)>0.5;
    x = 2*x-1;

    %generate the AWGN noise
    n = 1/sqrt(2)*(randn(1,nb_pilots) + 1i*randn(1,nb_pilots))*10^(-Es_No_dB/20);

    %received pilot symbols
    recieved_pilots = rayleigh_coeff(pilots_instances).*x+n;

    %estimate the rayleigh coefficients for pilot symbols
    rayleigh_est_pilots = recieved_pilots./x;
    
    %**********************************************%
    %From here is the channel estimation and plotting
    hold on
	
	%define x-axis
    t = (0:1:total_coeff-1)*Ts;
    
    %plot pilot positions vs received pilot rayleigh coefficients
    %NOTE: These received rayleigh coefficeints have noise factor
    pilotx = t(pilots_instances);
    piloty = real(rayleigh_est_pilots);
    plot(pilotx, piloty,  '.', 'color', 'm', 'markersize',7);
    
    if (method_sel(1))
        
        %estimate the channel coefficients by linear interpolation
        est_li = interp1(pilots_instances, rayleigh_est_pilots, 1:total_coeff);
        
        %MSE calculation
        MSE_li = ((real(rayleigh_coeff - est_li))*(real(rayleigh_coeff - est_li))' + (imag(rayleigh_coeff - est_li))*(imag(rayleigh_coeff - est_li))')/total_coeff;
        MSE_li_db = 10*log10(MSE_li);
        
        %Plot Linear Interpolation
        txt1 = ['MSE = ', num2str(MSE_li_db, 3), 'dB]'];
        plot(t, real(est_li), 'linewidth', 5,  'color', red(7,:));
		txt_legend = 'LI Estimation' ;
		txt2 = 'Channel Estimation By Linear Interpolation [';
        title(strcat(txt2,txt1));
        
    elseif (method_sel(2))
        
        %estimate the channel coefficients by Least Squares
        est_ls = ls(pilots_instances, pilot_spacing, total_coeff, rayleigh_est_pilots, Ts);
    
        %MSE calculation
        MSE_ls = ((real(rayleigh_coeff - est_ls))*(real(rayleigh_coeff - est_ls))' + (imag(rayleigh_coeff - est_ls))*(imag(rayleigh_coeff - est_ls))')/total_coeff;
        MSE_ls_db = 10*log10(MSE_ls);
    
        %Plot Least Squares
        txt1 = ['MSE = ', num2str(MSE_ls_db, 3), 'dB]'];
        plot(t, real(est_ls), 'linewidth', 5, 'color', red(7,:));
		txt_legend = 'Least Squares Estimation' ;
		txt2 = 'Channel Estimation By Least Squares Interpolation [';
        title(strcat(txt2,txt1));
        
    elseif (method_sel(3))
        
        %calculate the maximum Doppler shift fD
        if (speed_en)
            fD = fc*(speed_v_kmh*1000/3600.0)/3e8;
        end
    
		Es_No_dB_off = Es_No_dB + snr_est_off;
        %estimate the channel coefficients by lmmse
        est_lmmse_real = lmmse(pilots_instances, total_coeff, real(rayleigh_est_pilots), Es_No_dB_off, fD, Ts);
		est_lmmse_imag = lmmse(pilots_instances, total_coeff, imag(rayleigh_est_pilots), Es_No_dB_off, fD, Ts);

        %MSE calculation
        MSE_lmmse = ((real(rayleigh_coeff) - est_lmmse_real)*(real(rayleigh_coeff) - est_lmmse_real)' + (imag(rayleigh_coeff) - est_lmmse_imag)*(imag(rayleigh_coeff) - est_lmmse_imag)')/total_coeff;
        MSE_lmmse_db = 10*log10(MSE_lmmse);
        
        %Plot LMMSE
        txt1 = ['MSE = ', num2str(MSE_lmmse_db, 3), 'dB]'];
        plot(t, est_lmmse_real, 'linewidth', 5, 'color', red(7,:));
		txt_legend = 'LMMSE Estimation' ;
		txt2 = 'Channel Estimation By LMMSE [';
        title(strcat(txt2,txt1));
        
    end
	
	%plot time vs original rayleigh coefficients
    y = real(rayleigh_coeff);
    plot(t, y, 'b');
	
	%insert legend
	h = legend('received noisy pilots', txt_legend, 'original Rayleigh coefficients');
    set(h, 'fontsize', 7);
	
	%set limits for axis
	xlim([0 t(end)]);
	ylim([min(piloty)-0.2 max(piloty)+1]);
	
    %Add a text block on the plot displaying the MSE in db
    %text(1.5,-2.5,txt);
    xlabel('time in sec.');
    ylabel('real part of Rayleigh coefficients');
	
	
	% Create an image file. This image is displayed on the webpage. Its size
	% must be 640x480.
	print(filename, '-dpng',output_format);
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN: Least Squares Estimation methods
function est_ls = ls(pilots_instances, pilot_spacing, total_coeff, est_pilots, Ts)

    %The window length is fixed
    window_length = 8;
    
    est_ls = zeros(1,total_coeff);
    est_ls(pilots_instances) = est_pilots;

    for  i = 1:length(pilots_instances)-1
        
        instances_vector = get_instances_vector(i, pilots_instances);
        pilot_vector = est_ls(instances_vector);
        poly_matrix = get_poly_matrix(instances_vector, Ts);
        weights_ls = poly_matrix \ pilot_vector.';
        
        for i_2 = 1:pilot_spacing
            coeff = pilot_spacing*(i-1) + i_2;
            est_ls(coeff) = (weights_ls).' * get_poly_vector(coeff*Ts);
        end
        
    end
    
end

%This function returns the sliding window of pilots
function instances_vector = get_instances_vector(i, pilots_instances)
    
    if i < 4
        instances_vector = pilots_instances(1:8);
    elseif i > length(pilots_instances)-4
        instances_vector = pilots_instances(length(pilots_instances)-7:end);
    else
        instances_vector = pilots_instances(i-3:i+4);
    end

end

%This function returns the matrix with the values of the polynomials
%evaluated
function poly_matrix = get_poly_matrix(instances_vector, Ts)
    
    %The matrix is of the size window length * polynomial degree+1
    poly_matrix = zeros(8,4);
    
    for i = 1:8
        poly_var = instances_vector(i)*Ts;
        poly_matrix(i,:) = (get_poly_vector(poly_var))';
    end
    
end

%This function returns a vector of polynomials
function poly_vector = get_poly_vector(poly_var)
    
    poly_vector = zeros(4,1);
    poly_vector(1) = 1;
    poly_vector(2) = poly_var;
    poly_vector(3) = 0.5*(3*poly_var^2 - 1);
    poly_vector(4) = 0.5*(5*poly_var^3 - 3*poly_var);
    
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN: LMMSE estimation
function est_lmmse = lmmse(pilots_instances, total_coeff, rayleigh_est_pilots, Es_No_dB_off, fD, Ts)
    est_lmmse = zeros(1,total_coeff);
    
    auto_cor_inv = calc_auto_cor(pilots_instances, Es_No_dB_off, fD, Ts );
    
    for i = 1:total_coeff
         cross_cor = besselj(0,2*pi*fD*(pilots_instances'-i)*Ts);
         weiner_coeff = auto_cor_inv*cross_cor;
         est_lmmse(i) = weiner_coeff'*rayleigh_est_pilots';
    end
end


function auto_cor_inv = calc_auto_cor(pilots_instances, Es_No_dB_off, fD, Ts )
    diff_fun = repmat(pilots_instances,length(pilots_instances),1) - repmat(pilots_instances',1,length(pilots_instances));
    auto_cor = besselj(0, 2*pi*fD*diff_fun*Ts);
    auto_cor_inv = inv(auto_cor + 10^(-Es_No_dB_off/20)*eye(length(pilots_instances)));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% BEGIN: Rayleigh fading generation.
%These functions are extracted from the webdemo: http://webdemo.inue.uni-stuttgart.de/webdemos/02_lectures/communication_3/rayleigh_fading/
function h_ray = rayleigh_fading_chest(nb_samples, Ts,fD,speed_v_kmh,fc, speed_en,  t_offset, N0)

  if (speed_en)
    fD = fc*(speed_v_kmh*1000/3600.0)/3e8;
  else % take fD from commandline
    speed_v_kmh = fD/fc/1000*3600*3e8;
  end

  % note: only Ts*fD relevant; Ts, fD separate only for labeling of plot axes
  TsfD = Ts*fD;
  
  if (t_offset<0)
    t_offset = floor(rand(1)*1e6);   % random time offset
  end

   h_ray = get_time_sequence(TsfD, nb_samples, N0, t_offset);

end


function y = get_time_sequence(TsfD, nb_samples, N0, t_offset)

  % some constants
  wm = 2*pi*TsfD;

  is_vec = zeros(1,nb_samples);
  qs_vec = zeros(1,nb_samples);

  t = t_offset;

  for n=1:nb_samples
    [is, qs] = get_coeff(N0, wm, t);

    is_vec(n) = is;
    qs_vec(n) = qs;

    t=t+1;    % discrete time index
  end

  % check power normalization
%  mean_is = sum(is_vec .* is_vec)/nb_samples;  % checked 2013/11/27, E[is^2]=0.5, ok
%  mean_iq = sum(qs_vec .* qs_vec)/nb_samples;  % checked 2013/11/27, E[qs^2]=0.5, ok

  y = is_vec + j*qs_vec;
end


function [is, qs] = get_coeff(N0, wm, t)

  is = 0;
  qs = 0;

  for n=1:N0
    beta_n = pi*n/N0;                   % modified Jakes model ('Dent model')
    tmp = cos(wm*t*cos(pi*(n-0.5)/2/N0 ));   % note:  wm = 2*pi*TsfD

    is = is + sin(beta_n)*tmp;
    qs = qs + cos(beta_n)*tmp;
  end

  is = is*sqrt(2/N0);
  qs = qs*sqrt(2/N0);

end