%%
% 1. По передаточной функции построить теоретические графикию
% ЛАЧХ и ЛФЧХ.
% 2. Получить формулы АЧХ 
% 3. Построить на вычислителе графики ЛАЧХ и ЛФЧХ по полученным
% формулам. 
% 4. Рассчитать разностную формулу цифрового звена, реализующего
% заданную передаточную характеристику, рассчиать минимальный период
% тактирования.
% Согласно разностной формуле построить на компьютере ЛАЧХ цифрового звена.

clear; close all;

%% 
K = 10;
T1 = 0.001;
ksi = 0.4;



N = 200;
W10min = 0;
W10max = 5;


%% АЧХ по ПФ

w = logspace(W10min, W10max, N);

LMR = K./sqrt( (1 - T1.^2.*(w.^2)).^2 + 4*ksi.^2*(T1).^2.*w.^2 );
LMR = 20*log10(LMR);


%% Расчет ЦЗ
Td = 2*pi/10 * T1 * 0.2;          % Период дискретизации.
fD = 1/Td ;                 % Частота дискретизации.

duration = 200;              % Длительность сигнала.
nSamples = duration * fD;   % Число отсчётов в сигнале. 

Ldf = zeros(size(w));       % Точки АЧХ цифрового фильтра. 
for i = 1:length(w)-1
    s = (1:nSamples);
    x = sin( w(i)/fD *s );   % Синус текущей частоты. 
    
    % Цифровая фильтрация
    yf = x;
    for j = 3:length(x)
%         k1 = (2*T1*T2)/Td/Td + (T1*T2)/Td ;
%         k2 = -T1*T2/Td;
%         kx1 = K/Td;
%         kx2 = -K/Td;
         k1 = 2*T1*T1*Td+2*ksi*T1 ;
         k2 = -T1*T1;
         kx1 = Td*Td;
      %   kx2 = -K/Td;
        
        yf(j) = k1*yf(j-1) + K*k2*yf(j-2) + K*kx1*x(j);
        yf(j) = yf(j) / T1*T1*Td+2*ksi*T1 ;
    end
    
    margin = floor(0.1*length(x));
    xs = x(margin:end);
    yfs = yf(margin:end);
    
  
    Ldf(i) = 20*log10(pi/2*sum(abs(yfs))/length(yfs) ); 
    
end

semilogx(w, LMR);
title('ЛАЧХ');
ylabel('L, дБ');
xlabel('{\omega}, Гц');
grid;

hold all;

semilogx(w, Ldf, 'r--');
ylabel('L, дБ');
xlabel('{\omega}, Гц');
title('АЧХ цифрового звена');
legend('ЛАЧХ 1 звена', 'ЛАЧХ 2 звена');