% Практическая работа 1 ТСиУПД Источники случ.сообщений пар-ры, энтр. 
% Для ист-ка B - нормальное распределение
clc;
close all;
clear all;
% Блок ввода данных
disp('Введите номер варианта [1 . . . 30]');
n=input('n=');
if n<=0;
  n=1;
elseif n>30;
  n=30;
end
disp('Исходные данные варианта');disp(n);
disp('Введите Количество отсчетов N [1000...6000]');
N=input('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('Количество отсчетов N = ');disp(N);
disp('Частота дискретизации dF = ');disp(dF);
disp('МИН знач.пар.ист.сообщ. AMin = ');disp(AMin);
disp('МАКС знач.пар.ист.сообщ. AMax = ');disp(AMax);
disp('Интервал дискр.пар.ист.сообщ. dA = ');disp(dA);
%N=2500;  % Количество отсчетов
%dF=1000;
%AMin=-1.7;
%AMax=19.2;
%dA=0.2;
%
t=1:N; % Формируем источник сообщений B
ARan=randn(1,N); dT=1/dF;
MaxARan=max(ARan); MinARan=min(ARan);
ARan1=(ARan.-MinARan)./(MaxARan-MinARan);
A=(ARan1.*abs(AMax-AMin)).+AMin; 
Es=ones(N);
%
Amean=mean(A);
Astd=std(A);
%
hf11=figure(); % Рисуем на графике результат (400 точек)
hold on; title('First 400 samples of message source B','fontsize',14);
plot(t(1:400),A(1:400),'b','linewidth',2); xlabel('samples A, n=1...400'); ylabel('Source B parameter value');
plot(t(1:400),Es(1:400).*Amean,'r','linewidth',3);
plot(t(1:400),Es(1:400).*Amean.+Astd,'g','linewidth',2);
plot(t(1:400),Es(1:400).*Amean.-Astd,'g','linewidth',2);
grid on;
hf21=figure(); % Рисуем на графике фрагмент (100 точек)
hold on; title('First 100 samples of message source B','fontsize',14);
plot(t(1:100),A(1:100),'b','linewidth',2); xlabel('samples A, n=1...100'); ylabel('Source B parameter value');
plot(t(1:100),Es(1:100).*Amean,'r','linewidth',3);
plot(t(1:100),Es(1:100).*Amean.+Astd,'g','linewidth',3);
plot(t(1:100),Es(1:100).*Amean.-Astd,'g','linewidth',3);
grid on;
%
% Строим гистограмму данных и функцию распределения
Stat_A_1=[min(A), median(A), max(A), mean(A), std(A), skewness(A), kurtosis(A), var(A)]; % Модуль статистики по А
Stat_A_2=quantile(A);
K_1_hist=ceil(1+3.32*log10(N)); % Число колонок гистограммы (правило Штюргеса ГОСТ)
Mhist=K_1_hist; % Число колонок гистограммы::=> выбрать(но, не более N/5-рек.ГОСТ)
[NZSig, IDZSig]=hist(A,Mhist); Pv_NZSig=NZSig./N;
[NZSig3, IDZSig3]=hist(A,(3*Mhist)); Pv_NZSig3=NZSig3./N;% гистограмма на 3*K-колонок
% Считаем эмпирическую функцию распределения(ЭФР)- сглаженную по гистограмме 
FZSig(1)=0;
for l=2:length(NZSig3);
FZSig(l)=FZSig(l-1)+NZSig3(l);
endfor
FxZSig=[AMin, IDZSig3, AMax];
FyZSig=[0, (FZSig(1)+1)/N, FZSig(2:1:length(NZSig3))./N, 1];
% или второй вариант ЭФР- полную по ВХ
FxZSig_1=sort(A); FyZSig_1=1/N:1/N:1;
% Рисуем на графике гистограмму данных
hf31=figure(); hold on; title('Message source histogram B','fontsize',14);
bar(IDZSig,NZSig,'c'); plot(IDZSig,NZSig,'r','linewidth',2); xlabel('A'); ylabel('K');
bar(IDZSig3,NZSig3,'g'); plot(IDZSig3,NZSig3,'b','linewidth',2); grid on;
hf41=figure(); hold on; title('The probability density function of the source of messages B','fontsize',14);
bar(IDZSig,Pv_NZSig,'c'); plot(IDZSig,Pv_NZSig,'r','linewidth',2); xlabel('A'); ylabel('P');
bar(IDZSig3,Pv_NZSig3,'g'); plot(IDZSig3,Pv_NZSig3,'b','linewidth',2); grid on;
% Рисуем на графике функцию распределения
hf51=figure(); hold on; title('Distribution function of message source B','fontsize',14); xlabel('A'); ylabel('P');
plot(FxZSig,FyZSig,'g','linewidth',3,FxZSig_1,FyZSig_1,'b','linewidth',2); grid on;

% ОПСЬ....
% Определяемем разрядность квантования Nk
% нормируем и квантуем по 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-но полярн. полож.
%
% рисуем графики сравнения квант.сигнал, амплитуд. распр., ф-ю ошибки.;  
hf61=figure();
subplot(311);
plot(t,ZSig,"b;signal;",
     t,((ZSig_NK-1).*Kpr).+ZSig_min,";signal-Diskr-Kvant;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
subplot (312);
plot(t,sort((ZSig)),'b;sign-sort;',
     t,sort(((ZSig_NK-1).*Kpr).+ZSig_min), ';sign-DKsort;');
xlabel ("Rang");
ylabel ("Amplitude");
subplot (313);
plot(t,sort(ZSig)-sort(((ZSig_NK-1).*Kpr).+ZSig_min),'b;sign-error;');
xlabel ("Rang");
ylabel ("Error");
% Ошибка при восстановлении
Zosh=ZSig.-(((ZSig_NK-1)*Kpr).+ZSig_min);
Zoshsr=mean(Zosh);
Zoshproc=(abs(Zosh)./abs(ZSig));
hf71=figure();
plot(t,ZSig,"b;signal;",
     t,(Zosh*1),";signal oshibki;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
hf81=figure();    
plot(t,Zoshproc,"b;oshibka kvant %;");
% определяем энтропию источника, как для непрерывного процесса 
Hnrav=ceil(log2(sqrt(12*e)*Astd)); Hdpr=ceil(log2((AMax-AMin)/dA));
% определяем инф. энтропию источника, как для дискр. процесса (1-го дискрета)
Hdrav=log2(2^Nk);
% определяем кол-во байт для одного дискрета и битовую скорость передачи ист-ка
Kbt=ceil(Hdrav/8); Bbit=Kbt*8/dT;
% Блок сохранения результатов
Rez=0;
disp('Желаете сохранить рез-ты в файлы Praktic_1_istB.txt и графики в ***.pdf?');
disp('[ 1 - ДА, желаю; 0 - НЕТ, не желаю (по умолчанию) ]'); 
Rez=input('Rez=');
if Rez<=0;
  Rez=0;
elseif Rez>0;
  Rez=1;
% Запись результатов в файл 
f_stud_rezy=fopen('Praktic_1_istB_TSiUPD.txt','wt');
fprintf(f_stud_rezy,"Исходные данные варианта \t%u\n\n",n);
fprintf(f_stud_rezy,"Количество отсчетов N = \t%u\n\n",N);
fprintf(f_stud_rezy,"Частота дискретизации 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,"  \t\n\n");
fprintf(f_stud_rezy,"Результаты вычислений параметров ист.сообщ. B \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,"Гистограмма построена по правилу Штюргеса(ГОСТ)+3K  K = \t%u\n\n",Mhist);
fprintf(f_stud_rezy,"Энтропия непр.ист.сообщ. B Hnrav = \t%u\n\n",Hnrav);
fprintf(f_stud_rezy,"Энтропия информ.ист.сообщ. B Hdrav = \t%u\n\n",Hdrav);
fprintf(f_stud_rezy,"Кол-во уровней дискрет.ист.сообщ. B = \t%u\n\n",2^Hdrav);
fprintf(f_stud_rezy,"Битовая скорость передачи ист.сообщ. B Bbit = \t%u\n\n",Bbit);
fclose(f_stud_rezy);
% сохраняем графики
print(hf11, "Praktic_1_istB_Graf_400N_M_SKO.pdf", "-dpdf","-S600,300");
print(hf21, "Praktic_1_istB_Graf_100N_M_SKO.pdf", "-dpdf","-S600,300");
print(hf41, "Praktic_1_istB_Graf_Plotn_ver.pdf", "-dpdf","-S600,400");
print(hf51, "Praktic_1_istB_Graf_Funk_rasp.pdf", "-dpdf","-S600,400");
print(hf61, "Praktic_1_istB_Graf_Diskr_ist.pdf", "-dpdf","-S600,600");
print(hf81, "Praktic_1_istB_Graf_Err_kvant.pdf", "-dpdf","-S600,300");  
end
% конец программы

