hybridmodel.m --- Mixed Modeling demo.
This function demonstrates a hybrid waveguide model consisting of N/2 K-nodes,
N/2 W-nodes, and the KW converter, as shown in Figure 5. The excitation signal
Uext is connected to the K-node N1 with appropriate filtering. Y1 and YN are
termination admittances, and Y2 is the admittance of both the K-pipes and
W-lines, The junction pressure PJ is returned. 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:
PH = hybridmodel(Uext,.10,1,50,.10); % Non-inverting terminations.
PH = hybridmodel(Uext,5,1,50,10); % Inverting terminations.
PH = hybridmodel(Uext,.10,1,50,1); % Admittance-matched termination (no
reflection from the right).0001 function PJ = hybridmodel(Uext,Y1,Y2,N,YN,fIDX); 0002 % hybridmodel.m --- Mixed Modeling demo. 0003 % 0004 % This function demonstrates a hybrid waveguide model consisting of N/2 K-nodes, 0005 % N/2 W-nodes, and the KW converter, as shown in Figure 5. The excitation signal 0006 % Uext is connected to the K-node N1 with appropriate filtering. Y1 and YN are 0007 % termination admittances, and Y2 is the admittance of both the K-pipes and 0008 % W-lines, The junction pressure PJ is returned. Optional <fIDX> flag shows an 0009 % animation in figure(fIDX). 0010 % 0011 % Excitation Examples: 0012 % Uext = [1 zeros(1,399)]; % Impulse 0013 % Uext = [hamming(10)' zeros(1,390)]; % A smooth excitation signal 0014 % 0015 % Demo Examples: 0016 % PH = hybridmodel(Uext,.10,1,50,.10); % Non-inverting terminations. 0017 % PH = hybridmodel(Uext,5,1,50,10); % Inverting terminations. 0018 % PH = hybridmodel(Uext,.10,1,50,1); % Admittance-matched termination (no 0019 % reflection from the right). 0020 0021 %%%%%%%%%%%%%%%%%%%%%%%%%%%%% -*- Mode: Matlab -*- %%%%%%%%%%%%%%%%%%%%%%%%%%%% 0022 %% Author : Cumhur.Erkut@erase.hut.fi 0023 %% Created On : Wed Jun 30 12:33:22 2004 0024 %% Last Modified By: Cumhur.Erkut@erase.hut.fi 0025 %% Last Modified On: Thu Jul 01 12:10:02 2004 0026 %% Update Count : 208 0027 %% Reference : Digital Waveguides versus Finite Difference Structures: 0028 %% Equivalence and Mixed Modeling, Matti Karjalainen and 0029 %% Cumhur Erkut. EURASIP J. of Applied Signal Processing, 0030 %% Volume 2004, Number 7, pp. 978-989, 15 June 2004. 0031 %% http://asp.hindawi.com/volume-2004/S1110865704401176.html 0032 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0033 % ERROR CHECK 0034 error(nargchk(5,6,nargin)); 0035 if nargin == 5 0036 fIDX = 0; 0037 end 0038 %%%%%%%%%%%%%%%%%%%% SIMULATION PARAMETERS %%%%%%%%%%%%%%%%%%%% 0039 MAXSTEP = length(Uext); 0040 NK = fix(.5*N); 0041 NW = N - NK; 0042 Umax = max(Uext); 0043 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0044 0045 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0046 % ----------------- INITIALIZATION -------------------------- 0047 % Excitation filtering (Figure 3) 0048 Uext = filter([1 0 -1],1, Uext); 0049 0050 % INITIALLY RELAXED STATES 0051 P = zeros(3,NK); 0052 PJ = zeros(1,N); 0053 0054 % DELAY LINES 0055 P2 = 0; PN = 0; P1m = 0; PNp = 0; 0056 Pm = zeros(1,NW-1); 0057 Pp = zeros(1,NW-1); 0058 0059 % KW-Converter memory 0060 kmem0 = 0; 0061 wmem0 = 0; 0062 wmem1 = 0; 0063 % ------------------------------------------------------------- 0064 0065 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0066 % ----------------- RECURSION -------------------------- 0067 for n = 1:MAXSTEP 0068 0069 % CALCULATE THE K-structure 0070 P(1,1) = 1/(Y1+Y2)*( 2* Y1 * P(3,1) + Uext(n) + 2* Y2 * P(2,2)) - P(3,1); 0071 P(1,2:NK-1) = P(2,1:NK-2) + P(2,3:NK) - P(3,2:NK-1); 0072 P(1,NK) = P(2,NK-1) + kmem0 - P(3,NK); 0073 0074 % CALCULATE THE W-nodes 0075 P2 = (P(2,NK) - wmem1 + Pp(1)); 0076 PN = 1/(Y2+YN) * (2* Y2 * Pm(NW-1)); 0077 0078 % CALCULATE THE SCATTERED WAVES 0079 P1m = P2 - Pp(1); 0080 PNp = PN - Pm(NW-1); 0081 KWin = Pp(1); 0082 0083 % CALCULATE THE JUNCTION PRESSURES 0084 PJ = [P(1,:) P2 Pm(1:NW-2)+Pp(2:NW-1) PN]; 0085 0086 % UPDATE THE KW-Converter Memory 0087 kmem0 = P2; 0088 wmem1 = wmem0; 0089 wmem0 = KWin; 0090 0091 % UPDATE THE STATES 0092 P(3,:) = P(2,:); 0093 P(2,:) = P(1,:); 0094 0095 % UPDATE DELAY LINES 0096 Pm = [P1m Pm(1:NW-2)]; 0097 Pp = [Pp(2:NW-1) PNp]; 0098 0099 0100 % % ------------- PLOT ------------------- 0101 if fIDX 0102 figure(fIDX);clf; 0103 set(gcf,'Renderer','OpenGL'); 0104 0105 % Junction Pressure 0106 h = plot(1:N,PJ);set(h,'LineWidth',2); 0107 grid on; 0108 set(gca,'XLim',[1 N],... 0109 'YLim',[-1.1*Umax 1.1*Umax],... 0110 'FontSize',14,... 0111 'FontName','TimesNewRoman'); 0112 h = xlabel('Position k'); 0113 ylabel('Pressure P_{J,k}'); 0114 set(h,'VerticalAlignment','middle') 0115 title('Hybrid Model') 0116 drawnow; 0117 hold on; plot([NW NW],.75*[-Umax Umax],'r-.'); hold off; 0118 h=text(0.4*NK,.5*Umax,'K-model'); 0119 h=text(NK+0.4*NW,.5*Umax,'W-model'); 0120 h=text(0.85*NK,.85*Umax,'KW-Converter'); 0121 end 0122 end 0123