clear all; close all;
rhoM = [0.5 0.99];
figure; hold on;
for k = 1:length(rhoM);
    N = 8; n = 10000;
    %---------------------------
    % Message Construction
    m_Re = 2*[rand(N,n)>0.5]-1; m_Im = 2*[rand(N,n)>0.5]-1;
    m = m_Re + m_Im*1j;
    %---------------------------
    %Channel Construction
    omega = 2*pi/N*(0:N-1);
    rho = rhoM(k);
    z1 = rho*exp(1j*2*pi/3); z2 = 1/rho; N = 8;
    Eh = 1 + (abs(z1+z2))^2 + (abs(z1)*abs(z2))^2;
    z = exp(j*omega);
    H = (1-z1*z.^(-1)).*(1-z2*z.^(-1));
    H = H/sqrt(Eh);
    h = [1 -(z1+z2) z1*z2]; h = h/sqrt(Eh);
    v = length(h)-1; %Memory length
    %---------------------------
    % OFDM
    x = m;
    xt = (reshape(x,N,n)).';
    xt = (sqrt(N))*ifft(xt,[],2);
    xt = [xt(:,(N-(v-1)):N) xt]; % Add Cyclic Prefix
    xt = (reshape(xt.',((N+v)*n),1)).';
    %---------------------------
    % Convolve with channel 
    y = conv(xt,h);
    % Noise Construction
    % Q3
    % noise = 0;
    % Q4
    N0 = 1; noise = sqrt(N0/2)*(randn(size(y))+1j*randn(size(y)));
    % Add noise
    y = y+noise;
    %---------------------------
    % At receiver
    y = y(1:((N+v)*n));
    yt = reshape(y,(N+v),n).';
    yt = yt(:,v+1:v+N); %Eliminate cyclic prefix
    yt = (1/sqrt(N))*fft(yt,[],2); 
    r = yt;
    for i=1:length(H)
        r(:,i) = yt(:,i)/H(i);
    end
    % ML Decoder
    m_hat_Re = [2*[real(r)>0]-1]; m_hat_Im = [2*[imag(r)>0]-1];
    % SER Calculation
    correct = sum([m_hat_Re == m_Re.'].*[m_hat_Im == m_Im.'])/n;
    SER = 1-correct
    plot(abs(H),SER,'k*')
    % Theoretical values
    p = Q(abs(H)*sqrt(2/N0));
    SER_Theoretical = (2*p) - (p.^2)
end
grid on
xlabel('|H_k|'); ylabel('SER');
