% Лабораторная работа 2 ТСиУПД
% Исследование характеристик АМ сигналов
clc;
close all;
clear all;
% Блок ввода данных
disp('Введите номер варианта [1 . . . 30]');
n=input('n=');
n=round(n);
if n<=0;
  n=1;
elseif n>30;
  n=30;
end
disp('Исходные данные варианта');disp(n);
disp('Введите Количество отсчетов N [1000...6000]');
N=input('N=');
N=round(N);
if N<=1000;
  N=1000;
elseif N>6000;
  N=6000;
end
disp('Введите Частоту дискретизации dF [400...5000]');
dF=input('dF=');
if dF<=400;
  dF=400;
elseif dF>5000;
  dF=5000;
end
disp('Введите МИН значение амплитуды сигнала AMin [-250...250]');
AMin=input('AMin=');
if AMin<=-250;
  AMin=-250;
elseif AMin>250;
  AMin=250;
end
disp('Введите МАКС значение амплитуды сигнала AMax [-250...250]');
AMax=input('AMax=');
if AMax<=-250;
  AMax=-250;
elseif AMax>250;
  AMax=250;
end
if AMax<=AMin;
  AMax=AMin+20;
elseif AMax>270;
  AMax=270;
end
disp('Введите амплитудный интервал дискретизации сигнала dA [0.01...0.5]');
dA=input('dA=');
if dA<0.01;
  dA=0.01;
elseif dA>0.5;
  dA=0.5;
end
disp('Введите форму исх.сигнала A [1-синус,2-меандр(def),3-треуг.сиг-л,4-калибр.]');
Aform=input('Aform=');
Aform=floor(Aform);
if ((Aform>=1) & (Aform<=4))
  Aform;
else
  Aform=2;
endif
disp('Введите индекс модуляции MA [0.05...2] 1-(default)');
MA=input('MA=');
if ((MA>=0.05) & (MA<=2));
  MA;
else
  MA=1;
end
disp('Введите отношение сигнал/шум в разах SN [0.05...2000] 2000-(default)');
SN=input('SN=');
if ((SN>=0.05) & (SN<=2000))
  SN;
else
  SN=2000;
end
disp('Введите форму перекр.помехи APF(0-без помехи-def)[1-синус,2-меандр,3-треуг.сиг-л,4-калибр.]');
APF=input('APF=');
APF=round(APF);
if ((APF>=1) & (APF<=4))
  APF;
else
  APF=0;
end
disp('Введите уровень перекр.помехи UPF(0-def) [0.01...1]');
UPF=input('UPF=');
if ((UPF>=0.01) & (UPF<=1))
  UPF;
else
  UPF=0;
end
disp('Количество отсчетов N = ');disp(N);
disp('Частота дискретизации dF = ');disp(dF);
disp('МИН знач.пар.ист.сообщ. AMin = ');disp(AMin);
disp('МАКС знач.пар.ист.сообщ. AMax = ');disp(AMax);
disp('Интервал дискр.ампл.сигнала dA = ');disp(dA);
disp('Форма исх. сигнала Aform = ');disp(Aform);
disp('Индекс модуляции MA = ');disp(MA);
disp('Отношение сигнал/шум в разах SN = ');disp(SN);
disp('Форму перекр.помехи APF = ');disp(APF);
disp('Уровень перекр.помехи UPF = ');disp(UPF);
%
t=0:1:N; % Формируем сигнал
Fmax=dF/2.4; 
if ((Aform>=1) & (Aform<=3)) 
  k=[1,2,3,4,6,8,10];
