clear, close all, clc;
Fs = 2000;            % Frekvence vzorkování                    
Fa = 100;             % Frekvence signálu
Amp = 20;             % Amplituda signálu
Amp = Amp/2;
%T = 1/Fs;             % Vzorkovací perioda       
Ts = 1/Fs;            % Vzorkovací perioda
L = 3000;             % Délka signálu
t = (0:L-1)*Ts;        % Časový vektor
AD = 2^16;            % AD rozlišení

%%%%%%%%%%%%%% Sampling bez chyby

S = Amp*sin(2*pi*Fa*t);    % Analogový vstup
    p = min(S);
    S = S - p;
    
%%%%%%%%%%%%%% Sampling s chybou

for i  = 1:L
    
    chyba = randi([1 75])/100;
    %chyba = 75/100;
    vyskyt = randi([0 8]);
    if vyskyt == 0
        te(i)=(i-1)*Ts + Ts*chyba;
    elseif vyskyt == 1
        te(i)=(i-1)*Ts - Ts*chyba;
    else
        te(i) = (i-1)*Ts;
    end
end

Sae = Amp*sin(2*pi*Fa*te);
    p = min(Sae);
    Sae = Sae - p;
    
%%%%%%%%%%% Kvantování
%%%%%%%%%%% Urovně
          
rozptyl = max(S)-min(S);
j=1;
for i=min(S):rozptyl/(AD-1):max(S)
    Kvant(j)=i;
    j=j+1;
end

%%%%%%%%%%% Přiřazení úrovní
            % signal Sa
          
for i = 1:length(S)
   for j =1:length(Kvant)
       if S(i) - Kvant(j) <= 0
           if j<=1
               j=2;
           end
           if j == 1024
               Sa(i) = Kvant(j);
           else
               Sa(i) = Kvant(j-1);
           end
           break;
       end
   end
end

            % signal Se
            
for i = 1:length(Sae)
   for j =1:length(Kvant)
       if Sae(i) - Kvant(j) <= 0
           if j<=1
               j=2;
           end
           if j == 1024
               Sae(i) = Kvant(j);
           else
               Sae(i) = Kvant(j-1);
           end
           break;
       end
   end
end
%%%%%%%%%%% Vykreslení do grafu
j=1;
for i = 1:length(Sa)-1
    gSa(j)=Sa(i);
    gSae(j)=Sae(i);
    gx(j)=t(i);
    
    j=j+1;
  
    gSa(j)=Sa(i);
    gSae(j)=Sae(i);
    
    if i==length(Sa)
        gx(j)=t(j-1)+t(1);
    else
        gx(j)=t(i+1);
    end
    j=j+1;
end
figure(4);
hold on;

plot(gx,gSae,'r-');
plot(gx,gSa,'k-');
axis([0 2/Fa min(Sae) max(Sae)+1]);
title('Digitalizovaný signál')
xlabel('t (s)')
ylabel('I (mA)')
legend('Yp - signál s nepřesností','Y - signál digitální')

%%%%%%%%%% FFT

figure(1);
Y = fft(Sa);
P2 = abs(Y/L*(50/24.96));
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
hold on;

plot(f,P1); 
%axis([0 0.035 -1 1])
title('Frekvenční spektrum')
xlabel('f (Hz)')
ylabel('A (dB)')

figure(2);
Y = fft(Sa);
P2 = abs(Y/L);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
AMP2 =  sort(P1);
k=max(S)/AMP2(length(AMP2)-1);
    Y = fft(Sa);
    P2 = abs(Y/L*k);
    P1 = P2(1:L/2+1);
    P1(2:end-1) = 2*P1(2:end-1);
    f = Fs*(0:(L/2))/L;
    AMP2 =  sort(P1);
hold on;
FFTorP1 = P1;
FFTorf = f;
plot(f,P1); 
title('Frekvenční spektrum')
xlabel('f (Hz)')
ylabel('I (mA)')
for i = 1:length(f)
    if f(i)>Fa
        index = i-1;
        break;
    end
end


figure(3);


axis([0 max(f) 0 Amp/10]);
Y = fft(Sae);
P2 = abs(Y/L*k);
P1 = P2(1:L/2+1);
P1(2:end-1) = 2*P1(2:end-1);
f = Fs*(0:(L/2))/L;
hold on;
plot(f,P1,'r-'); 
plot(FFTorf,FFTorP1,'k-');
title('Frekvenční spektrum')
xlabel('f (Hz)')
ylabel('U (V)')
(legend('Yp - signál s nepřesností','Y - signál digitální'))
AMP = sort(P1);
avg = mean(AMP(1:length(AMP)-1));
mod = mode(round(AMP,3));
pole = [P1(index) round(max(P1(index+30:length(P1)-10)),3) round(avg,3) round(mod,3) ]
%close all;
% %%%%% wav write

% audiowrite('Sa.wav',Sa,Fs);
% audiowrite('Sae.wav',Sae,Fs);