%% ECE403 reflection coefficient filter
% Author: Roger D. Serwy
% email:  serwy@illinois.edu
% Date: 2010-03-31


% This file shows how to set up a reflection coefficient simulation using
% the bilinear transform.


% The reflection coefficient is:
% Gamma = (Z_load - Z_src) / (Z_load + Z_src)

% Do a reflection filter based on a TX line terminated into
% an inductor.

%    _____________    L
%   (___z0_______(---nnnn---+
%                           |
%                          GND
%


% Set up reflection
% G = (s*L - z0) / (s*L + z0)

% **these questions are for your own understanding**
% if L == 0, then it is a short circuit. What should gamma be?
% if L == inf, then it is an open circuit. What should gamma be?

% Does the simulated response behave as expected?

z0 = 100;       % characteristic impedance
L = 10000;      % inductance value

% correction: the coeffs vector are for decreasing powers of s.
num = [L, -z0];   % numerator of gamma, coeffs of s
den = [L, z0];    % denominator of gamma, coeffs of s
Fs = 100;          % sample rate

[B, A] = bilinear(num, den, Fs);

% compute a response, sample by sample

% The filter command maintains an internal state. This state needs to be
% saved in between each invocation of filter. See "doc filter".


Nsteps = 30;
h = zeros(Nsteps, 1);

drive = zeros(Nsteps, 1);
drive(1) = 1;       % send in an impulse

% initialize filter state to zeros (this is normal)
filt_state = zeros(max(length(A),length(B))-1, 1); % see filter documentation

for n=1:Nsteps
    % grab the input sample
    input_sample = drive(n);

    % compute the output sample, and update the filt_state
    [out_sample, filt_state] = filter(B, A, input_sample, filt_state);

    h(n) = out_sample;    % save the sample
end


figure(1)
subplot(2,1,1)
plot(h)
title('sample-by-sample computation of reflection response')


subplot(2,1,2);
% The filter command can operate on a long vector of samples, and
% automatically keeps track of the internal state. 
h_quick = filter(B, A, drive);  
plot(h_quick);
title('correct result')


figure(2)
freqz(B, A);    % plots the frequency response of a filter. 
title('Frequency response of the reflection filter')




