%% Matlab demo: Frequency response
% Illustrate frequency response concept and equations

%% Start

close all
clear

%% Define a system

% The system y(n) = 0.2 x(n) + 0.3 x(n-1) + 0.8 y(n-1)
% is represented by its coefficient vectors 'a' and 'b':
b = [0.2 0.3];
a = [1 -0.8];

%% Define a cosine

om1 = 0.09*pi;   % radians/sample

n = 0:100;
x = cos(om1 * n);

figure(1)
stem(n, x, '.')
ylim([-2 2])

% The frequency om1 is f1 = 0.045 cycles/sample.  f1 = om1/(2*pi).
% So, one cycle is 1/0.045 = 22.2 samples,
% which can be seen in the plot.

%% Compute system output

y = filter(b, a, x);

figure(1)
stem(n, y, '.')

%% Evaluate H(om1)

z1 = exp(1j * om1);
Hom1 = (0.2*z1 + 0.3)/(z1 - 0.8);   % Value of frequency response at om1

%%

abs(Hom1)       % This is |H(om1)|

%%

angle(Hom1)      % This is the angle of H(om1)


%% Formula for output
% The output, as predicted by the frequency response formula, is:

p = abs(Hom1) * cos(om1 * n + angle(Hom1));

figure(1)
stem(n, p, '.')

%%
% Plot both the true output and the predicted output signal to see if they
% agree

figure(1)
plot(n, y, n, p, '--')
legend('True output', 'Predicted output')

% Note: the predicted output is correct except for the transient.
% The formula is correct once _steady_state_ is reached.

%% Use freqz
% The frequency response can be computed using the 'freqz' command.

[H, om] = freqz(b, a);

figure(1)
plot(om, abs(H))
xlabel('Frequency (\omega) in radians/sample')

%%
% The value we computed for H(om1) should be consistent.
% We can check that H(om1) lies on this plot.

figure(1)
plot(om, abs(H), om1, abs(Hom1), 'ro')
xlabel('Frequency (\omega) in radians/sample')

%%
% Plot using frequency in units of cycles/sample

figure(1)
plot(om/(2*pi), abs(H), om1/(2*pi), abs(Hom1), 'ro')
xlabel('Freqency (f) in cycles/sample')
title('Frequency response |H(f)|')

%% Pole/zero plot

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

%% Two frequencies
% Create a signal consisting of two cosines
% and compute the output of the system.

om2 = 0.6*pi;

x = cos(om1 * n) + 0.5 * cos(om2 * n);

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

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

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

%% Evaluate H(om2)
% H(om1) is found already above. Now find H(om2)

z2 = exp(1j * om2);
Hom2 = (0.2*z2 + 0.3)/(z2 - 0.8);   % Value of frequency response at om2

%%
% Compare H(om1) and H(om2)

abs(Hom1)
% Note: |H(om1)| > 1 so the first cosine is _amplified_ by the system.

%%

abs(Hom2)
% Note: |H(om2)| < 1 so the second cosine is _attenuated_ by the system.

%% 
% Predict the output using the frequency response formula

p = abs(Hom1) * cos(om1 * n + angle(Hom1)) + 0.5 * abs(Hom2) * cos(om2 * n + angle(Hom2));

figure(1)
clf
subplot(2,1,1)
plot(n, p, n, y, '--')
legend('True output', 'Predicted output')

% The predicted output is valid in the steady state.

%% The transient?
% Q: How long does it take for the system to reach steady state?  
%
% A: As long as it takes for the impulse response to diminish.
% Which depends on how close the pole is to the unit circle.
% The closer the pole is to the unit circle, the longer it takes
% for the impulse response to decay to zero.