else
k=1;
end
difARan=(abs(AMax-AMin))/N; dT=1/dF;
Adelta=abs((AMax-AMin)/2); Asdv=AMax-Adelta;
ARan1=ones(1,N);
Es=ones(1,N);
nFi=0; % нач. фаза несущей частоты
nPHi=0; % нач.фаза для синуса;
% массив для ошибок
FrSig=ones(1,length(k));
ErrSig=ones(1,length(k));
% Ехаем ту-ту
for j=1:length(k)
A11=Adelta*sin(k(j)*Fmax*(dT*t)/10+nPHi)+Asdv; % Sinus
A12=Adelta*square((k(j)*Fmax*(dT*t)/10),50)+Asdv; % Meandr
A13=Adelta*sawtooth((k(j)*Fmax*(dT*t)/10),0.5)+Asdv; % Triangl.Wave
A14=AMin+difARan*t; % калибровка
%
switch Aform
case 1
disp('Sinus');
A=A11;
Fsignala=k(j)*Fmax/10
case 2
disp('Meandr');
A=A12;
Fsignala=k(j)*Fmax/10
case 3 
disp('Triangl');
A=A13;
Fsignala=k(j)*Fmax/10
case 4
disp('Kalibrovka');
A=A14;
Fsignala=1/(2*pi);
otherwise
end
%

% ОПСЬ....
%Nk=8; % Задаем разрядность квантования
% нормируем и квантуем по Nk
ZSig=A;
ZSig_max=max(ZSig);
ZSig_min=min(ZSig);
Nk=ceil(log2((ZSig_max-ZSig_min)/dA));
% определяем поправочн. коэфф.
Kpr=(abs(ZSig_max-ZSig_min))/(2^(Nk)-2);
% NorKvant
ZSig_NK=fix(((ZSig.-ZSig_min)/(ZSig_max-ZSig_min))*(2^Nk-2))+1; % 1-но полярн. полож.
ZSig_NK_Polar=ZSig_NK-2^(Nk-1); % 2-x полярн. симм.
ZSig_NKmax=max(ZSig_NK); ZSig_NKmin=min(ZSig_NK);
ZSig_NKmaxP=max(ZSig_NK_Polar); ZSig_NKminP=min(ZSig_NK_Polar);
%
MZSig_1=mean(ZSig); % опред. пост. составл.
MZSig_1NK=mean(ZSig_NK); % опред. пост. составл.
MZSig_1NK2=MZSig_1NK-2^(Nk-1);
MZSig_1NKP=mean(ZSig_NK_Polar); % опред. пост. составл.
% запись в файл данные с квантованием по уровню двоичн. 1...(2^Nk-1);
%dlmwrite('ZSig_NK.txt',ZSig_NK,',');
% запись в файл данные с квантованием по уровню двуполярные
%dlmwrite('ZSig_NKP.txt',ZSig_NK_Polar,',');
% готовим сигнал к модуляции
fnA=dF*20; % частота модуляции (несущей)
dTn=1/fnA; % период несущей
zNN=1:1:N*(fnA/dF); % новое число отсчетов
MfdA=2^(Nk-1);
ZSig_NKins=interp1(t,ZSig_NK_Polar,0:1/(fnA/dF):N,'spline');
% Модулируем АМ
zNN=1:1:N*(fnA/dF);
Akan=ones(1,length(zNN)); %zf(1)=0;
for i=1:length(zNN)
Akan(i)=(MfdA.+MA*ZSig_NKins(i))*cos(fnA.*(dTn.*zNN(i)).+nFi);
endfor
% готовим помеху
AP11=Adelta*sin(Fmax*(dT*t)/12)+Asdv; % Sinus
AP12=Adelta*square((Fmax*(dT*t)/12),50)+Asdv; % Meandr
AP13=Adelta*sawtooth((Fmax*(dT*t)/12),0.5)+Asdv; % Triangl.Wave
AP14=AMin+difARan*t; % калибровка
%
switch APF
case 1
disp('Pomeha-Sinus');
AP=AP11;
Fpomehi=Fmax/12
case 2
disp('Pomeha-Meandr');
AP=AP12;
Fpomehi=Fmax/12
case 3 
disp('Pomeha-Triangl');
AP=AP13;
Fpomehi=Fmax/12
case 4
disp('Pomeha-Kalibrovka');
AP=AP14;
Fpomehi=1/(2*pi);
otherwise
disp('KANAL BEZ Pomehi');
AP=zeros(1,length(t));
end
%
ZSigP=AP;
ZSig_NKP=fix(((ZSigP.-ZSig_min)/(ZSig_max-ZSig_min))*(2^Nk-2))+1; % 1-но полярн.пом.
ZSig_NKP_Polar=ZSig_NKP-2^(Nk-1); % 2-x полярн.симм.помеха
% интерп. для модуляции
ZSig_NKPinp=interp1(t,ZSig_NKP_Polar,0:1/(fnA/dF):N,'spline');
% модулируем помеху
Apom=ones(1,length(zNN)); %zf(1)=0;
for i=1:length(zNN)
Apom(i)=(MfdA.+MA*ZSig_NKPinp(i))*cos(fnA.*(dTn.*zNN(i)).+nFi);
endfor
% Моделируем в канале
AkanSN=awgn(Akan,SN,(2^(Nk)),'linear').+UPF.*Apom; % добавляем шум и помеху
%
Aprin=abs(AkanSN); % детектируем АМ балансным дет-ром
%
Aidet=decimate(Aprin,(fnA/dF)); % фильтруем НЧ и прореживаем до исходного
Ur1A=mean(Aidet);
Acalc=schtrig(Aidet,Ur1A,rst=1);
Acalc_Nk=Acalc*(2^Nk-2)+1;
Adet_max=max(Aidet);
Adet_min=min(Aidet);  % Нормируем и приводим к уровню и вып. АЦП
Adet_NK=fix(((Aidet.-Adet_min)/(Adet_max-Adet_min))*(2^Nk-2))+1;
%  восстановление исходного сигнала;
AdetVOST=((Adet_NK-1).*Kpr).+ZSig_min;
AcalcVOST=((Acalc_Nk-1).*Kpr).+ZSig_min;
% Ошибка при восстановлении
Zosh=ZSig(1:N).-AdetVOST;
Zoshcalc=ZSig(1:N).-AcalcVOST;
disp('Енергия ошибки восстановленного сигнала');
Zoshsr=meansq(Zosh)
Zoshproc=(abs(Zosh)./abs(ZSig(1:N)));
% ту
% Рисуем сигнальчики в канале
Amean=mean(Akan);
Astd=std(Akan);
%
hf11=figure(); % Рисуем на графике результат (6000 точек)
hold on; title('First 6000 samples of KANAL OUT PD A','fontsize',14);
plot(zNN(1:6000),Akan(1:6000),'b','linewidth',2); xlabel('samples A, n=1...6000'); ylabel('KANAL OUT PD (A)');
plot(zNN(1:6000),ones(1,6000).*Amean,'r','linewidth',3);
plot(zNN(1:6000),ones(1,6000).*Amean.+Astd,'g','linewidth',2);
plot(zNN(1:6000),ones(1,6000).*Amean.-Astd,'g','linewidth',2);
grid on;
hf21=figure(); % Рисуем на графике фрагмент (6000 точек)
hold on; title('First 6000 samples of KANAL IN PR A+N+P','fontsize',14);
plot(zNN(1:6000),AkanSN(1:6000),'b','linewidth',2); xlabel('samples A, n=1...100'); ylabel('KANAL IN PR (A+N+P)');
plot(zNN(1:6000),ones(1,6000).*Amean,'r','linewidth',3);
plot(zNN(1:6000),ones(1,6000).*Amean.+Astd,'g','linewidth',3);
plot(zNN(1:6000),ones(1,6000).*Amean.-Astd,'g','linewidth',3);
grid on;
%

