%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Timothy J. Peters
%   John Waterston
%
%   EE 360
%   Final Project
%
%   Blind adaptive closed loop power control
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


clear all
close all
warning off

avg = 5;

time_plot =0;
big_plot = 0;
little_plot = 0;

vel = [ 100 ];                        %  Vehicle velocity [km/h]
increment = [1];
fc = 2000;                                              %  Carrier frequency [fc]
Rs = 1500;                                              %  Power update rate [Hz]
fs = 4;                                                 %  Oversample ratio
rate = [8 4 2 1 0.5 0.25].*fs;
updates = 1400;                                          %  Number of updates calculated
frame = 700;                                            %  Updates per frame
pwr_target = 0;                                         %  Target power [dB]

N = Rs*fs;                                              %  Number of samples
time_step = 1/(N);                                      %  Time step [s]


for i = 1:avg
    %randn('state',i-1);
    
    for k = 1:length(rate)
        rate(k)
    %for k = 1:length(vel)
        v = vel;
     %   v = vel(k)
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  Generate Rayleigh channel 
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  Note:  Code to generate Rayleigh channel written by
        %  Mark A. Wickert at the University of Colorado
        %  in Colorado Springs in July, 1994
        %  and applied by John Waterston in his senior project
        %  entitled "Development and Implementation of an Adaptive
        %  Error Correction Coding Scheme for a Full Duplex Comms
        %  Channel."   
        
        [rw] = rayleigh_channel(v,fc,Rs,fs, N);
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        pwr = 10*log10(abs(rw).^2);                             %  Normalized power received at the base station [dB]
        pwr = pwr(1:updates*fs);                                %  Only look at the desired number of updates
        
        for j = 1:length(increment)
            inc = increment(j);
            dec = -increment(j);
            
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %  UMTS and IS95 power control algorithms
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            
            [pc1,error1] = UMTS(pwr,inc,dec,pwr_target,rate(k));    %  Calulate power control, error, and std of error

            %[pc1,error1] = UMTS(pwr,inc,dec,pwr_target,fs);          %  Calulate power control, error, and std of error
            %[pc2,error2] = IS95(pwr,inc,dec,pwr_target,fs);          %  Calulate power control, error, and std of error
            
            error1_avg(i,j,k,:) = error1;
            %error2_avg(i,j,k,:) = error2;
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        end
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %  Adaptive power control algorithm
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        
        [pc3,error3] = adaptive(pwr,pwr_target,rate(k));                  %  Calulate power control, error, and std of error
        %[pc3,error3] = adaptive(pwr,pwr_target,fs);                        %  Calulate power control, error, and std of error
       
        error3_avg(i,k,:) = error3;
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    end
end   



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Average CDMA_pc data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for j = 1:length(increment)
        for k = 1:length(rate)
        %for k = 1:length(vel)
        sum_std1 = 0;
        %sum_std2 = 0;
        for i = 1:avg
            for m = 1:(updates/frame)
                sum_std1 = sum_std1 + std(error1_avg(i,j,k,1+(m-1)*(frame*fs):m*frame*fs));
        %        sum_std2 = sum_std2 + std(error2_avg(i,j,k,1+(m-1)*(frame*fs):m*frame*fs));
            end
        end
        avg_std_dev1(j,k) = sum_std1/(avg*(updates/frame));
        %avg_std_dev2(j,k) = sum_std2/(avg*(updates/frame));
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Average linear pc data
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:length(rate)
%for k = 1:length(vel)
    sum_std3 = 0;
    for i=1:avg
        for m = 1:(updates/frame)
            sum_std3 = sum_std3 + std(error3_avg(i,k,1+(m-1)*(frame*fs):m*frame*fs));
        end
    end
    avg_std_dev3(k) = sum_std3/(avg*(updates/frame));
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Generate Plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


time = 0:time_step:(length(rw)-1)*time_step;            %  Total time [s]

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Plot of STD vs update rate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

up_per_sec = (fs./rate).*1500;
plot(up_per_sec,avg_std_dev1(1,:),'rv-',up_per_sec,avg_std_dev3(:),'kh-')
title({['Standard Deviation of Power Control Error as a Function of Update Rate'],['Velocity = 10 km/h']})
xlabel('Update Rate [updates/sec]')
ylabel('Standard Deviation of Power Control Error')
grid on
legend('UMTS/IMT2000','Our Algorithm')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Plot of STD vs Vel for 1 increment
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (little_plot ==1)

plot(vel,avg_std_dev2(1,:),'gx-',vel,avg_std_dev1(1,:),'rv-',vel,avg_std_dev3(:),'kh-')
xlabel('Velcocity [km/h]')
ylabel('Standard Deviation of Power Control Error [dB]')
title('Performance of Power Control Algorithms')
legend('IS-95 = 1 dB inc','UMTS/IMT2000 = 1 dB inc','Our Algorithm')
grid on
end
%print -djpeg90 umts_is95_us.jpg
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Plot of power control and error vs time
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (time_plot == 1)

subplot(8,1,4:5)

plot(time(1:updates*fs), pwr(1:updates*fs), time(1:updates*fs),pc1(1:updates*fs),'r', ...
    time(1:updates*fs),error1(1:updates*fs),'g')
xlabel('Time [s]')
ylabel('Normalized Power [dB]')
title('UMTS/IMT2000 Power Control Algorithm')
legend('Normalized Received Power','Normalized Power Control','Power Control Error')
grid on

subplot(8,1,1:2)

plot(time(1:updates*fs), pwr(1:updates*fs), time(1:updates*fs),pc2(1:updates*fs),'r', ...
    time(1:updates*fs),error2(1:updates*fs),'g')
xlabel('Time [s]')
ylabel('Normalized Power [dB]')
title('IS-95 Power Control Algorithm')
legend('Normalized Received Power','Normalized Power Control','Power Control Error')
grid on

subplot(8,1,7:8)

plot(time(1:updates*fs), pwr(1:updates*fs), time(1:updates*fs),pc3(1:updates*fs),'r', ...
    time(1:updates*fs),error3(1:updates*fs),'g')
xlabel('Time [s]')
ylabel('Normalized Power [dB]')
title('Our Power Control Algorithm')
legend('Normalized Received Power','Normalized Power Control','Power Control Error')
grid on
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  Plot of STD vs Vel for various increments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (big_plot == 1)
    

plot(vel,avg_std_dev1(1,:),'mx-',vel,avg_std_dev1(2,:),'r*-',vel,avg_std_dev1(3,:),...
    'gs-',vel,avg_std_dev1(4,:),'bd-',vel,avg_std_dev1(5,:),'cv-',vel,avg_std_dev3(:),'kh-')
legend('Step = +/- 0.25 dB','Step = +/- 0.5 dB','Step = +/- 1 dB','Step = +/- 2 dB','Step = +/- 3 dB','Our Algorithm')
grid on
xlabel('Velocity of Mobiles [km/h]')
ylabel('Standard Deviation of Power Control Error [dB]')
title({['Standard Deviation of Power Control Error'],['Using UMTS/IMT2000 Power Control Algorithm']})
end
%print -djpeg90 umts_vs_inc