pkg load communications;
clc;
close all;
clear all;
nSym = 1000; %число точек данных

M=16;%задать число точек PAM
Nfaza=pi/21; %(0.99*pi);%задать угол вращения созвездия точек PAM
k=log2(M);

%-----------------Signal-to-noise ratio----------
pmin=-6; pmax=16; dp=2; %пределы и шаг шума
p=pmin:dp:pmax;
Nzcap=10; %число точек усреднения
%-----------------Transmitter--------------------
Spwr=3; %мощность сигнала дБ
d=ceil(M.*rand(1,nSym))-1;

s=pammod(d,M,Nfaza,"gray");%модуляция PAM кодом Грея,"bin" для двоичн.
%-------------- Channel ----------------
errcap=zeros(1,length(p));
for ip1=1:length(p);
  er1=0;
  for n=1:Nzcap;
z=awgn(s,p(ip1),Spwr);% помеха
ZkanWithawgn=real(z)*cos(2*pi)+imag(z)*sin(2*pi);
ZkanWithoutAwgn=real(s)*cos(2*pi)+imag(s)*sin(2*pi);


dCap=pamdemod(s,M,Nfaza,"gray");
dCap2=pamdemod(z,M,Nfaza,"gray");

%figure;plot(dCap);
%figure; plot(dCap2);

sp=pammod(dCap2,M,Nfaza,"gray");
KanalVost = real(sp)*cos(2*pi)+imag(sp)*sin(2*pi);
%figure; plot(KanalVost,'linewidth',1); xlabel('time, s'); ylabel('Amplitude A');


BitError = 0;
 for i =1:length(d);
  if(d(i)!=dCap2(i));
  BitError=BitError+1;
  endif;
 endfor;
 er1=er1+BitError;
 endfor;
 er2=(er1/Nzcap)/nSym;
 errcap(ip1)=er2;
end;
%
%созвездия
ht1=figure(); plot(s,'rx','linewidth',2); xlabel('I'); ylabel('Q');
ht2=figure(); plot(z,'gx','linewidth',2); xlabel('I'); ylabel('Q');
%канальные сигналы
ht3=figure(); plot(ZkanWithoutAwgn,'b','linewidth',1); xlabel('time, s'); ylabel('Amplitude A');
%ht13=figure(); periodogram(ZkanWithoutAwgn) %периодограмма (мгновенный спектр)
%ht23=figure(); pwelch(ZkanWithoutAwgn,256,[],[],'no-strip') %спектр мощности
ht4=figure(); plot(ZkanWithawgn,'c','linewidth',1); xlabel('time, s'); ylabel('Amplitude A');
%ht14=figure(); periodogram(ZkanWithawgn) %периодограмма (мгновенный спектр)
%ht24=figure(); pwelch(ZkanWithawgn,256,[],[],'no-strip') %спектр мощности
ht5=figure(); plot(KanalVost,'m','linewidth',1); xlabel('time, s'); ylabel('Amplitude A');
%ht15=figure(); periodogram(KanalVost) %периодограмма (мгновенный спектр)
%ht25=figure(); pwelch(KanalVost,256,[],[],'no-strip') %спектр мощности
%график ошибок
ht6=figure(); semilogy(p,errcap,'k','linewidth',2); xlabel('Signal-to-noise ratio, dB'); ylabel('BitError');
%конец
 