%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Part a  from report6_DSP.pdf
% To build upon Part I, we are to implement a LMS equalizer using Matlab. This first part is initial
% timing acquisition. To do this, we need to create a data sequence. We use an m-sequence for
% this with state vector of 8, a connection polynomial of 37 and a length of 127. This gives us a
% PN code of length 127. The length of 127 was purely arbitrary. The function mseq was used
% from a previous lab.
clear all
close all
data_seq = mseq(8,37,127);
%For our signal, we are required to have 8 samples per symbol. To do this, we
%basically over-sampled the signal data_seq.
samples=8;
current_chip=1;
xmt_index=1;
for i=1:length(data_seq)
    current_bit=data_seq(i);
    for j=1:samples
        xmt_bpsk(xmt_index) = current_chip*current_bit;
        xmt_index = xmt_index+1;
    end
end

% We then pass the signal xmt_bpsk through a channel filter of specifications [1 0 0 0 0.7 0 0 -0.3].
% We also have to define a matched filter of length 8 because we have 8 samples per symbol.

h_channel=[1 0 0 0 0.7 0 0 -0.3];
chan_out_shift=filter(h_channel,1,xmt_bpsk);
h_mf=[1 1 1 1 1 1 1 1];
% Next we insert a timing delay by shifting the data sequence to where it loses an arbitrary number
% of samples in the beginning.
delay=0;
data_seq=data_seq(1,delay+1:length(data_seq));
xmt_bpsk=xmt_bpsk(1,delay+1:length(xmt_bpsk));

buffer=[];
acq_stat=[];
T=8;
acq_add=0;
cnt=0;
sym=10;
for j=1:sym
    index=1;
    a=1+cnt;
    b=(T)+cnt;
    for k=a:b
        buffer=chan_out_shift(1,k:(T-1)+k);
        temp=sum(buffer.*h_mf);
        acq_stat(index) = abs(temp);
        index=index+1;
        cnt=cnt+1;
    end
%     figure(j)
%     plot(acq_stat);
    acq_add=acq_add + acq_stat;
    %[acq_max,acq_index]=max(acq_add/j)
end
[acq_max,acq_index]=max(acq_add/sym)

%%%%%%%%%%%%%part a acquisition over


%%%%%%%%%%%%%%part b tracking begin
% For part b, we are required to do symbol tracking. We use the same code as for part a to do the
% acquisition and to set up our data sequence. The only difference is we have a data sequence of
% length 2005 for this part. Once we acquire the signal and find the index of where the peaks were
% found, we must shift our data sequence by the index so that we can have maximum correlation at
% the start of each 8 samples. To accomplish the tracking, we use the Least Mean Squared method
% with blind equalization. This means that our equalizer has no prior information about our
% incoming data. We set up a filter to be used for the equalizer. This filter consists of 8 ones and 7
% zeros. We divide the filter by 8 so we normalize the 8 ones.

taps=15;
h_eq=[1 1 1 1 1 1 1 1 0 0 0 0 0 0 0]./8;
zi=zeros(1,14);
mu=0.005;
%fill tap delay line
sym_index = 1;
for i = 1:15
    new_samp = chan_out_shift(i);  %
    buffer = [new_samp zi];
    [y(sym_index),zi] = my_filter(h_eq,new_samp,zi);
end
%--------------------%
% Blind Equalization %
%--------------------%
samp_cnt = 0;
for i=16:length(chan_out_shift)  %
    samp_cnt=samp_cnt+1;
    new_samp=chan_out_shift(i);   %chan_out_shift
    buffer=[new_samp zi];
    [y(sym_index),zi]=my_filter(h_eq,new_samp,zi);
    if samp_cnt == T
        d=sign(y(sym_index));
        error(sym_index)=d-y(sym_index);
        grad=error(sym_index).*buffer;
        h_eq = h_eq + mu.*grad(length(grad):-1:1);
        sym_index = sym_index + 1;
        samp_cnt = 0;
    end
end
plot(error)
grid;
title('Error Signal');



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%part c
% sym_index=0;
% n_sym=200*samples;
% i=1;
% while sym_index < length(chan_out)
%     sym_index=sym_index+1+(1/n_sym);
%     chan_out_loss(i)=chan_out(round(sym_index));
%     i=i+1;
% end


%In doing tracking, we have to keep up with an on-time, an early, and a late statistic. So each
%statistic from part b now has 3 statistics. If the sample is early, we must take the sample above it
%and if it’s late, we take the one before it. This will constantly make sure our peak is actually a
%peak.

samp_cnt = 0;
zi_on=zi;
zi_early=zi;
zi_late=zi;
i=15;
while sym_index <= 2000
        i=i+1;
        samp_cnt=samp_cnt+1;
        new_samp_on=chan_out_shift(i);
        new_samp_late=chan_out_shift(i-1);
        new_samp_early=chan_out_shift(i+1);
        [y_on(sym_index),zi_on]=my_filter(h_eq,new_samp_on,zi_on);
        [y_early(sym_index),zi_early]=my_filter(h_eq,new_samp_early,zi_early);
        [y_late(sym_index),zi_late]=my_filter(h_eq,new_samp_late,zi_late);
        if samp_cnt == T
                d_on=sign(y_on(sym_index));
                d_early=sign(y_early(sym_index));
                d_late=sign(y_late(sym_index));
                error_on(sym_index)=d_on-y_on(sym_index);
                error_early(sym_index)=d_early-y_early(sym_index);
                error_late(sym_index)=d_late-y_late(sym_index);
                [error_stat(sym_index),error_index]=min([abs(error_on(sym_index)),...
                abs(error_early(sym_index)), abs(error_late(sym_index))]);
                if error_index==1
                    new_samp=new_samp_on;
                    output_bit(sym_index) = d_on;
                    error_stat(sym_index) = error_on(sym_index);
                    T=8;
                    zi=zi_on;
                elseif error_index==2
                    new_samp=new_samp_early;
                    output_bit(sym_index) = d_early;
                    error_stat(sym_index) = error_early(sym_index);
                    T=9;
                    zi=zi_early;
                elseif error_index==3
                    new_samp=new_samp_late;
                    output_bit(sym_index) = d_late;
                    error_stat(sym_index) = error_late(sym_index);
                    T=7;
                    zi=zi_late;
                end
                grad=error_stat(sym_index).*buffer;
                h_eq = h_eq + mu.*grad(length(grad):-1:1);
                buffer=[new_samp zi];
                sym_index = sym_index + 1;
                samp_cnt = 0;
                continue;
        end
        buffer=[new_samp_on zi_on];
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%part d
if (F == 31 & sym_index > 100)
    F_sym_index=sym_index;
    if blind == 1 %finished blind, enter training
          blind=0;
    else %finished training, enter blind
          blind=1;
    end
end



%%%%%%%%%%%%%%%%%%%%%%part e
% Part e consists of creating a BPSK signal and adding noise and passing it through a channel.
% Then we will demodulate using the Costas Loop from last lab and use our equalizer to acquire
% and track the signal to give us our output bits. To do this, we created a modulation function and
% a demodulation function. The modulation function takes in an array and multiplies it with a
% carrier. In addition, we also added noise at this stage as well. This function was called via the
% main equalizer code as:
xmt_bpsk=mod(bpsk);


% The demodulation function is called in a similar fashion and performs carrier recovery using a
% Costas Loop. The function takes in the modulated signal with channel corruption and removes
% the carrier. The resulting signal is fed into the acquisition and tracking algorithms to make
% decisions on the bits as in part d.


% chan_output=demod(chan_out_loss);