% Практическая работа 4 ТСиУПД
% Оценка эффекта маскирования сигналов при дискретизации и квантовании
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('Введите закон распред.случ.сигн. Ap [1-равн.,2-норм.(def),3-экспон.]');
Azakrasp=input('Azakrasp=');
Azakrasp=floor(Azakrasp);
if ((Azakrasp>=1) & (Azakrasp<=3))
  Azakrasp;
else
  Azakrasp=2;
endif
disp('Введите форму детерм.сигнала As [1-синус,2-меандр(def),3-треуг.сиг-л,4-лин.нар.сиг-л]');
Aform=input('Aform=');
Aform=floor(Aform);
if ((Aform>=1) & (Aform<=4))
  Aform;
else
  Aform=2;
endif
%KAs=1;
%KAp=1;
disp('Введите коэф.уровня детерм.сигнала KAs [0.1 ... 1]');
KAs=input('KAs=');
if KAs<0.1;
  KAs=0.1;
elseif KAs>1;
  KAs=1;
end
disp('Введите коэф.уровня случ.сигнала KAp [0.1 ... 1]');
KAp=input('KAp=');
if KAp<0.1;
  KAp=0.1;
elseif KAp>1;
  KAp=1;
end
disp('Количество отсчетов N = ');disp(N);
disp('Частота дискретизации dF = ');disp(dF);
disp('МИН знач.пар.ист.сообщ. AMin = ');disp(AMin);
disp('МАКС знач.пар.ист.сообщ. AMax = ');disp(AMax);
disp('Интервал дискр.ампл.сигнала dA = ');disp(dA);
disp('Форма закон распред.случ.сигн. Azakrasp = ');disp(Azakrasp);
disp('Коэф.уровня случ.сигнала KAp = ');disp(KAp);
disp('Форма детерм.сигн. Aform = ');disp(Aform);
disp('Коэф.уровня детерм.сигнала KAs = ');disp(KAs);
%
t=1:1:N; % Формируем сигнал
Fmax=dF/2.4; k=1; %[1,2,3,4,6,8,10]; %,12,15,18,21,24,30,36,42,48];
difARan=(abs(AMax-AMin))/N; dT=1/dF;
Adelta=abs((AMax-AMin)/2); Asdv=AMax-Adelta;
ARan1=ones(1,N);
Es=ones(1,N);
nPHi=0; % нач. фаза информационного сигнала-только для SINUSa
%
% равномерная ПРВ
AaRan=rand(1,N); 
% нормальная ПРВ
AbRan=randn(1,N); 
% экспоненциальная ПРВ
AcRan=rande(1,N); 
%
switch Azakrasp
case 1
disp('Uniform distribution');
ARan=AaRan;
MaxARan=max(ARan); MinARan=min(ARan);
ARan1=(ARan.-MinARan)./(MaxARan-MinARan);
Ap=(ARan1.*abs(AMax-AMin)).+AMin; 
case 2
disp('Normal distribution');
ARan=AbRan;
MaxARan=max(ARan); MinARan=min(ARan);
ARan1=(ARan.-MinARan)./(MaxARan-MinARan);
Ap=(ARan1.*abs(AMax-AMin)).+AMin; 
case 3 
disp('Exponential distribution');
ARan=AcRan;
MaxARan=max(ARan); MinARan=min(ARan);
ARan1=(ARan.-MinARan)./(MaxARan-MinARan);
Ap=(ARan1.*abs(AMax-AMin)).+AMin; 
otherwise
end
%
A1=Adelta*sin(k*Fmax*(dT*t)/10+nPHi)+Asdv; % Sinus
A2=Adelta*square((k*Fmax*(dT*t)/10),50)+Asdv; % Meandr
A3=Adelta*sawtooth((k*Fmax*(dT*t)/10),0.5)+Asdv; % Triangl.Wave
A4=AMin+difARan*t;
%
switch Aform
case 1
disp('Sinus');
As=A1;
case 2
disp('Meandr');
As=A2;
case 3 
disp('Triangl');
As=A3;
case 4
disp('Linearly growing signal');
As=A4;
otherwise
end
%
A=(KAp.*Ap).+(KAs.*As);
%
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); % опред. пост. составл.
% Исходный сигнал
% Аналитический сигнал
ZSig_H=hilbert(ZSig);
ZSig_HF=fft(ZSig_H);
ZSig_HFmod=abs(ZSig_HF)./(N/2);
f_HF=[0,(1:floor((N-1)/2))].*(1/(N*(1/dF)));
hf10=figure(); % Рисуем на графике результат 
hold on; title('Amplitude analitic A (Hilbert)','fontsize',14);
plot(t,ZSig,'b','linewidth',5,t,real(ZSig_H),'g','linewidth',2, imag(ZSig_H),'m','linewidth',2);
xlabel('sample A, [1...N]'); ylabel('Amplitude A');
grid on;
% Вещественный спектр FFT после Гилберта;
hf104=figure(); % Рисуем на графике результат 
hold on; title('Spectrum analitic A (Hilbert)','fontsize',14);
plot(f_HF,ZSig_HFmod(1:length(f_HF)),'b','linewidth',2);
xlabel('frequency A, [Hz]'); ylabel('Amplitude A');
grid on;
% ДПФ
ZSig_F=fftshift(fft(ZSig));   %'Amplitude spectrum A'
ZSig_Fmod=abs(ZSig_F)./(N/2);
ZSig_Farg=arg(ZSig_F);
%f=[0:(dF/2)/(length(ZSig_F)-1):dF/2];
f=[-(ceil((N-1)/2):-1:1),0,(1:floor((N-1)/2))].*(1/(N*(1/dF)));
%
hf11=figure(); % Рисуем на графике результат 
hold on; title('Amplitude spectrum A (DPF)','fontsize',14);
plot(f,ZSig_Fmod,'b','linewidth',2); xlabel('frequency A, [Hz]'); ylabel('Amplitude A');
grid on;
hf12=figure(); % Рисуем на графике результат 
hold on; title('Phase spectrum A','fontsize',14);
plot(f,ZSig_Farg,'b','linewidth',2); xlabel('frequency A, [Hz]'); ylabel('Phase A, [rad]');
grid on;
% СПМ
[ZSig_F1,freq]=pwelch(ZSig,N,0.5,N,dF,'onesided','plot','none','sloppy');
hf111=figure(); % Рисуем на графике результат 
hold on; title('Power spectral density A','fontsize',14);
plot(freq,ZSig_F1,'b','linewidth',2); xlabel('frequency A, [Hz]'); ylabel('Power A');
grid on;
% Двуполярный дискретный квантованный сигнал
% Аналитический сигнал
ZSig_H_NP=hilbert(ZSig_NK_Polar);
ZSig_HF_NP=fft(ZSig_H_NP);
ZSig_HFmod_NP=abs(ZSig_HF_NP)./(N/2);
ZSig_HFang_NP=arg(ZSig_HF_NP);
f_HF=[0,(1:floor((N-1)/2))].*(1/(N*(1/dF)));
hf14=figure(); % Рисуем на графике результат 
hold on; title('Amplitude analitic AD (Hilbert)','fontsize',14);
plot(t,ZSig_NK_Polar,'b','linewidth',5,t,real(ZSig_H_NP),'g','linewidth',2, imag(ZSig_H_NP),'m','linewidth',2);
xlabel('sample AD, [1...N]'); ylabel('Amplitude AD');
grid on;
% Вещественный спектр FFT после Гилберта;
hf108=figure(); % Рисуем на графике результат 
hold on; title('Spectrum analitic AD (Hilbert)','fontsize',14);
plot(f_HF,ZSig_HFmod_NP(1:length(f_HF)),'b','linewidth',2);
xlabel('frequency AD, [Hz]'); ylabel('Amplitude AD');
grid on;
hf1081=figure(); % Рисуем фазу AD_H
hold on; title('Phase Spectrum analitic AD (Hilbert)','fontsize',14);
plot(f_HF,ZSig_HFang_NP(1:length(f_HF)),'b','linewidth',2);
xlabel('frequency AD, [Hz]'); ylabel('Phase AD, [rad]');
grid on;
% ДПФ
ZSig_F_NP=fftshift(fft(ZSig_NK_Polar));   %'Amplitude spectrum A'
ZSig_Fmod_NP=abs(ZSig_F_NP)./(N/2);
ZSig_Farg_NP=arg(ZSig_F_NP);
f=[-(ceil((N-1)/2):-1:1),0,(1:floor((N-1)/2))]*(1/(N*(1/dF)));
%
hf15=figure(); % Рисуем на графике результат 
hold on; title('Amplitude spectrum AD (DPF)','fontsize',14);
plot(f,ZSig_Fmod_NP,'b','linewidth',2); xlabel('frequency AD, [Hz]'); ylabel('Amplitude AD');
grid on;
hf16=figure(); % Рисуем на графике результат 
hold on; title('Phase spectrum AD','fontsize',14);
plot(f,ZSig_Farg_NP,'b','linewidth',2); xlabel('frequency AD, [Hz]'); ylabel('Phase AD, [rad]');
grid on;
% СПМ
[ZSig_F1_NP,freq_NP]=pwelch(ZSig_NK_Polar,N,0.5,N,dF,'onesided','plot','none','sloppy');
hf112=figure(); % Рисуем на графике результат 
hold on; title('Power spectral density AD','fontsize',14);
plot(freq,ZSig_F1_NP,'b','linewidth',2); xlabel('frequency AD, [Hz]'); ylabel('Power AD');
grid on;
% gradient и discrete Laplace operator
ZSig_L=del2(ZSig);
ZSig_G=gradient(ZSig);
hf114=figure(); % Рисуем на графике результат 
hold on; title('Amplitude, gradient &  discrete Laplace operator A','fontsize',14);
plot(t,ZSig./10,'b','linewidth',4,t,ZSig_G,'g','linewidth',3,ZSig_L,'m','linewidth',2);
xlabel('sample A, [1...N]'); ylabel('Value A');
grid on;
%
ZSig_L_NK=del2(ZSig_NK_Polar);
ZSig_G_NK=gradient(ZSig_NK_Polar);
hf116=figure(); % Рисуем на графике результат 
hold on; title('Amplitude, gradient & discrete Laplace operator AD','fontsize',14);
plot(t,ZSig_NK_Polar./10,'b','linewidth',4,t,ZSig_G_NK,'g','linewidth',3,ZSig_L_NK,'m','linewidth',2);
xlabel('sample AD, [1...N]'); ylabel('Value AD');
grid on;
%
AsS=abs(As);
ApP=abs(Ap);
disp('Мощность исходного детерминированного сигнала');
PA_s=(sqrt(sumsq(AsS(1:N))).*KAs)/N
disp('Мощность исходного случайного шумоподобного сигнала');
PA_p=(sqrt(sumsq(ApP(1:N))).*KAp)/N
% SN_ish=PA_s/PA_p; SN_dB_ish=10*log10(PA_s/PA_p);
disp('Отношение сигнал/помеха исходного сигнала реализации');disp('в разах');
SN_ish=PA_p/PA_s
disp('в децибелах');
SN_dB_ish=10*log10(PA_p/PA_s)
disp('Мощность случайного шумоподобного AD сигнала в канале');
PAD_p=sum(ZSig_HFmod_NP((N/2-102):(N/2-2)))/2
disp('Мощность детерминированного AD сигнала в канале'); 
PAD_s=max(ZSig_HFmod_NP(2:dF/10))
disp('Отношение сигнал/помеха для AD в канале');disp('в разах');
SND_kan=PAD_p/PAD_s
disp('в децибелах');
SND_dB_kan=10*log10(PAD_p/PAD_s)
% Блок сохранения результатов
NameFile=strcat('PR_4SP_TSiUPD_','Var_',num2str(n,2),'_ZakRasp_',num2str(Azakrasp,1),'_Aform_',num2str(Aform,1),'.txt');
%
disp('Желаете сохранить рез-ты в файлы PR_4SP_TSiUPD_Var_ZakRasp_Aform.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,"Закон распр.случ.сигнала Azakrasp = \t%u\n\n",Azakrasp);
fprintf(f_stud_rezy,"Форма детерм.сигнала Aform = \t%u\n\n",Aform);
fprintf(f_stud_rezy,"Коэф.уровня случ.сигнала KAp = \t%u\n\n",KAp);
fprintf(f_stud_rezy,"Коэф.уровня детерм.сигнала KAs = \t%u\n\n",KAs);
fprintf(f_stud_rezy,"  \t\n\n");
fprintf(f_stud_rezy,"Результаты вычислений параметров сигнала. А \t\n\n");
fprintf(f_stud_rezy,"Кол-во уровней квант. дискрет.сигнала AD = \t%u\n\n",2^Nk);
fprintf(f_stud_rezy,"Мощность исходного детерминированного сигнала PA_s = \t%u\n\n",PA_s);
fprintf(f_stud_rezy,"Мощность исходного случайного шумоподобного сигнала PA_p = \t%u\n\n",PA_p);
fprintf(f_stud_rezy,"Отношение [раз] сигнал/помеха исходного сигнала реализации SN_ish = \t%u\n\n",SN_ish);
fprintf(f_stud_rezy,"Отношение [дБ] сигнал/помеха исходного сигнала реализации SN_dB_ish = \t%u\n\n",SN_dB_ish);
fprintf(f_stud_rezy,"Мощность случайного шумоподобного AD сигнала в канале PAD_p = \t%u\n\n",PAD_p);
fprintf(f_stud_rezy,"Мощность детерминированного AD сигнала в канале PAD_s = \t%u\n\n",PAD_s);
fprintf(f_stud_rezy,"Отношение [раз] сигнал/помеха для AD в канале SND_kan = \t%u\n\n",SND_kan);
fprintf(f_stud_rezy,"Отношение [дБ] сигнал/помеха для AD в канале SND_dB_kan = \t%u\n\n",SND_dB_kan);
fclose(f_stud_rezy);
% сохраняем графики
print(hf10,strcat(NameFile,'_analita_A','.pdf'), "-dpdf","-S600,400");
print(hf104,strcat(NameFile,'_SP_analit_A','.pdf'), "-dpdf","-S600,400");
print(hf11,strcat(NameFile,'_amp_DPF_A','.pdf'), "-dpdf","-S600,400");
print(hf12,strcat(NameFile,'_phase_DPF_A','.pdf'), "-dpdf","-S600,400");
print(hf111,strcat(NameFile,'_SPM_A','.pdf'), "-dpdf","-S600,400");
%
print(hf14,strcat(NameFile,'_analita_AD','.pdf'), "-dpdf","-S600,400");
print(hf108,strcat(NameFile,'_ampSP_analit_AD','.pdf'), "-dpdf","-S600,400");
print(hf1081,strcat(NameFile,'_phaseSP_analit_AD','.pdf'), "-dpdf","-S600,400");
print(hf15,strcat(NameFile,'_amp_DPF_AD','.pdf'), "-dpdf","-S600,400");
print(hf16,strcat(NameFile,'_phase_DPF_AD','.pdf'), "-dpdf","-S600,400");
print(hf112,strcat(NameFile,'_SPM_AD','.pdf'), "-dpdf","-S600,400");
%
print(hf114,strcat(NameFile,'_grad_Laplac_A','.pdf'), "-dpdf","-S600,400");
print(hf116,strcat(NameFile,'_grad_Laplac_AD','.pdf'), "-dpdf","-S600,400");
else
Rez=0
endif
% конец программы