%The received signal 
clear
%Duration 
t = 0 : 1/1e3 : 0.999; 
 
%Source 
a = randint(1, 20, 2); 
a=[1,1,1,0,0,0,0,1,1,0,0,1,1,1,0,1,0,1,1,1];
%Serial to Parallel 
Ai = 2*a(1:2:20) - 1; 
Aq = 2*a(2:2:20) - 1; 
 
%QAM Modulation 
g = ones(1, 100); 
g = [g, g, g, g, g, g, g, g, g, g]; 
s = Ai(ceil(10*t+0.01)).*g.*cos(2*pi*100*t)+ Aq(ceil(10*t+0.01)).*g.*sin(2*pi*100*t); 
 
 
Y=fftseq(s,512);
Pyy = Y.* conj(Y) / 512;
f = 1000*(0:256)/512;

%Drawing 



r = awgn(s, 0); 
 
%Two correlators 
T = 0:1/1e3:0.099; 
Icorrelator = cos(2*pi*100*T); 
Qcorrelator = sin(2*pi*100*T); 
 
%QAM demodulation 
%The result is two dimensions variant 
for n = 0:9 
    result(n+1, :) = [trapz(T, r(n*100+1:n*100+100).*Icorrelator),trapz(T, r(n*100+1:n*100+100).*Qcorrelator)]; 
end 
 
%Calculate the signal points in vector space 
%The s0, s1, s2, s3 are two dimensions variant 
s0 = [trapz(T, -1*ones(1, 100).*cos(2*pi*100*T).*Icorrelator),trapz(T, -1*ones(1, 100).*sin(2*pi*100*T).*Qcorrelator)]; 
s1 = [trapz(T, -1*ones(1, 100).*cos(2*pi*100*T).*Icorrelator),trapz(T, 1*ones(1, 100).*sin(2*pi*100*T).*Qcorrelator)]; 
s2 = [trapz(T, 1*ones(1, 100).*cos(2*pi*100*T).*Icorrelator),trapz(T, -1*ones(1, 100).*sin(2*pi*100*T).*Qcorrelator)]; 
s3 = [trapz(T, 1*ones(1, 100).*cos(2*pi*100*T).*Icorrelator),trapz(T, 1*ones(1, 100).*sin(2*pi*100*T).*Qcorrelator)]; 
 
%Detection 
for n = 0:9 
    [C, I(n+1)] = min([norm(result(n+1,:)-s0), norm(result(n+1,:)-s1), norm(result(n+1,:)-s2),norm(result(n+1,:)-s3)]); 
end 
 
%The input bit stream 
input = a; 
 
%The ouput stream 
output(1:2:20) = floor((I-1)/2); 
output(2:2:20) = mod((I-1), 2); 
 
%Drawing 
figure(1)
subplot(4, 1, 1); 
plot(t, input(ceil((100*t+0.1)/5))); 
axis([0 1 -0.2 1.2]); 
subplot(4, 1, 2);
plot(t, s); 
subplot(4,1,3);
plot(f,Pyy(1:257))
subplot(4, 1, 4); 
plot(t, output(ceil((100*t+0.1)/5))); 
axis([0 1 -0.2 1.2]); 
figure(1)
scatterplot(Ai+j*Aq);