% Лабораторная работа 1 ТСиУПД
% Исследование характеристик дискр-квантование треуг. сигнала
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);
%
t=1:1:N; % Формируем сигнал
Fmax=dF/2.4; k=[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);
% массив для ошибок
FrSig=ones(1,length(k));
ErrSig=ones(1,length(k));
ErrSigSS=ones(1,length(k));
ErrSigARSS=ones(1,length(k));
for i=1:length(k)
A=Adelta*sawtooth((k(i)*Fmax*(dT*t)/10),0.5)+Asdv;
Fsignala=k(i)*Fmax/10
%
Amean=mean(A);
Astd=std(A);
%
hf11=figure(); % Рисуем на графике результат (400 точек)
hold on; title('First 400 samples of A','fontsize',14);
plot(t(1:400),A(1:400),'b','linewidth',2); xlabel('samples A, n=1...400'); ylabel('Source A 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 A','fontsize',14);
plot(t(1:100),A(1:100),'b','linewidth',2); xlabel('samples A, n=1...100'); ylabel('Source A 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;
%

% ОПСЬ....
%Nk=8; % Задаем разрядность квантования
% нормируем и квантуем по Nk
ZSig=A;
ZSig_max=max(ZSig);
ZSig_min=min(ZSig);
Nk=ceil(log2((ZSig_max-ZSig_min)/dA));
% определяем поправочн. коэфф.
%Kpmax=abs(ZSig_max)/(2^(Nk-1)-1);
%Kpmin=abs(ZSig_min)/(2^(Nk-1)-1);
%Kpn=(Kpmax+Kpmin)/2;
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); % опред. пост. составл.
%Kvar=MZSig_1/MZSig_1NKP;
% запись в файл данные с квантованием по уровню двоичн. 1...(2^Nk-1);
%dlmwrite('ZSig_NK.txt',ZSig_NK,',');
% запись в файл данные с квантованием по уровню двуполярные
%dlmwrite('ZSig_NKP.txt',ZSig_NK_Polar,',');
% восстановление исходного сигнала;
ZSigVOST=((ZSig_NK-1).*Kpr).+ZSig_min;
hf61=figure();
subplot(311);
plot(t,ZSig,"b;signal;",
     t,((ZSig_NK-1).*Kpr).+ZSig_min,";signal-Norm-Kvant;");
ylabel ("Amplitude");
xlabel ("Time (ms)");
subplot (312);
plot(t,sort((ZSig)),'b;sign-sort;',
     t,sort(((ZSig_NK-1).*Kpr).+ZSig_min), ';sign-NKsort;');
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=meansq(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;otnosit oshibka kvant;",
     t,Es.*Zoshsr,";energiya oshibki kvant;");

% Блок фильтрации(сглаживания) исходного сигнала
ZSigW=ZSigVOST;
% Максимальный период (для определения нижней границы частоты условного сигнала)
Tmax_ZSigW=N/(2^Nk); Fmin_ZSigW=1/Tmax_ZSigW;
% Находим функцию скользящего среднего ZSigW-аналогично FIR-фильтрации без весов
wZSW=(fix(Tmax_ZSigW/5)); % задаем размер окна КИХ фильтра
if(mod(wZSW,2)==0)
  wZSW=wZSW+1;
endif
hwZSW=(wZSW-1)/2; % интервал влево и вправо от середины окна
MZSwW=zeros(1,N); MZSwW(1)=ZSigW(1);
for jW=2:N-1 % Задаем размер индексного окна в массиве
      MZSanW=0; zwW=0;
     if(jW<=hwZSW)
     k1W=1; k2W=jW; zwW=k2W-k1W+1;
     elseif(jW+hwZSW>N)
     k1W=jW-hwZSW; k2W=jW; zwW=k2W-k1W+1;
     else
     k1W=jW-hwZSW; k2W=jW+hwZSW; zwW=wZSW;
     endif
     for lW=k1W:k2W
       MZSanW=MZSanW+ZSigW(lW);
     endfor
  MZSwW(jW)=MZSanW/zwW;
endfor
MZSwW(N)=ZSigW(N);
M_ZSigW=MZSwW;

