kmodel

 kmodel.m --- Terminated K-delay line demo.

 This function demonstrates the terminated K-delay line consisting of N
 K-nodes, connected by K-pipes of equal admittance Y2, as shown in Figure 3
 (the third port of K-node N1 is omitted). The excitation signal Uext is
 connected to N1 with appropriate filtering. Y1 and YN are termination
 admittances. The junction pressure PJ is returned. Except in nodes N1 and NN,
 the admittance matched chain is calculated as in Figure 7. Optional <fIDX>
 flag shows an animation in figure(fIDX). 

 Excitation Examples: 
 Uext = [1 zeros(1,399)]; % Impulse
 Uext = [hamming(10)' zeros(1,390)]; % A smooth excitation signal

 Demo Examples:
 PK = kmodel(Uext,.10,1,50,.10); % Non-inverting terminations.
 PK = kmodel(Uext,5,1,50,10);    % Inverting terminations.
 PK = kmodel(Uext,.10,1,50,1);   % Admittance-matched termination (no
                                   reflection from the right).

0001 function PJ = kmodel(Uext,Y1,Y2,N,YN,fIDX);
0002 % kmodel.m --- Terminated K-delay line demo.
0003 %
0004 % This function demonstrates the terminated K-delay line consisting of N
0005 % K-nodes, connected by K-pipes of equal admittance Y2, as shown in Figure 3
0006 % (the third port of K-node N1 is omitted). The excitation signal Uext is
0007 % connected to N1 with appropriate filtering. Y1 and YN are termination
0008 % admittances. The junction pressure PJ is returned. Except in nodes N1 and NN,
0009 % the admittance matched chain is calculated as in Figure 7. Optional <fIDX>
0010 % flag shows an animation in figure(fIDX).
0011 %
0012 % Excitation Examples:
0013 % Uext = [1 zeros(1,399)]; % Impulse
0014 % Uext = [hamming(10)' zeros(1,390)]; % A smooth excitation signal
0015 %
0016 % Demo Examples:
0017 % PK = kmodel(Uext,.10,1,50,.10); % Non-inverting terminations.
0018 % PK = kmodel(Uext,5,1,50,10);    % Inverting terminations.
0019 % PK = kmodel(Uext,.10,1,50,1);   % Admittance-matched termination (no
0020 %                                   reflection from the right).
0021 
0022 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -*- Mode: Matlab -*- %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0023 %% Author          : Cumhur.Erkut@erase.hut.fi
0024 %% Created On      : Wed Jun 30 12:33:22 2004
0025 %% Last Modified By: Cumhur.Erkut@erase.hut.fi
0026 %% Last Modified On: Wed Jun 30 13:13:05 2004
0027 %% Update Count    : 43
0028 %% Reference       : Digital Waveguides versus Finite Difference Structures:
0029 %%                   Equivalence and Mixed Modeling, Matti Karjalainen and
0030 %%                   Cumhur Erkut. EURASIP J. of Applied Signal Processing,
0031 %%                   Volume 2004, Number 7, pp. 978-989, 15 June 2004.
0032 %%                   http://asp.hindawi.com/volume-2004/S1110865704401176.html
0033 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0034 % ERROR CHECK
0035   error(nargchk(5,6,nargin));
0036   if nargin == 5
0037       fIDX = 0;
0038   end
0039 %%%%%%%%%%%%%%%%%%%% SIMULATION PARAMETERS %%%%%%%%%%%%%%%%%%%%
0040  MAXSTEP = length(Uext);
0041  Umax = max(Uext);
0042 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0043 
0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 % -----------------  INITIALIZATION  --------------------------
0046  % Excitation filtering (Figure 3)
0047    Uext = filter([1 0 -1],1, Uext);
0048 
0049  % INITIALLY RELAXED STATES
0050    P = zeros(3,N);
0051 % -------------------------------------------------------------
0052 
0053 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
0054 % -----------------    RECURSION     --------------------------
0055   for n = 1:MAXSTEP 
0056     
0057     % CALCULATE THE FDTD
0058     P(1,1) = 1/(Y1+Y2)*( 2* Y1 * P(3,1) + Uext(n) + 2* Y2 * P(2,2)) - P(3,1);
0059     P(1,2:N-1) = P(2,1:N-2) + P(2,3:N) -  P(3,2:N-1);
0060     P(1,N) = 1/(YN+Y2)*( 2* YN * P(3,N) + 2* Y2 * P(2,N-1))  - P(3,N);
0061     PJ = P(1,:);
0062 
0063     % % -------------    PLOT    -------------------
0064     if fIDX
0065       figure(fIDX);clf;
0066       set(gcf,'Renderer','OpenGL');
0067 
0068       % Junction Pressure
0069       h = plot(1:N,PJ);set(h,'LineWidth',2);
0070       grid on;
0071       set(gca,'XLim',[1 N],...
0072               'YLim',[-1.1*Umax 1.1*Umax],...
0073               'FontSize',14,...
0074               'FontName','TimesNewRoman');
0075       h = xlabel('Position k'); 
0076           ylabel('Pressure P_{J,k}');
0077           set(h,'VerticalAlignment','middle')
0078           title('K-model')
0079       drawnow;
0080     end
0081 
0082     % UPDATE THE STATES
0083     P(3,:) = P(2,:);
0084     P(2,:) = P(1,:);
0085   end
0086