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

clear; close all;

%% 
K = 10;
T1 = 1;
T2 = 0.1;

% K/(T1p + 1)/(T2p + 1)

N = 1000;
W10min = -1;
W10max = 3;

%% Построение эталонной модели

% filter = tf([K], [T1*T2, T1+T2, 1]);
% figure('name', 'АЧХ и ФЧХ непрерывного звена');
% bode(filter);
% grid;

%% АЧХ и ФЧХ по ПФ
w = logspace(W10min, W10max, N);

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

LPR = -atan( ((T1 + T2)*w)./(1 - T1*T2*w.^2) ); 
LPR( LPR > 0 ) = LPR( LPR > 0) - pi;
LPR = LPR * 180/pi;

%% Расчёт ЧХ цифрового звена

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

duration = 300;              % Длительность сигнала.
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/Td;
        kx = K;
        
        yf(j) = k1*yf(j-1) + k2*yf(j-2) + kx*x(j);
        yf(j) = yf(j) / ( T1*T2/Td/Td + (T1 + T2)/Td + 1 ) ;
    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) );   

%     plot(xs);
%     hold all;
%     plot(yf, '--');
%     hold off

end

%%
figure('name', 'Самодельные характеристики');
hold all;

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

% subplot(2, 2, 1);
hold all
semilogx(w, Ldf, '*');
ylabel('L, дБ');
xlabel('{\omega}, Гц');
title('АЧХ цифрового звена');
legend('ЛАЧХ непрерывного звена', 'ЛАЧХ дискретного звена');

subplot(2, 1, 2);
semilogx(w, LPR);
title('ЛФЧХ');
ylabel('{\phi}, {^\circ}');
xlabel('{\omega}, Гц');
grid;