% demo.m 
% Script for the example in Section 5.1.2, pp 278
% Author: Dr. David Pan, for EE 648, UAH

clear all
close all

% Generate the filter
b1 = [1 -0.98*exp(j*0.8*pi)];
b2 = [1 -0.98*exp(-j*0.8*pi)];

a1 = [1 -0.8*exp(j*0.4*pi)];
a2 = [1 -0.8*exp(-j*0.4*pi)];

H1_b = conv(b1, b2);
H1_a = conv(a1, a2);

%{
figure;
freqz(H1_b, H1_a)
figure;
zplane(H1_b, H1_a)
figure;
[h, w] = freqz(H1_b, H1_a);
plot(w/pi,abs(h));grid
%}

H2_b = 1;
H2_a = 1;

for k = 1: 4
    c = 0.95*exp(j*(0.15*pi + 0.02*pi*k));
    cc = conj(c);
    
    const = cc * c; 
    b1 = [1 -1/cc];
    b2 = [1 -1/c];
    b_c = const*conv(b1, b2);
    
    a1 = [1 -c];
    a2 = [1 -cc];
    a_c = conv(a1, a2);
    
    % Take the square
    bk = conv(b_c, b_c);
    ak = conv(a_c, a_c);

    H2_b = conv(H2_b, bk);
    H2_a = conv(H2_a, ak);
end

%{
figure;
freqz(H2_b, H2_a)
figure;
zplane(H2_b, H2_a)
figure;
[h, w] = freqz(H2_b, H2_a);
plot(w/pi,abs(h));grid
%}

b = conv(H1_b, H2_b);
a = conv(H1_a, H2_a);

% Plot the frequency response
freqz(b, a);

% Zero-pole plot
figure;
zplane(b, a);

% Group delay
figure;
grpdelay(b, a);

% Magnitude without dB
figure;
[h, w] = freqz(b, a);
plot(w/pi,abs(h));grid

% Generate the input sequence
M = 60;

x = zeros(1, 300);
for n = 1: 300
    s = 0;
    if (n >= 0 && n <= M)
        w = 0.54 - 0.46*cos(2*pi*n/M);
    else
        w = 0;
    end
    s = s + w * cos(0.8*pi*n + pi/5);
    
    if ((n-M-1) >= 0 && (n-M-1) <= M)
        w = 0.54 - 0.46*cos(2*pi*(n-M-1)/M);
    else
        w = 0;
    end
    s = s + w * cos(0.2*pi*(n-M-1)); 
    
    if ((n-2*M-2) >= 0 && (n-2*M-2) <= M)
        w = 0.54 - 0.46*cos(2*pi*(n-2*M-2)/M);
    else
        w = 0;
    end
    s = s + w * cos(0.4*pi*(n-2*M-2)- pi/2);
    
    x(n) = s;
end
figure;
plot(x)
figure;
freqz(x,1);
    
% Apply the filter on the input sequence
y = filter(b,a,x);
figure;
plot(y)


    
    