n = 100;
data=rand(1,n)<=0.5;
sigA = 0.1;
sigPhi = 0.1;

figure (1)
axis([-2 2 -2 2]) % Wyznaczenie limitu osi
hold on;
title('Wykres 16QAM dla wspolrzednych biegunowych, amplitudy i fazy')
xlabel('A*cos[phi]')
ylabel('Y*sin[phi]')

n=length(data);
A=ones(1,n)+normrnd(0,sigA,1,n) % Tworzenie macierzy z szumem
phi=zeros(1,n);
f=70; % czestotliwosc
t=0:1/f:n-1/f;
y=zeros(1,length(t));%
hold on;
for i=1:n  % Poczatek petli glownej w ktorej ustawiamy faze+szum i dane sinusa

  if(data(i)==0 && data(i+1)==0 && data(i+2)==0 && data(i+3)==0)
      phi(i)=pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),A(i)*sin(phi(i)),'bo') %blue circles

  elseif(data(i)==0 && data(i+1)==0 && data(i+2)==0 && data(i+3)==1)
      phi(i)=3*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),A(i)*sin(phi(i)),'ro') %red circles

  elseif(data(i)==0 && data(i+1)==0 && data(i+2)==1 && data(i+3)==0)
        phi(i)=5*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),A(i)*sin(phi(i)),'go') %green circles

  elseif(data(i)==0 && data(i+1)==0 && data(i+2)==1 && data(i+3)==1)
       phi(i)=7*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),A(i)*sin(phi(i)),'yo') %yellow circles

  elseif(data(i)==0 && data(i+1)==1 && data(i+2)==0 && data(i+3)==0)
      phi(i)=pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),A(i)*sin(phi(i)),'y+') %blue circles

  elseif(data(i)==0 && data(i+1)==1 && data(i+2)==0 && data(i+3)==1)
      phi(i)=3*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),A(i)*sin(phi(i)),'g+') %red circles

  elseif(data(i)==0 && data(i+1)==1 && data(i+2)==1 && data(i+3)==0)
        phi(i)=5*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),A(i)*sin(phi(i)),'b+') %green circles

  elseif(data(i)==0 && data(i+1)==1 && data(i+2)==1 && data(i+3)==1)
       phi(i)=7*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),A(i)*sin(phi(i)),'r+') %yellow circles

  elseif(data(i)==1 && data(i+1)==0 && data(i+2)==0 && data(i+3)==0)
      phi(i)=pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'rx') %blue circles

  elseif(data(i)==1 && data(i+1)==0 && data(i+2)==0 && data(i+3)==1)
      phi(i)=3*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'bx') %red circles

  elseif(data(i)==1 && data(i+1)==0 && data(i+2)==1 && data(i+3)==0)
        phi(i)=5*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'yx') %green circles

  elseif(data(i)==1 && data(i+1)==0 && data(i+2)==1 && data(i+3)==1)
       phi(i)=7*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'gx') %yellow circles

  elseif(data(i)==1 && data(i+1)==1 && data(i+2)==0 && data(i+3)==0)
      phi(i)=pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'gs') %blue circles

  elseif(data(i)==1 && data(i+1)==1 && data(i+2)==0 && data(i+3)==1)
      phi(i)=3*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'rs') %red circles

  elseif(data(i)==1 && data(i+1)==1 && data(i+2)==1 && data(i+3)==0)
        phi(i)=5*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'ys') %green circles

  elseif(data(i)==1 && data(i+1)==1 && data(i+2)==1 && data(i+3)==1)
       phi(i)=7*pi/4+normrnd(0,sigPhi);
      y(1,70*(i-1)+1:70*i)=A(i)*cos(2*pi*t(1,70*(i-1)+1:70*i));
      plot(2*A(i)*cos(phi(i)),2*A(i)*sin(phi(i)),'bs') %yellow circles
  end

end

%amp=[0:0.05:5]; %badamy amplitude i faze dla wszystkich mozliwosci od 0 d0 5 co 0.05
%faza=[0:0.05:5];
%bity=1000;
%[AMP,FAZA]=meshgrid(amp, faza);
%BER=zeros(length(amp), length(faza);
%for i=1:length(amp)
  %for j=1:length(faza)
    %BER(i,j)=PSK(bity,amp(i), faza(j));
    %kazdy element BER wytwarzamy sprawdzajac cala
    %funkcje PSK dla podanych amplitudy i fazy
  %end
%endmesh(AMP,FAZA,BER)