|
tp@0
|
1 function [ir,ninit] = EDB1wedgeN(cylS,cylR,cylE,ncylrows,nyveclist,edgelengthlist,dzvec,method,...
|
|
tp@0
|
2 pathalongplane,nedgeelcombs,R_irstart,bc,cair,fs,BIGEDGE1stvalue)
|
|
tp@0
|
3 % EDB1wedgeN - Gives the Nth-order diffraction impulse response.
|
|
tp@0
|
4 % ir The Nth order diffraction impulse response for a combination of edges
|
|
tp@0
|
5 % ninit The initial delay that has been cut away before the
|
|
tp@0
|
6 % beginning of ir. Note that the time zero of the ir (including ninit)
|
|
tp@0
|
7 % corresponds to R_irstart.
|
|
tp@0
|
8 %
|
|
tp@0
|
9 % with input parameters
|
|
tp@0
|
10 % cylS List containing [rS,thetaS,zS] = cyl. coord. of the source rel. to the first edge
|
|
tp@0
|
11 % cylR List containing [rR,thetaR,zR] = cyl. coord. of the receiver rel. to the last edge
|
|
tp@0
|
12 % cylE Matrix containing consecutive blocks of
|
|
tp@0
|
13 % [rstartE2_r1,thetastartE2_r1,zstartE2_r1;rendE2_r1,thetaendE2_r1,zendE2_r1]
|
|
tp@0
|
14 % [rstartE1_r2,thetastartE1_r2,zstartE1_r2;rendE1_r2,thetaendE1_r2,zendE1_r2]
|
|
tp@0
|
15 % etc (E3_r2,E2_r3)
|
|
tp@0
|
16 % ncylrows The number of rows in cylE
|
|
tp@0
|
17 % nyveclist List of the wedge indices of the N edges
|
|
tp@0
|
18 % edgelengthlist List of the lengths of the N edges
|
|
tp@0
|
19 % dzvec List of the N edge element sizes
|
|
tp@0
|
20 % method 'n' for the New method or 'v' for Vanderkooys method
|
|
tp@0
|
21 % pathalongplane List of 1 or 0, depending on if the (N-1) paths run along a plane or not
|
|
tp@0
|
22 % nedgeelcombs The value prod(ndiv)
|
|
tp@0
|
23 % R_irstart A distance to which the ir time zero will correspond
|
|
tp@0
|
24 % bc Matrix containing [refl11 refl12 refl21 refl22] where the values are +1 or -1
|
|
tp@0
|
25 % indicating the reflection coefficients of the four wedge planes involved.
|
|
tp@0
|
26 % +1 means rigid and -1 means soft (pressure-release)
|
|
tp@0
|
27 % refl11 is edge 1, plane 1 (refplane)
|
|
tp@0
|
28 % refl12 is edge 1, plane 2
|
|
tp@0
|
29 % refl21 is edge 2, plane 1 (refplane)
|
|
tp@0
|
30 % refl22 is edge 2, plane 2
|
|
tp@0
|
31 % default is all-rigid
|
|
tp@0
|
32 % cair The speed of sound
|
|
tp@0
|
33 % fs The sampling frequency
|
|
tp@0
|
34 % BIGEDGE1stvalue The constant value that should be in the first column
|
|
tp@0
|
35 % of BIGEDGESTEPMATRIX
|
|
tp@0
|
36 %
|
|
tp@0
|
37 % BIGEDGESTEPMATRIX (global)
|
|
tp@0
|
38 % a matrix, size [Nelem,N], where Nelem = nedge1*nedge2*nedge3*....
|
|
tp@0
|
39 % nedge1 is the number of edge
|
|
tp@0
|
40 % elements that edge1 is sub-divided into and nedge2 is the
|
|
tp@0
|
41 % number of edge elements that edge2 is sub-divided into, etc. The
|
|
tp@0
|
42 % matrix contains values between 0 and 1, representing all
|
|
tp@0
|
43 % combinations of relative positions along the two edges for the
|
|
tp@0
|
44 % mid-points of these edge elements. They are organized like
|
|
tp@0
|
45 % below, if nedge1 = 10 and nedge2 = 10 and nedge3 = 25:
|
|
tp@0
|
46 % BIGEDGESTEPMATRIX =
|
|
tp@0
|
47 % [ 0.05 0.05 0.02;
|
|
tp@0
|
48 % 0.05 0.05 0.06;
|
|
tp@0
|
49 % 0.05 0.05 0.10;
|
|
tp@0
|
50 % ...
|
|
tp@0
|
51 % 0.05 0.05 0.98;
|
|
tp@0
|
52 % 0.05 0.15 0.02;
|
|
tp@0
|
53 % ...
|
|
tp@0
|
54 % 0.95 0.95 0.98];
|
|
tp@0
|
55 % These values are multiplied with the lengths of the respective
|
|
tp@0
|
56 % edges to give the actual position, in meters along each edge.
|
|
tp@0
|
57 % New version: The first column isn't included in the matrix, but
|
|
tp@0
|
58 % transferred separately as BIGEDGE1stvalue (see above).
|
|
tp@0
|
59 %
|
|
tp@0
|
60 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
61 % This file is part of the Edge Diffraction Toolbox by Peter Svensson.
|
|
tp@0
|
62 %
|
|
tp@0
|
63 % The Edge Diffraction Toolbox is free software: you can redistribute it and/or modify
|
|
tp@0
|
64 % it under the terms of the GNU General Public License as published by the Free Software
|
|
tp@0
|
65 % Foundation, either version 3 of the License, or (at your option) any later version.
|
|
tp@0
|
66 %
|
|
tp@0
|
67 % The Edge Diffraction Toolbox is distributed in the hope that it will be useful,
|
|
tp@0
|
68 % but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
tp@0
|
69 % FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
|
|
tp@0
|
70 %
|
|
tp@0
|
71 % You should have received a copy of the GNU General Public License along with the
|
|
tp@0
|
72 % Edge Diffraction Toolbox. If not, see <http://www.gnu.org/licenses/>.
|
|
tp@0
|
73 % ----------------------------------------------------------------------------------------------
|
|
tp@0
|
74 % Peter Svensson (svensson@iet.ntnu.no) 20061101
|
|
tp@0
|
75 %
|
|
tp@0
|
76 % [ir,ninit] = EDB1wedgeN(cylS,cylR,cylE,ncylrows,nyveclist,edgelengthlist,dzvec,method,...
|
|
tp@0
|
77 % pathalongplane,nedgeelcombs,R_irstart,bc,cair,fs,BIGEDGE1stvalue);
|
|
tp@0
|
78
|
|
tp@0
|
79 global BIGEDGESTEPMATRIX
|
|
tp@0
|
80
|
|
tp@0
|
81 Ndifforder = length(edgelengthlist);
|
|
tp@0
|
82
|
|
tp@0
|
83 multfac = prod(1./(pathalongplane + 1));
|
|
tp@0
|
84
|
|
tp@0
|
85 %-----------------------------------------------------------------------
|
|
tp@0
|
86 % Create matrices with the z-, r- and theta-values
|
|
tp@0
|
87 %
|
|
tp@0
|
88 % BigB is a matrix that contains relative steps, ]0,1[, along each edge
|
|
tp@0
|
89 %
|
|
tp@0
|
90 % All the big matrices will have one column for each "edge n rel. to edge
|
|
tp@0
|
91 % m" combination:
|
|
tp@0
|
92 % Column 1 E2_RE1
|
|
tp@0
|
93 % Column 2 E1_RE2
|
|
tp@0
|
94 % Column 3 E3_RE2
|
|
tp@0
|
95 % Column 4 E2_RE3
|
|
tp@0
|
96 % Column 5 E4_RE3
|
|
tp@0
|
97 % Column 6 E3_RE4
|
|
tp@0
|
98 % Column 7 E5_RE4
|
|
tp@0
|
99 % Column 8 E4_RE5
|
|
tp@0
|
100 % etc
|
|
tp@0
|
101 % For Ndifforder = 2, there are 2 columns
|
|
tp@0
|
102 % For Ndifforder = 3, there are 4 columns
|
|
tp@0
|
103 % For Ndifforder = 4, there are 6 columns
|
|
tp@0
|
104 % etc
|
|
tp@0
|
105
|
|
tp@0
|
106 onesvec = uint8(ones(nedgeelcombs,1));
|
|
tp@0
|
107
|
|
tp@0
|
108 colvec = [ [2:Ndifforder].' [1:Ndifforder-1].'];
|
|
tp@0
|
109 colvec = reshape(colvec.',1,(Ndifforder-1)*2);
|
|
tp@0
|
110 ivecwithhole = [1 [3:length(colvec)]];
|
|
tp@0
|
111 colvec = colvec(ivecwithhole)-1;
|
|
tp@0
|
112
|
|
tp@0
|
113 %-----------------------------------------------------------------------
|
|
tp@0
|
114 % The z-values for each edge, relative to each own edge can be calculated
|
|
tp@0
|
115 % from BIGEDGESTEPMATRIX.
|
|
tp@0
|
116
|
|
tp@0
|
117 zEn_REn = BIGEDGESTEPMATRIX.*edgelengthlist(onesvec,2:end);
|
|
tp@0
|
118 zEn_REn1stvalue = BIGEDGE1stvalue.*edgelengthlist(1);
|
|
tp@0
|
119
|
|
tp@0
|
120 % The z-values for each edge, relative to the two neighbour edges
|
|
tp@0
|
121
|
|
tp@0
|
122 zreSTA = cylE(1:2:ncylrows,3).';
|
|
tp@0
|
123 zreDIF = cylE(2:2:ncylrows,3).' - zreSTA;
|
|
tp@0
|
124
|
|
tp@0
|
125 zEn_REm = zreSTA(onesvec,ivecwithhole) + BIGEDGESTEPMATRIX(:,colvec).*zreDIF(onesvec,ivecwithhole);
|
|
tp@0
|
126 zEn_REm2ndvalue = zreSTA(2) + BIGEDGE1stvalue.*zreDIF(2);
|
|
tp@0
|
127
|
|
tp@0
|
128 %-----------------------------------------------------------------------
|
|
tp@0
|
129 % The S2E and E2R distances
|
|
tp@0
|
130
|
|
tp@0
|
131 % S2E
|
|
tp@0
|
132 S2Edist = sqrt( cylS(1)^2 + ( zEn_REn1stvalue - cylS(3) ).^2);
|
|
tp@0
|
133
|
|
tp@0
|
134 % E2R
|
|
tp@0
|
135 E2Rdist = sqrt( cylR(1)^2 + ( zEn_REn(:,Ndifforder-1) - cylR(3) ).^2);
|
|
tp@0
|
136
|
|
tp@0
|
137 %-----------------------------------------------------------------------
|
|
tp@0
|
138 % The r-values for each edge, relative to the two neighbour edges
|
|
tp@0
|
139
|
|
tp@0
|
140 rreSTA = cylE(1:2:ncylrows,1).';
|
|
tp@0
|
141 rreDIF = cylE(2:2:ncylrows,1).' - rreSTA;
|
|
tp@0
|
142
|
|
tp@0
|
143 rEn_REm = rreSTA(onesvec,ivecwithhole) + BIGEDGESTEPMATRIX(:,colvec).*rreDIF(onesvec,ivecwithhole);
|
|
tp@0
|
144 rEn_REm2ndvalue = rreSTA(2) + BIGEDGE1stvalue.*rreDIF(2);
|
|
tp@0
|
145
|
|
tp@0
|
146 %-----------------------------------------------------------------------
|
|
tp@0
|
147 % The edge-to-edge distances
|
|
tp@0
|
148
|
|
tp@0
|
149 % E2Edist will be "edge2 re 1","edge3 re 2" etc
|
|
tp@0
|
150 modcols = [1 [3:2:Ndifforder*2-3]-1];
|
|
tp@0
|
151 E2Edist = sqrt( ([zEn_REn1stvalue*double(onesvec) zEn_REn(:,1:Ndifforder-2)] - zEn_REm(:,modcols)).^2 + rEn_REm(:,modcols).^2 );
|
|
tp@0
|
152
|
|
tp@0
|
153 %-----------------------------------------------------------------------
|
|
tp@0
|
154 % Calculate the coshnyeta values first, to use them in the directivity
|
|
tp@0
|
155 % functions
|
|
tp@0
|
156
|
|
tp@0
|
157 if method(1) == 'n'
|
|
tp@0
|
158
|
|
tp@0
|
159 % CH will get the coshnyeta values for the directivity functions
|
|
tp@0
|
160
|
|
tp@0
|
161 CH = zeros(nedgeelcombs,Ndifforder-1);
|
|
tp@0
|
162
|
|
tp@0
|
163 %-----------------------------------------------------------------------
|
|
tp@0
|
164 % For the coshnyeta values we need the quantity
|
|
tp@0
|
165 %
|
|
tp@0
|
166 % ( ( ze1_RE1 -zs_RE1 )*( ze1_RE1 - ze2_RE1 ) + n*m )/rs/re2_RE1
|
|
tp@0
|
167
|
|
tp@0
|
168 CH(:,1) = (( zEn_REn1stvalue - cylS(3) ).*( zEn_REn1stvalue - zEn_REm(:,1) ) + S2Edist.*E2Edist(:,1))/cylS(1)./rEn_REm(:,1);
|
|
tp@0
|
169 CH(:,1) = ( real(sqrt( CH(:,1).^2-1)) + CH(:,1) ).^nyveclist(1);
|
|
tp@0
|
170 CH(:,1) = ( CH(:,1) + 1./CH(:,1))/2;
|
|
tp@0
|
171
|
|
tp@0
|
172 for ii = 2:Ndifforder-1
|
|
tp@0
|
173 mcol1 = ii*2-3;
|
|
tp@0
|
174 mcol2 = ii*2-2;
|
|
tp@0
|
175 if ii == 2
|
|
tp@0
|
176 CH(:,ii) = ( ( zEn_REn(:,ii-1) - zEn_REm2ndvalue ).*( zEn_REn(:,ii-1) - zEn_REm(:,mcol2) ) + E2Edist(:,ii-1).*E2Edist(:,ii) )./rEn_REm2ndvalue./rEn_REm(:,mcol2);
|
|
tp@0
|
177 else
|
|
tp@0
|
178 CH(:,ii) = ( ( zEn_REn(:,ii-1) - zEn_REm(:,mcol1) ).*( zEn_REn(:,ii-1) - zEn_REm(:,mcol2) ) + E2Edist(:,ii-1).*E2Edist(:,ii) )./rEn_REm(:,mcol1)./rEn_REm(:,mcol2);
|
|
tp@0
|
179
|
|
tp@0
|
180 end
|
|
tp@0
|
181 CH(:,ii) = ( sqrt( CH(:,ii).^2-1) + CH(:,ii) ).^nyveclist(ii);
|
|
tp@0
|
182 CH(:,ii) = real( CH(:,ii) + 1./CH(:,ii))/2;
|
|
tp@0
|
183 end
|
|
tp@0
|
184
|
|
tp@0
|
185 CHN = ( ( zEn_REn(:,Ndifforder-1) - zEn_REm(:,2*Ndifforder-3) ).*( zEn_REn(:,Ndifforder-1) - cylR(3) ) + E2Edist(:,Ndifforder-1).*E2Rdist )./rEn_REm(:,2*Ndifforder-3)./cylR(1);
|
|
tp@0
|
186 CHN = ( sqrt( CHN.^2-1) + CHN ).^nyveclist(Ndifforder);
|
|
tp@0
|
187 CHN = real( CHN + 1./CHN)/2;
|
|
tp@0
|
188
|
|
tp@0
|
189 elseif method(1) == 'v'
|
|
tp@0
|
190 CH = ones(1,Ndifforder-1);
|
|
tp@0
|
191 CHN = 1;
|
|
tp@0
|
192 else
|
|
tp@0
|
193 error('Method not allowed')
|
|
tp@0
|
194 end
|
|
tp@0
|
195
|
|
tp@0
|
196 %-----------------------------------------------------------------------
|
|
tp@0
|
197 % The theta-values for each edge, relative to the two neighbour edges
|
|
tp@0
|
198
|
|
tp@0
|
199 thetareSTA = cylE(1:2:ncylrows,2).';
|
|
tp@0
|
200 thetareDIF = cylE(2:2:ncylrows,2).' - thetareSTA;
|
|
tp@0
|
201
|
|
tp@0
|
202 if all(thetareDIF==0)
|
|
tp@0
|
203 thetaEn_REm = thetareSTA(ivecwithhole);
|
|
tp@0
|
204 else
|
|
tp@0
|
205 % CHECK: Here we need to implement something like
|
|
tp@0
|
206 thetaEn_REm = acos(rEn_REm./E2Edist);
|
|
tp@0
|
207 disp(['ERROR: Not implemented the proper calculation of theta-values for non-parallel edge pairs yet!'])
|
|
tp@0
|
208 % old version:
|
|
tp@0
|
209 % thetaEn_REm = thetareSTA(onesvec,ivecwithhole) + BIGEDGESTEPMATRIX(:,colvec).*thetareDIF(onesvec,ivecwithhole);
|
|
tp@0
|
210 end
|
|
tp@0
|
211
|
|
tp@0
|
212 thetaEn_REm2ndvalue = thetareSTA(2) + BIGEDGE1stvalue.*thetareDIF(2);
|
|
tp@0
|
213
|
|
tp@0
|
214 %------------------------------------------------------
|
|
tp@0
|
215 % Calculate the product of all directivity factors and inverse distances
|
|
tp@0
|
216
|
|
tp@0
|
217 % First "leg" from source, via the first edge, to the second edge, but
|
|
tp@0
|
218 % including the DF only for the first edge
|
|
tp@0
|
219
|
|
tp@0
|
220 if pathalongplane(1) == 0
|
|
tp@0
|
221 AMP = multfac*prod(nyveclist)*prod(dzvec)*(-1/4/pi)^Ndifforder*(sin(nyveclist(1)*(pi + cylS(2) + thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi + cylS(2) + thetaEn_REm(:,1) ))) + ...
|
|
tp@0
|
222 sin(nyveclist(1)*(pi + cylS(2) - thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi + cylS(2) - thetaEn_REm(:,1) ))) + ...
|
|
tp@0
|
223 sin(nyveclist(1)*(pi - cylS(2) + thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi - cylS(2) + thetaEn_REm(:,1) ))) + ...
|
|
tp@0
|
224 sin(nyveclist(1)*(pi - cylS(2) - thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi - cylS(2) - thetaEn_REm(:,1) ))))./S2Edist./E2Edist(:,1);
|
|
tp@0
|
225 else
|
|
tp@0
|
226 if cylS(2) == 0
|
|
tp@0
|
227 AMP = 4*multfac*prod(nyveclist)*prod(dzvec)*(-1/4/pi)^Ndifforder*(sin(nyveclist(1)*(pi + thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi + thetaEn_REm(:,1) ))))./S2Edist./E2Edist(:,1);
|
|
tp@0
|
228 else
|
|
tp@0
|
229 AMP = 2*multfac*prod(nyveclist)*prod(dzvec)*(-1/4/pi)^Ndifforder*(sin(nyveclist(1)*(pi + cylS(2) + thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi + cylS(2) + thetaEn_REm(:,1) ))) + ...
|
|
tp@0
|
230 sin(nyveclist(1)*(pi - cylS(2) + thetaEn_REm(:,1) ))./(CH(:,1)- cos(nyveclist(1)*(pi - cylS(2) + thetaEn_REm(:,1) ))))./S2Edist./E2Edist(:,1);
|
|
tp@0
|
231 end
|
|
tp@0
|
232 end
|
|
tp@0
|
233
|
|
tp@0
|
234 % Second "leg" form edge 1, via the second edge, to the third edge
|
|
tp@0
|
235 % including the DF only for the second edge. We keep this out of the
|
|
tp@0
|
236 % for-loop because the matrix thetaEn_REm2ndvalue gives a special
|
|
tp@0
|
237 % formulation.
|
|
tp@0
|
238
|
|
tp@0
|
239 for ii = 2:2
|
|
tp@0
|
240 iv2 = 2;
|
|
tp@0
|
241 if pathalongplane(ii-1)*pathalongplane(ii) == 0
|
|
tp@0
|
242 AMP = AMP.*(sin(nyveclist(ii)*(pi + thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
243 sin(nyveclist(ii)*(pi + thetaEn_REm2ndvalue - thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm2ndvalue - thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
244 sin(nyveclist(ii)*(pi - thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi - thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
245 sin(nyveclist(ii)*(pi - thetaEn_REm2ndvalue - thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi - thetaEn_REm2ndvalue - thetaEn_REm(:,iv2)))))./E2Edist(:,2);
|
|
tp@0
|
246 else
|
|
tp@0
|
247 AMP = 4*AMP.*(sin(nyveclist(ii)*(pi + thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm2ndvalue + thetaEn_REm(:,iv2)))))./E2Edist(:,2);
|
|
tp@0
|
248 end
|
|
tp@0
|
249 end
|
|
tp@0
|
250
|
|
tp@0
|
251 % All the following "legs", except the last one that reaches the receiver.
|
|
tp@0
|
252
|
|
tp@0
|
253 for ii = 3:Ndifforder-1
|
|
tp@0
|
254 iv1 = ii*2-3;
|
|
tp@0
|
255 iv2 = iv1+1;
|
|
tp@0
|
256 if pathalongplane(ii-1)*pathalongplane(ii) == 0
|
|
tp@0
|
257 AMP = AMP.*(sin(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
258 sin(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) - thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) - thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
259 sin(nyveclist(ii)*(pi - thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi - thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))) + ...
|
|
tp@0
|
260 sin(nyveclist(ii)*(pi - thetaEn_REm(:,iv1) - thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi - thetaEn_REm(:,iv1) - thetaEn_REm(:,iv2)))))./E2Edist(:,ii);
|
|
tp@0
|
261 else
|
|
tp@0
|
262 AMP = 4*AMP.*(sin(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))./(CH(:,ii)-cos(nyveclist(ii)*(pi + thetaEn_REm(:,iv1) + thetaEn_REm(:,iv2)))))./E2Edist(:,ii);
|
|
tp@0
|
263 end
|
|
tp@0
|
264 end
|
|
tp@0
|
265
|
|
tp@0
|
266 % The last "leg" that reaches the receiver, and including the last
|
|
tp@0
|
267 % directivity factor.
|
|
tp@0
|
268
|
|
tp@0
|
269 if pathalongplane(Ndifforder-1) == 0
|
|
tp@0
|
270 AMP = AMP.*(sin(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))) + ...
|
|
tp@0
|
271 sin(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))) + ...
|
|
tp@0
|
272 sin(nyveclist(Ndifforder)*(pi - thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi - thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))) + ...
|
|
tp@0
|
273 sin(nyveclist(Ndifforder)*(pi - thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi - thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))))./E2Rdist;
|
|
tp@0
|
274 else
|
|
tp@0
|
275 AMP = 2*AMP.*(sin(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) + cylR(2)))) + ...
|
|
tp@0
|
276 sin(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))./(CHN-cos(nyveclist(Ndifforder)*(pi + thetaEn_REm(:,2*Ndifforder-3) - cylR(2)))))./E2Rdist;
|
|
tp@0
|
277 end
|
|
tp@0
|
278
|
|
tp@0
|
279 %-------------------------------------------------------------------------
|
|
tp@0
|
280 % Determine in which sample slots the amplitude contributions should
|
|
tp@0
|
281 % be added, based on the total distance.
|
|
tp@0
|
282
|
|
tp@0
|
283 % B2 is number of sample slots as non-integers
|
|
tp@0
|
284
|
|
tp@0
|
285 if Ndifforder == 2
|
|
tp@0
|
286 B2 = ( S2Edist + E2Edist + E2Rdist - R_irstart)/(cair/fs)+1;
|
|
tp@0
|
287 else
|
|
tp@0
|
288 B2 = ( S2Edist + sum(E2Edist.').' + E2Rdist - R_irstart)/(cair/fs)+1;
|
|
tp@0
|
289 end
|
|
tp@0
|
290
|
|
tp@0
|
291 % B1 is the sampleslot numbers in integer number
|
|
tp@0
|
292
|
|
tp@0
|
293 B1 = floor(B2);
|
|
tp@0
|
294
|
|
tp@0
|
295 % B2 is a value between 0 and 1 stating how much of dh that should be added
|
|
tp@0
|
296 % to the first sample slot, i.e. the one given on B1. The rest of dh should
|
|
tp@0
|
297 % be added to the following sample slot.
|
|
tp@0
|
298
|
|
tp@0
|
299 B2 = B2 - B1;
|
|
tp@0
|
300
|
|
tp@0
|
301 % New approach: simply state that the amplitudevalue, mult. by 1-B2, should
|
|
tp@0
|
302 % be placed in the integer slot given in B1. In addition, the
|
|
tp@0
|
303 % amplitudevalues, mult by B2, should be placed in the integer slot given
|
|
tp@0
|
304 % by B1+1.
|
|
tp@0
|
305
|
|
tp@0
|
306 B1 = [B1;B1+1];
|
|
tp@0
|
307
|
|
tp@0
|
308 AMP = [AMP.*(1-B2);AMP.*B2];
|
|
tp@0
|
309
|
|
tp@0
|
310 ir = zeros(max(B1),1 );
|
|
tp@0
|
311
|
|
tp@0
|
312 [B1,sortvec] = sort(B1);
|
|
tp@0
|
313 AMP = cumsum(AMP(sortvec));
|
|
tp@0
|
314 ivstep = (find(diff([(B1);(B1(end))+1])));
|
|
tp@0
|
315 ir(B1(ivstep)) = diff([0;AMP(ivstep)]);
|
|
tp@0
|
316
|
|
tp@0
|
317 % Here we could have the possibility to save space by cutting out initial zeros and
|
|
tp@0
|
318 % set the value ninit correspondingly to the number of removed zeros. However, the rest of
|
|
tp@0
|
319 % the EDB1 functions have not implented support for this.
|
|
tp@0
|
320
|
|
tp@0
|
321 ninit = 0;
|
