%% ECE403 Homework 6 startup code
% Author: Roger D. Serwy
% email:  serwy@illinois.edu
% Date: 2010-03-31; 2013-03-23(jba)

% This code will simulate an impulse traveling into a tube terminated in a
% short circuit (Z_load = 0). If you want to slow down the simulation,
% uncomment the "pause" instruction. This requires you to press a key
% between each simulation step.

% Each discrete section along the ear canal has its own 
% left and right boundaries. The right boundary of index 1 is the 
% left boundary of index 2, etc...

% A depiction of two sections along a TX line
%
%       |                |                |
% INPUT |          U_p(1)|        ->U_p(2)|
%       |                |                |
%       |U_n(1)<-        |U_n(2)          | 
%       |                |                |
%       B1               B2               B3
%  

% For each simulation time step, a plane wave will travel a distance dx.
% Each section has a length dx.

% At the next time step, U_p(2) gets U_p(1)'s value.
% Likewise, U_n(1) gets U_n(2)'s value.

% U_p(2) is at B3. The wave may reflect from the boundary, if the boundary conditions
% require it, and this reflect wave travels to U_n(2) in one simulation step.

% At boundary B2, the net volume velocity is  U = U_p(1) - U_n(2). 
% The total pressure is  P = z0 * (U_p(1) + U_n(2)).

% The computed next state buffers, denoted by U_p_next and U_n_next
% contain the shifted values. Thus the total U at B2 can 
% be expressed as U = U_p(1) - U_p_next(1)

% This file uses "Cell Mode". In between the %% delimiters, press
% Ctrl+Enter to run the code. See matlab help index, under "cell mode"

%%
% physical parameters (pick you own units if you want, but be consistent)
z0 = 50;   % characteristic impedance of the TX line
c = 100;    % m/s
L = 1;      % m

N = 4;     % number of sections
dx = L / N;   % length of each section of the ear canal
T = (L / c) / N;   % time for wave to travel whole length, div by # of sections
Fs = 1 / T;         % sample rate required

xpos = (1/2 + [0:N-1]) * dx;   % center position of each TX line section


%% A simulation of a wave propagating, using vol. velocity variables
figure(1);  
clf();          % clear the figure

U_p = zeros(1, N);      % forward-traveling wave
U_n = zeros(1, N);      % backward-traveling wave

U_p_next = zeros(1, N); % temporary buffer 
U_n_next = zeros(1, N); % temporary buffer

Nsteps = 12;            % Number of simulation steps to make
t = [0:Nsteps-1] / Fs;  % physical time, given a simulation step index

P_save = zeros(Nsteps, 1);   % a buffer for saving simulation results
U_save = zeros(Nsteps, 1);   

Gamma_R = -1;       % reflection coefficient at the end

%drive = cos(2*pi*t*50);    % drive with a cosine wave (don't use for impulse response)

drive = zeros(Nsteps, 1);   % driving function for the tube, [volume velocity]
drive(1) = 1;       


for n=1:Nsteps

    % compute next state of the TX line
    U_p_next = [0, U_p(1:end-1)];   % right-shift the forward-traveling wave
    U_n_next = [U_n(2:end), 0];     % left-shift the backward-traveling wave

    % grab wave fronts that will hit left and right boundaries
    U_out_R = U_p(end);       % incident wave pressure at the right boundary
    U_out_L = U_n(1);         % incident wave pressure at the left boundary

    % compute reflections
    U_ref_R = Gamma_R * U_out_R;   % reflected wave from right side

    % place reflections into other direction
    U_n_next(end) = U_ref_R;  % place reflected wave into backward path
    %U_p_next(1) = U_ref_L;   % commented out on purpose.

    % place driving wavefront into the system
    U_p_next(1) = drive(n);          % inbound wave front into the tube

if (1)  % change to (0) if you don't want simulation plots
    % display current state of the TX line

    subplot(2,2,1);
    stem(U_p);   % current vector in blue
    xlim([0,N+1]); ylim([-2, 2])  
    title('forward vector') 
    ylabel('volume velocity')

    
    subplot(2,2,3);
    stem(U_n);
    xlim([0,N+1]); ylim([-2, 2])
    title('backward vector')
    xlabel('index')
    ylabel('volume velocity')
    
    
    subplot(2,2,2); 
    stem(xpos + dx/2, U_p - U_n_next); 
    % draw lines representing the right side of each section
    for i=1:N
        line((i) * dx * [1, 1], [-2, 2], 'color', [.5, .5, .5]);
        text((i-1) * dx + dx/2, -1, sprintf('%i', i));  % label each segment
    end 
    xlim([-dx, L+dx]); ylim([-2, 2]);
    title('total U at right side of each section')
    xlabel('position (m)')

    
    subplot(2,2,4);
    stem(xpos - dx/2, U_p_next - U_n); 

    % draw lines representing the left side of each section
    for i=1:N
        line((i-1) * dx * [1, 1], [-2, 2], 'color', [.5, .5, .5]);
        text((i-1) * dx + dx/2, -1, sprintf('%i', i));  % label each segment
    end 
    
    xlim([-dx, L+dx]); ylim([-2, 2]);
    title('total U at left side of each section')
    drawnow(); 

    %pause();    % uncomment to require key presses between sim steps.
                 % (You must have the figure window in focus to accept key presses.)
end

    % save pressure and volume velocity of left side of canal
    P_save(n) = z0 * (U_p_next(1) + U_n(1));
    U_save(n) = (U_p_next(1) - U_n(1));

    % clock it
    U_p = U_p_next;
    U_n = U_n_next;
end


%% plot the input and the measured response

figure(2)
subplot(3,1,1);
stem(t, drive)
xlim([t(1) - T, t(end)+T])
title('Driving vol. vel. waveform at left of ear canal')
ylabel('volume velocity')

subplot(3,1,2);
stem(t, U_save);
xlim([t(1) - T, t(end)+T])
title('Measured volume velocity at left of ear canal')
ylabel('volume velocity')

subplot(3,1,3);
stem(t, P_save);
xlim([t(1) - T, t(end)+T])
title('Measured pressure at of ear canal')
ylabel('pressure')
xlabel('time')

%% measure the input impedance to the tube
figure(3)
% compare to theory from Wikipedia
% see http://en.wikipedia.org/wiki/Transmission_line#Input_impedance_of_lossless_transmission_line

% We need to measure the input impedance of the tube, which is taken
% as the frequency domain voltage divided by the frequency domain current:
% Z = P/U

% The frequency domain can be found by FFT of the time domain signal.
% To get better resolution, use a longer point FFT.

npoints = 2^12;

P_freq = fft(P_save, npoints);
U_freq = fft(U_save, npoints);
Z = P_freq ./ U_freq;

subplot(2,1,1);
plot(abs(Z));
ylim([0, 1000]);  % actual impedances go to infinity, truncate view
xlim([0, npoints/2]);  % Only the first half of the FFT is necessary
title('Input impedance from simulation results')
ylabel('|Z|');


subplot(2,1,2);
% Theory for open circuit input impedance, where Z_load = 0.

f = linspace(0, Fs, npoints);  % frequency vector
B = 2*pi*f / c;   % beta parameter
Z = sqrt(-1) * z0 * tan(B*L);

plot(abs(Z))
ylim([0, 1000])
xlim([0, npoints/2])
title('Input impedance according to theory')
xlabel('frequency')
ylabel('|Z|')

% these plots should match






