close all hidden
clear
clc

set(groot, 'DefaultTextInterpreter', 'LaTeX');
set(groot, 'DefaultAxesTickLabelInterpreter', 'LaTeX');
set(groot, 'DefaultAxesFontName', 'LaTeX');
set(groot, 'DefaultAxesFontSize', 14);
set(groot, 'DefaultLegendInterpreter', 'LaTeX'); 
cmap = lines;

%**************************************************************************
% Basic parameters
%**************************************************************************
f_s = 48e3;          % Sampling freq in Hz

%**************************************************************************
% System characterisation
%**************************************************************************
systemSel = 4;

% All systems are real-valued (partly complex-conj)
switch(systemSel)
    case 1
        % Stable FIR highpass filter
        % Minimal phase
        z = [0.5 + 0.5j, 0.5 - 0.5j]';
        p = [0 0]';
        [b, a] = zp2tf(z, p, 1);
     case 2
        % Instable IIR lowpass/bandpass filter
        % Maximal phase
        z = [-1 + 1i, -1 - 1i]';
        p = [1/sqrt(2) * (1.2 + 1i), 1/sqrt(2) * (1.2 - 1i)]';
        [b, a] = zp2tf(z, p, 1);
    case 3
        % Stable IIR lowpass (Butterworth) filter
        % Mixed phase (implementation-dependent)
        N = 4; 
        [b, a] = butter(N, 0.5);
        [z, p, ~] = tf2zpk(b, a);
    case 4
        % Stable FIR highpass filter
        % Linear phase
        z = [1/(2 * sqrt(2)) *[1 + 1i, 1 - 1i], sqrt(2) *[1 + 1i, 1 - 1i]]';
        p = [0 0 0 0]';
        [b, a] = zp2tf(z, p, 1);       
end

% Get system from coefficients
H = tf(b, a, f_s);            

% Get frequency response from coefficients
[H_f, w] = freqz(b, a, 2^10);
% Norm frequency vector
w = w / pi;

% Get impulse response from system
[h, t_h] = impulse(H);

%**************************************************************************
% Plots
%**************************************************************************
% Evaluate section by section to get plots after each other

figure
% set(gcf, 'Position', [2000 0 1400 1350]) % Set window size

phi = linspace(0, 2*pi, 1000);
x = cos(phi);
y = sin(phi);

% Pole zero plot
subplot(4, 1, 1);
plot(real(z), imag(z), 'LineStyle', 'none', 'Marker', 'o', 'Color', cmap(1, :))
hold on
plot(real(p), imag(p), 'LineStyle', 'none', 'Marker', 'x', 'Color', cmap(1, :))
plot(x, y, '--', 'Color', 'Black')
plot([-4 4], [0 0], '--', 'Color', 'Black')
plot([0 0], [-4 4], '--', 'Color', 'Black')
xlim([-2 2])
ylim([-2 2])
grid on
xlabel('$\mathrm{Re}\big(z\big)$')
ylabel('$\mathrm{Im}\big(z\big)$')
title('Pole-zero plot')

%%
% Magnitude of frequency response
subplot(4, 1, 2)
plot(w, 20*log10(abs(H_f)))
grid on
xlabel('$\frac{\Omega}{\pi}$')
ylabel('$|H\big(e^{j \Omega}\big)|$ / dB')
title('Frequency response - magnitude')

%%
% Impulse response
subplot(4, 1, 3)
stem(t_h, h)
grid on
xlabel('$n \, T_s$ / s')
title('Impulse response')

%%
% Phase of frequency response
subplot(4, 1, 4)
plot(w, rad2deg(unwrap(angle(H_f)))) 
grid on
xlabel('$\frac{\mathrm{\Omega}}{\pi}$')
ylabel('$\phi$ / deg')
title('Frequency response - phase')