0
votes

I'm looking for a way to speed up some simple two port matrix calculations. See the below code example for what I'm doing currently. In essence, I create a [Nx1] frequency vector first. I then loop through the frequency vector and create the [2x2] matrices H1 and H2 (all functions of f). A bit of simple matrix math including a matrix left division '\' later, and I got my result pb as a [Nx1] vector. The problem is the loop - it takes a long time to calculate and I'm looking for way to improve efficiency of the calculations. I tried assembling the problem using [2x2xN] transfer matrices, but the mtimes operation cannot handle 3-D multiplications.

Can anybody please give me an idea how I can approach such a calculation without the need for looping through f?

Many thanks: svenr

% calculate frequency and wave number vector

f = linspace(20,200,400);
w = 2.*pi.*f;

% calculation for each frequency w

for i=1:length(w)

   H1(i,1) = {[1, rho*c*k(i)^2 / (crad*pi); 0,1]};
   H2(i,1) = {[1, 1i.*w(i).*mp; 0, 1]};
   HZin(i,1) = {H1{i,1}*H2{i,1}};
   temp_mat = HZin{i,1}*[1; 0];
   Zin(i,1) = temp_mat(1,1)/temp_mat(2,1);
   temp_mat= H1{i,1}\[1; 1/Zin(i,1)];
   pb(i,1) = temp_mat(1,1); Ub(i,:) = temp_mat(2,1);
end
1
why the curly parenthesis {}?bla
try to use bsxfun with permute to solve the 2x2xn dimension issue.bla
@Nishant: No. 1i is the complex number i.rayryeng
One thing I can suggest. If it's seriously just a 2 x 2 matrix, save overhead by using the 2 x 2 formula for calculating the inverse: invA = (1/(ad - bc))*([d -b; -c a]);, given that A = [a b; c d]; After this, multiply by your 2 x 1 vector to get the answer.rayryeng
Rayryengs option of using the 2x2 matrix formula is a good way to speed up code. Other than that, have you tried the cell function cellfun?patrik

1 Answers

0
votes

Assuming that length(w) == length(k) returns true , rho , c, crad, mp are all scalars and in the last line is Ub(i,1) = temp_mat(2,1) instead of Ub(i,:) = temp_mat(2,1);

 temp = repmat(eyes(2),[1 1 length(w)]);
 temp1(1,2,:) = rho*c*(k.^2)/crad/pi;
 temp2(1,2,:) = (1i.*w)*mp;
 H1 = permute(num2cell(temp1,[1 2]),[3 2 1]);
 H2 = permute(num2cell(temp2,[1 2]),[3 2 1]);
 HZin =  cellfun(@(a,b)(a*b),H1,H2,'UniformOutput',0);
 temp_cell = cellfun(@(a,b)(a*b),H1,repmat({[1;0]},length(w),1),'UniformOutput',0);
 Zin_cell = cellfun(@(a)(a(1,1)/a(2,1)),temp_cell,'UniformOutput',0);
 Zin = cell2mat(Zin);
 temp2_cell = cellfun(@(a)({[1;1/a]}),Zin_cell,'UniformOutput',0);
 temp3_cell = cellfun(@(a,b)(pinv(a)*b),H1,temp2_cell);
 temp4 = cell2mat(temp3_cell);
 p(:,1) = temp4(1:2:end-1);
 Ub(:,1) = temp4(2:2:end);