%% Project 2
clc, clear all, close all

% == Parameters == %
delay_spread = 5.4e-6;
fc = 2e9;
v = 1e-100; % speed of receiver
fD = (v/(physconst('LightSpeed')))*fc;
B = 1e6;
Noise = 2.07e-20*B; % in [W]
fs = 1e6;
Ts = 1/fs; % Sample time
fDTs = fD*Ts; % Normalized Doppler frequency
E = 1;
N = 100;
L = ceil(delay_spread/Ts);
Ncp = L-1;
tau = (0:1:L-1); % Path delays in samples
P = (1/length(tau))*ones(1,length(tau)); % Flat power delay profile
PL = 1; %10^(-10.1);
M = 30; % Number of time samples
Q = 50; % Number of Pilotsymbols

if (N + Ncp)*fDTs*10 > 1 % Check constraint
    error('N or Ncp to large')
end
if N*Ts > 1/(2*fD) % Check if N*Ts is smaller than coherence time
    error('N too large')
end

% == Transmitter  == %
symbs = zeros(N,M);
y = zeros(N,M);

%[symbs(1:Q,1), ~] = qpsk_sequence(Q,E);
symbs(1:Q,1) = (sqrt(E))*(1+1i)*ones(Q,1);

symbs = repmat(symbs(:,1),1,M);


for m = 1:M
    [symbs(Q+1:end,m), ~] = qpsk_sequence(N-Q,E);
    symbs_ifft = sqrt(N)*ifft(symbs(:,m),N);
    z = [symbs_ifft(end-Ncp+1:end);symbs_ifft];
    
    % == Channel  == %
    [z, ~] = Fading_Channel(z, tau, fDTs, P);
    z = z(1:end-L+1); % Remove delayed symbols
    
    %c_output = awgn(z,i,'Measured');
    z = z*sqrt(PL); % Pathloss
    c_output = z + (randn(length(z),1) +1i*rand(length(z),1))*sqrt(Noise/N);   % add noise
    
    % == Reciever  == %
    r2 = c_output(Ncp+1:end); %remove cp
%     h = (r2(1:Q)./symbs(1:Q,m));
%     ch = fft(h,N,2); % FFT of channel
%     ch = mean(ch,1);

    h = mean((r2(1:Q)./symbs(1:Q,m)),1);
    ch = fft(h,N); % FFT of channel
      
    y(:,m) = fft(r2,N)/sqrt(N); % FFT of received symbols
    y(:,m) =y(:,m)./ch.';  % zero forcing equalizer with the first tap of the channel frequency response.
    
end

scatterplot(y(1,:)), hold on; % scatterplot for the first subcarrier
title('Scatterplot n=1');
scatterplot(y(ceil(0.25*N),:)); 
scatterplot(y(ceil(0.5*N),:));
scatterplot(y(N,:));