% Выполним фильтрацию ZSigW-аналогично авторегресии(АР) со СС, но с весами
% Определим коэфф. IIR-фильтра БИХ
[b,a]=butter(6,0.5);
ARM_ZSigW=filtfilt(b,a,ZSigW);
% Определяем ошибки восстановления сигнала с учетом фильтрации FIR и IIR
ZoshfilSS=A.-M_ZSigW;
ZoshfilARSS=A.-ARM_ZSigW;
Izm_ZoshfilSS=meansq(ZoshfilSS);
Izm_ZoshfilARSS=meansq(ZoshfilARSS);
hf91=figure();
plot(t,A,"b;signal;",
     t,(M_ZSigW),";signal vost.SS(FIR);");
ylabel ("Amplitude");
xlabel ("Time (ms)");
hf911=figure();
plot(t,A,"b;signal;",
     t,(ARM_ZSigW),";signal vost.ARSS(IIR);");
ylabel ("Amplitude");
xlabel ("Time (ms)");
hf101=figure();    
plot(t,ZoshfilSS,"b;oshibka kvant SS(FIR);",
     t,Es.*Izm_ZoshfilSS,";energiya oshibki kvant SS;");
hf102=figure();    
plot(t,ZoshfilARSS,"b;oshibka kvant ARSS(IIR);",
     t,Es.*Izm_ZoshfilARSS,";energiya oshibki kvant ARSS;");
% 
FrSig(i)=FrSig(i)*Fsignala;
ErrSig(i)=ErrSig(i)*Zoshsr;
ErrSigSS(i)=ErrSigSS(i)*Izm_ZoshfilSS;
ErrSigARSS(i)=ErrSigARSS(i)*Izm_ZoshfilARSS;
% Блок сохранения результатов
NameFile=strcat('Lab_TSiUPD_1_TriangW_','Fsig_',num2str(round(Fsignala),5),'Hz','.txt');
RezS=0;
disp('Желаете сохранить рез-ты в файлы Lab_TSiUPD_1_TriangW.txt и графики в ***.pdf?');
disp('[ 1 - ДА, желаю; 0 - НЕТ, не желаю (по умолчанию) ]'); 
Rez=input('Rez=');
if Rez<=0;
  Rez=0;
elseif Rez>0;
  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,"  \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, "Laba_TSiUPD_1_TriangW_Graf_400N_M_SKO.pdf", "-dpdf","-S600,300");
%print(hf21, "Laba_TSiUPD_1_TriangW_Graf_100N_M_SKO.pdf", "-dpdf","-S600,300");
print(hf61,strcat('Lab_TSiUPD_1_TriangW_','Fsig_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,400");
print(hf71,strcat('Lab_TSiUPD_1_TriangW_','Graf_Sign_AbsolErr_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,400");
print(hf81,strcat('Lab_TSiUPD_1_TriangW_','Graf_OtnErr_EnjErr_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,600");
print(hf101,strcat('Lab_TSiUPD_1_TriangW_','Graf_Err_kvant_FIR_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,300");
print(hf102,strcat('Lab_TSiUPD_1_TriangW_','Graf_Err_kvant_IIR_',num2str(round(Fsignala),5),'Hz','.pdf'), "-dpdf","-S600,300");  
endif
close all;
continue
endfor
hf103=figure();    
plot(FrSig,ErrSig,"b;Energiya oshibki kvant;");
xlabel('frequency sign.A, [Hz]'); ylabel('Norm. energiya err.A');
hf104=figure();    
plot(FrSig,ErrSigSS,"b;Energiya oshibki kvant SS;",
     FrSig,ErrSigARSS,";Energiya oshibki kvant ARSS;"); 
xlabel('frequency sign.A, [Hz]'); ylabel('Norm. energiya err.A');     
print(hf103, "Laba_TSiUPD_1_TriangW_Graf_Energii_Err_Lin.pdf", "-dpdf","-S600,300");
print(hf104, "Laba_TSiUPD_1_TriangW_Graf_Energii_Err_FIR_IIR.pdf", "-dpdf","-S600,300");
% конец программы