% Картинка суммарная вход-выход инфо
hf61=figure();
subplot(311);
plot(t(1:N),ZSig(1:N),"b;signal;",
     t(1:N),((Adet_NK(1:N)-1).*Kpr).+ZSig_min,";signal-outPR;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
subplot (312);
plot(t(1:N),sort((ZSig(1:N))),'b;sign-sort;',
     t(1:N),sort(((Adet_NK(1:N)-1).*Kpr).+ZSig_min), ';signPR-sort;');
xlabel ("Rang");
ylabel ("Amplitude");
subplot (313);
plot(t(1:N),sort(ZSig(1:N))-sort(((Adet_NK(1:N)-1).*Kpr).+ZSig_min),'b;errPR-sort;');
xlabel ("Rang");
ylabel ("Error");
% Картинка Ошибка при восстановлении
hf71=figure();
plot(t(1:N),ZSig(1:N),"b;signal;",
     t(1:N),Zosh(1:N),";signal oshibki;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
%
hf81=figure();    
plot(t(1:N),Zoshproc(1:N),"b;otnosit oshibka kvant;",
     t(1:N),ones(1,N)*Zoshsr,";energiya oshibki kvant;");
%
hf91=figure();
plot(t(1:N),A(1:N),"b;signal;",
     t(1:N),AdetVOST(1:N),";signal vost;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
% ту кубышка
FrSig(j)=FrSig(j)*Fsignala;
ErrSig(j)=ErrSig(j)*Zoshsr;
% Блок сохранения результатов
NameFile=strcat('Lab_TSiUPD_2_AM_','Fsig_',num2str(round(Fsignala),5),'Hz','.txt');
%
disp('Желаете сохранить рез-ты в файлы Lab_TSiUPD_2_AM.txt и графики в ***.pdf?');
disp('[ 1 - ДА, желаю; 0 - НЕТ, не желаю (по умолчанию) ]'); 
Rez=input('Rez=');
if Rez==1
% Запись результатов в файл 
f_stud_rezy=fopen(NameFile,'wt');
fprintf(f_stud_rezy,"Исходные данные варианта \t%u\n\n",n);
fprintf(f_stud_rezy,"Количество отсчетов N = \t%u\n\n",N);
fprintf(f_stud_rezy,"Частота дискретизации [Hz] dF = \t%u\n\n",dF);
fprintf(f_stud_rezy,"МИН знач.ампл.сигнала AMin = \t%u\n\n",AMin);
fprintf(f_stud_rezy,"МАКС знач.ампл.сигнала AMax = \t%u\n\n",AMax);
fprintf(f_stud_rezy,"Интервал дискр.ампл.сигнала dA = \t%u\n\n",dA);
fprintf(f_stud_rezy,"Форма исх. сигнала Aform = \t%u\n\n",Aform);
fprintf(f_stud_rezy,"Индекс модуляции MA = \t%u\n\n",MA);
fprintf(f_stud_rezy,"Отношение сигнал/шум в разах SN = \t%u\n\n",SN);
fprintf(f_stud_rezy,"Форма перекр.помехи APF = \t%u\n\n",APF);
fprintf(f_stud_rezy,"Уровень перекр.помехи UPF = \t%u\n\n",UPF);
fprintf(f_stud_rezy,"  \t\n\n");
fprintf(f_stud_rezy,"Результаты вычислений параметров сигнала. А \t\n\n");
%fprintf(f_stud_rezy,"Матожидание (сред.арифм.) Amean = \t%u\n\n",Amean);
%fprintf(f_stud_rezy,"СКО Astd = \t%u\n\n",Astd);
fprintf(f_stud_rezy,"Кол-во уровней квант. дискрет.сигнала A = \t%u\n\n",2^Nk);
fprintf(f_stud_rezy,"Частота сигнала А [Hz] Fsignala = \t%u\n\n",Fsignala);
fprintf(f_stud_rezy,"Енергия ошибки восст.сигнала А Zoshsqr = \t%u\n\n",Zoshsr);
fclose(f_stud_rezy);
% сохраняем графики
print(hf11,strcat('Laba_TSiUPD_2_AM_Graf_6000N_OUT_Kanal_','Fsig_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,300");
print(hf21,strcat('Laba_TSiUPD_2_AM_Graf_6000N_IN_Kanal_','Fsig_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,300");
print(hf61,strcat('Lab_TSiUPD_2_AM_','Fsig_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,400");
print(hf71,strcat('Lab_TSiUPD_2_AM_','Graf_Sign_AbsolErr_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,400");
print(hf81,strcat('Lab_TSiUPD_2_AM_','Graf_OtnErr_EnjErr_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,600");
print(hf91,strcat('Lab_TSiUPD_2_AM_','Graf_SignA_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,300");
else
Rez=0
endif
close all;
continue
endfor
hf103=figure();    
plot(FrSig,ErrSig,"*b;Energiya oshibki;");
xlabel('frequency sign.A, [Hz]'); ylabel('Norm. energiya err.A');
print(hf103, "Laba_TSiUPD_2_AM_Graf_Energii_Err_Lin.pdf", "-dpdf","-S600,300");
% конец программы


