% HW#4v2RX: 64-QAM Rx with MUSIC DOA Estimation and MVDR Beamforming 
% DCChang  NCU/CE 2022/5/12
%
clear all;
load r_hw4      % MxN received signal matrix

Nsig=4;
DOA0range=[-26 -5];
d=1/2;
M=16;
p=[0:M-1]';
N=length(r);

%%%%%%%%%%%%%%%%%%%%%%% RX %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% DOA and A0 Estimation
% Compute covariance matrix Ruu
Ruu=Rcompute(r(:,1:N));    %<================================= HW#4
%
[doa,spec_music,specang] = musicdoa(Ruu,Nsig);
% Determine DOA0
doa0=determineDOA(doa,DOA0range);   %<======================== HW#4
%
A0e=exp(1j*p*2*pi*d*sin(doa0*pi/180));   
% MVDR Algorithm
W=mvdr(A0e,Ruu);                   %<========================= HW#4
%
for k=1:N
   qamd_mvdr(k)=W'*r(:,k); 
end    
phis = 1j*2*pi*d*sin(-pi/2:pi/400:pi/2);
for k=1:length(phis)
    A = exp(phis(k).*[0:M-1]');
    beam(k)= W'*A;   
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
plot(specang,10*log10(spec_music))
xlabel('Angle (deg)')
ylabel('Magnitude (dB)')
title('DOA Spectrum')
grid

figure(2)
plot(qamd_mvdr,'.b')
hold on
qam=qammod([0:63],64)/sqrt(42); 
plot(qam,'or')
hold off
legend('MVDR Beamforming')
xlabel('I-channel')
ylabel('Q-channel')
axis(1.5*[-1 1 -1 1])
title('Received 64-QAM constellation'); 

figure(3)
x=linspace(-pi/2,pi/2,length(phis))*180/pi;
plot(x,10*log10(abs(beam)),'b')
hold on
plot(doa0+[0 10^(-6)],[-80 10],'--r')
hold off
grid on
legend('MVDR Beamforming','DOA0')
xlabel('Angle (deg)'),
ylabel('Magnitude (dB)')
axis([-90 90 -30 5])
title('MVDR Beampattern')

figure(4)
polardb([-pi/2:pi/400:pi/2],10*log10(abs(beam).^2),-20,'b'); %<= HW#4
title('MVDR Polar Beampattern')