%% Frequency response demo
% Illustration of the frequency response concept for discrete-time LTI
% systems.

%% Define system
% The system is defined by the difference equation
%
% y(n) = 0.2 x(n) - 0.3 x(n-1) + 0.2 x(n-2) + 1.8 y(n-1) - 0.9 * y(n-2)

clear

b = [0.2 -0.3  0.2];

a = [1 -1.8 0.9];

%% Plot frequency response (radians/sample)
% Use the freqz function to calculate the frequency response.

[H, om] = freqz(b,a);       % om : frequency in units of radians/sample

figure(1)
clf
plot(om, abs(H))
title('|H^f(\omega)|')
xlabel('frequency (\omega) [radians/sample]')
xlim([0 pi])


%% Use cycles/sample
% Frequency axis in units of cycles/sample

f = om / (2*pi);            % f  : frequency in units of cycles/sample

figure(1)
clf
plot(f, abs(H))
title('Frequency response magnitude')
xlabel('frequency (f) [cycles/sample]')

%% The phase response
% Plot the phase of the frequency response.
% Why is there a discontinuity in the phase response?
% What is the size of the phase jump?
% At what frequency is the phase jump?

figure(1)
clf
plot(f, angle(H))
title('Frequency response phase')
xlabel('frequency (f)')

%% Define a cosine signal

om1 = 0.1*pi ;          % Frequency of cosine signal in units of radians/sample
f1 = om1/(2*pi);        % Cycles/sample

n = 0:150;
x = cos(om1*n);         % Input signal

%% Compute output signal
% Use 'filter' to implement the difference equation.

y = filter(b, a, x);    % Output signal

figure(1)
clf
subplot(2, 1, 1)
plot(n, x)
title('Input signal')
xlabel('Time (n)')
ylim([-1.5 1.5])

subplot(2,1,2)
plot(n, y)
title('Output signal')
xlabel('Time (n)')


%% Use 'polyval' to evaluate the frequency response
% Find the value of the frequency response at om1 by evaluating
% the transfer function at z1 = exp(j om1).
% Verify that the value agrees with 'freqz'.

I = sqrt(-1);
z1 = exp(I*om1);

H1 = polyval(b, z1) / polyval(a, z1);    % Value of frequency response at om1

abs(H1)
angle(H1)

%%
figure(1)
clf
subplot(2, 1, 1)
plot(f, abs(H), f1, abs(H1), 'ro', 'markerface', 'red')
title('Frequency response magnitude')
xlabel('frequency (f)')

subplot(2, 1, 2)
plot(f, angle(H), f1, angle(H1), 'ro', 'markerface', 'red')
title('Frequency response phase')
xlabel('frequency (f)')


%% Formula for steady-state output
% Use the value of the frequency response at om1 to compute the 
% steady-state output signal.
% The output signal converges to the steady-state signal.

s = abs(H1) * cos(om1 * n + angle(H1));  % Steady-state output signal

figure(1)
clf
subplot(2, 1, 1)
plot(n, y, n, s, '--');
title('Output and steady-state signals')
xlabel('Time (n)')
legend('output', 'steady-state output')

%% The transient response
% We identify the transient response as the difference
% between the output signal and the steady-state output signal.

figure(1)
subplot(2, 1, 2)
plot(n, s-y);
title('Transient response')
xlabel('Time (n)')


%% Compute output signal for another frequency
% Create cosine signal with frequency om2 = 0.23*pi
% and compute the new output signal.
% Why does the output signal decay to zero? Can this be predicted
% directly from the frequency and/or the pole-zero diagram? [Hint: yes]

om2 = 0.23*pi;          % [radians/sample]
f2 = om2/(2*pi);        % [cycles/sample]
x = cos(om2*n);
y = filter(b, a, x);

figure(1)
clf
subplot(2,1,1)
plot(n, x)
title('Input signal')
xlabel('Time (n)')
ylim([-1.5 1.5])

subplot(2,1,2)
plot(n, y)
title('Output signal')
xlabel('Time (n)')

%% Pole-zero diagram
% The nulls in the frequency response correspond to zeros on the unit
% circle.

figure(1)
clf
zplane(b, a)
title('Pole-zero diagram')

%%

angle(roots(b))/pi


%% Evaluate frequency response at second frequency.
% Find frequency response value at om2.
% [Evaluate transfer function at z2 = exp(j om2)].
% The value is almost zero. The value z2 is (almost) a
% zero of the transfer function.

z2 = exp(1j * om2);
H2 = polyval(b, z2) / polyval(a, z2);    % Value of frequency response at om2
abs(H2)

%% 

figure(1)
clf
subplot(2, 1, 1)
plot(f, abs(H), f2, abs(H2), 'ro', 'markerface', 'red')
title('Frequency response magnitude')
xlabel('frequency (f)')

subplot(2, 1, 2)
plot(f, angle(H), f2, angle(H2), 'ro', 'markerface', 'red')
title('Frequency response phase')
xlabel('frequency (f)')


%% Two frequencies
% The input consists of two cosine signal but the
% output signal consists of one cosine.
% The filter annilhates one of them. (Which one?)

x = 1.0 * cos(om1 * n) + 1.2 * cos(om2 * n);

y = filter(b, a, x);     % Output of filter

figure(1)
subplot(2, 1, 1)
plot(n, x) 
title('Input signal')
xlabel('Time (n)')

subplot(2, 1, 2)
plot(n, y)
title('Output signal')
xlabel('Time (n)')



