modal parameter code, run out of memory 1

[mkoijn]

The following code is part of a program I've cobbled together to extract modal parameters from multiple frf's, it extracts the polynomial coeffiients. It's adapted from the diamond modal toolbox and I'm running it in GNU Octave which is an open source matlab clone.

The problem I'm having is that I run out of memeory when I use more than about 40 frf's ( using synthesised data consisting of 5 modes, 190 points per frf), also at 40 frf's it's taking about 2 minutes to calculate the polynomial coefficients ( 1GHz processor,380Mb ram). The function rpcoeffs is pretty much as it came from diamond
Has anyone an idea why this is happening? Does it react similarly in Matlab? any help would be gratefully received

If anyone would like to try it,the input arguments for rpcoeffs are as follows:
h = matrix of frf's arranged in columns
wi= vector of frequencies for frf's
m= order of numerator polynomial ( I used number of modes*2-1+4 for extra terms)
n=order of denominator polynomial (number of modes*2)





function [A,B] = rpcoeffs(h,wi,m,n)

% Usage: [A,B] = rpcoeffs(h,wi,m,n)


% List of variables:
% m = order of numerator polynomial
% n = order of denomenator polynomial
% L = number of frequency bins to curve fit
% h = FRF matrix values at selected frequencies
% wi = vector of frequencies to use (indexed 1 to L)
% Pnum, Pden = matrices of complex orthogonal Forsythe Polynomials (one-sided)
% row i = ith frequency
% column j = (j-1)th order polynomial
% qnum = weights for numerator Forsythe poly. (m+1 x 1)
% qden = weights for denomenator Forsythe poly. (n+1 x 1)
%

%
% Normalize frequency variables
% and make FRF curve 2-sided

j = sqrt(-1);
L1 = length(wi);

wmax = wi(L1);

wi = wi/wmax;
wi = [-wi(L1:-1:1);wi];
xi = j*wi;

h = [conj(h(L1:-1:1,);h];

L = length(wi);
nh = size(h,2); % Total number of FRFs to Curve-fit

%
% Generate Orthogonal Forsythe Polynomials
%

% Numerator Polynomials

qnum = ones(L,1); % numerator poly. weights

[Pnum,unum,vnum] = forsythe(xi,qnum,m);


% Denomenator Polynomials

qden = zeros(L,1);

for i = 1:L,
qden(i) = max(svd(h(i,)); % Use CMIF as denomenator poly. weights
end

[Pden,uden,vden] = forsythe(xi,qden,n);


%
% Form matrices for Least-Squares solution
%

P = Pnum;

SOLROW = zeros(L*nh,(m+1)*nh + n);
WCUM = zeros(L*nh,1);

for hc = 1:nh,

T = zeros(L,n);

for i = 1:n,

T,i) = Pden,i) .* h,hc);

end

W = h,hc) .* Pden,n+1);

SOLROW((hc-1)*L + 1 : hc*L , (hc-1)*(m+1) + 1 : hc*(m+1)) = P;
SOLROW((hc-1)*L + 1 : hc*L , (m+1)*nh + 1 : (m+1)*nh + n) = -T;

WCUM((hc-1)*L + 1 : hc*L,1 ) = W;

end

%
% Solve for the Forsythe polynomial coefficients
%

SOL = SOLROW \ WCUM;

for hc = 1:nh,

C,hc) = SOL((hc-1)*(m+1) + 1: hc*(m+1));

end

D = [SOL((m+1)*nh + 1m+1)*nh + n);1];

%
% Convert Forsythe Poly. coefficients back to normal poly. coeffs.
%

bkm_num = bkm_comp(m,unum,vnum);
bkm_den = bkm_comp(n,uden,vden);

for hc = 1:nh,

% Numerator

for k = 0:m,

C1(k+1) = 0;

for mm = k:m,

C1(k+1) = C1(k+1) + C(mm+1,hc) * bkm_num(k+1,mm+1);

end

C1(k+1) = C1(k+1) / wmax^k;
end

A,hc) = real(C1(length(C1):-1:1))';

end

% Denomenator

for k = 0:n,

D1(k+1) = 0;

for mm = k:n,

D1(k+1) = D1(k+1) + D(mm+1) * bkm_den(k+1,mm+1);

end

D1(k+1) = D1(k+1) / wmax^k;

end

B = real(D1(length(D1):-1:1))';



return

_____________________________________________

function [P,u,v] = forsythe(xi,q,m)

L = length(xi);
P = zeros(L,m+1);
u = zeros(m+1,1);
v = zeros(m+1,1);

%
% Define zero-th order poly.
%

P,1) = ones(L,1);

%
% Define first order poly.
%

d(1) = sum( P,1).^2 .* q.^2);
u(2) = sum( xi .* P,1).^2 .* q.^2) / d(1);

P,2) = (xi - u(2)) .* P,1);

%
% Define higher-order polys.
%

for k = 3m+1),

d(k-1) = sum( P,k-1).^2 .* q.^2);

v(k-1) = sum( xi .* P,k-1) .* P,k-2) .* q.^2) / d(k-2);

u(k) = sum( xi .* P,k-1).^2 .* q.^2) / d(k-1);

P,k) = (xi - u(k)) .* P,k-1) - v(k-1) .* P,k-2);

end

d(m+1) = sum( P,m+1).^2 .* q.^2);

return

________________________________________________________________

function bkm = bkm_comp(mmax,u,v)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k = 0:mmax,

for m = k:mmax

% First, call bkmcheck to see if recursion formula applies

bv = bkmcheck(k,m);

if bv == -1,

bv1 = bkmcheck(k-1,m-1);
bv2 = bkmcheck(k,m-1);
bv3 = bkmcheck(k,m-2);

if bv1 == -1,

bkm1 = bkm(k,m);

else

bkm1 = bv1;

end

if bv2 == -1,

bkm2 = bkm(k+1,m);

else

bkm2 = bv2;

end

if bv3 == -1,

bkm3 = bkm(k+1,m-1);

else

bkm3 = bv3;

end

bkm(k+1,m+1) = bkm1 - u(m+1)*bkm2 - v(m)*bkm3;

else

bkm(k+1,m+1) = bv;

end

end

end

return

__________________________________________________________________

function bv = bkmcheck(k,m)
%
% Name: bkmcheck
%
% Usage: bv = bkmcheck(k,m)
%
% Description: Subroutine of bkm_comp -- checks
% value of k and m to see which
% formula applies
%


if k < 0 | k > m,

bv = 0;

elseif k == m,

bv = 1;

else

bv = -1; % Need to use full formula (or do table lookup)

end


return

 

[MikeyP]

I have run polyreference curve fits on data with >80 frfs at 4096 frequency points on a pc with a similar spec to yours without too much trouble. Try debugging the function line by line with a &quot;whos&quot; after each assignment operation and keep an eye on the memory situation. Maybe the function generates multiple copies of the frfs or something and doesn't free up memory as it goes.

M

 

[mkoijn]

Thanks for the help MikeyP. I followed your suggestion debugging the code using &quot;whos&quot;.
The problem with the memory is the following lines in function rpcoeffs see below

SOLROW is a complex matrix of size 11520 rows * 430 cols and WCUM is size 11520 rows * 1 cols
this eats about 180Mb of memory with 30 frf's and 190 points.

the line: SOL = SOLROW \ WCUM
takes 85 seconds to calculate. The whole function takes 89 seconds to calculate.
SOL is a complex matrix of size 430 rows * 1 cols

my problem now is I don't know what to do next!? I'm out of my depth trying to reformulate
it, if it's possible to reformulate at all. could you give me some clues?

many thanks again

finnigan

__________________________________________________________________
%
% Form matrices for Least-Squares solution
%

P = Pnum;

SOLROW = zeros(L*nh,(m+1)*nh + n);
WCUM = zeros(L*nh,1);

for hc = 1:nh,

T = zeros(L,n);

for i = 1:n,

T,i) = Pden,i) .* h,hc);

end

W = h,hc) .* Pden,n+1);

SOLROW((hc-1)*L + 1 : hc*L , (hc-1)*(m+1) + 1 : hc*(m+1)) = P;
SOLROW((hc-1)*L + 1 : hc*L , (m+1)*nh + 1 : (m+1)*nh + n) = -T;

WCUM((hc-1)*L + 1 : hc*L,1 ) = W;

end

%
% Solve for the Forsythe polynomial coefficients
%

SOL = SOLROW \ WCUM;

 

[MikeyP]

Without rewriting anything, the first thing to do is to reduce the order of the basis function polynomials. What values do you currently have for m and n? Could they be reduced while maintaining accuracy?

Secondly, can you identify the modes in sections? ie. curve fit bins 1-50 then 51-100 etc?

Thirdly, do you need to use the data from all 30 frfs? Could the data from just 10 of the positions be enough to idenify the system?

Finally, you could always buy more memory!

Like I said before, I have used the polyreference approach to curve fit frfs with 9 inputs 9 outputs and 4096 frequency bins and it takes a lot less memory. The Polyreference (called Least-Squares Complex Exponential LSCE for a single input case) fits basis functions of the form A*exp((a+ib)t) to the IRF (the inverse fft of the FRF).

M

 

[maruli]

I think you could make your program more efficient by using half-range polynoms which would yield to half size matrix P, T and W; and also you could solve D and C separatelly. I made the second modification resulting on this source code:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [A,B] = rpcoeffs2(h,wi,m,n)

% Usage: [A,B] = rpcoeffs(h,wi,m,n)


% List of variables:
% m = order of numerator polynomial
% n = order of denomenator polynomial
% L = number of frequency bins to curve fit
% h = FRF matrix values at selected frequencies
% wi = vector of frequencies to use (indexed 1 to L)
% Pnum, Pden = matrices of complex orthogonal Forsythe Polynomials (one-sided)
% row i = ith frequency
% column j = (j-1)th order polynomial
% qnum = weights for numerator Forsythe poly. (m+1 x 1)
% qden = weights for denomenator Forsythe poly. (n+1 x 1)
%

%
% Normalize frequency variables
% and make FRF curve 2-sided

j = sqrt(-1);
L1 = length(wi);

wmax = wi(L1);

wi = wi/wmax;
wi = [-wi(L1:-1:1);wi];
xi = j*wi;

h = [conj(h(L1:-1:1,);h];

L = length(wi);
nh = size(h,2); % Total number of FRFs to Curve-fit

%
% Generate Orthogonal Forsythe Polynomials
%

% Numerator Polynomials

qnum = ones(L,1); % numerator poly. weights

[Pnum,unum,vnum] = forsythe(xi,qnum,m);

for hc=1:nh
% Denomenator Polynomials
qden=h,hc).^2;
[Pden,uden,hc),vden,hc)]=forsythe(xi,qden,n);

%Form matrices for least Squares Solution.
for s=1:n
T,s)=Pden,s).*h,hc);
end
X,:,hc)=-real(Pnum'*T);
Y,:,hc)=real(Pnum'*Pnum);
Z,:,hc)=real(T'*T);

W=h,hc).*Pden,n+1);
G,hc)=real(Pnum'*W);

Vg((hc-1)*+1:hc*,1)=G,hc);
Ug((hc-1)*+1:hc*,=-Y,:,hc)*(X,:,hc)'\Z,:,hc))+X,:,hc);
end

%Solve for the Forsythe Polynomial coeficients.
D=(Ug'*Ug)\(Ug'*Vg);

for hc=1:nh
C,hc)=-(X,:,hc)'\Z,:,hc))*D;
end

D=[D;1];

%
% Convert Forsythe Poly. coefficients back to normal poly. coeffs.
%

bkm_num = bkm_comp(m,unum,vnum);
bkm_den = bkm_comp(n,uden,vden);

for hc = 1:nh

% Numerator

for k = 0:m

C1(k+1) = 0;

for mm = k:m

C1(k+1) = C1(k+1) + C(mm+1,hc) * bkm_num(k+1,mm+1);

end

C1(k+1) = C1(k+1) / wmax^k;
end

A,hc) = real(C1(length(C1):-1:1))'

end

% Denomenator

for k = 0:n

D1(k+1) = 0;

for mm = k:n

D1(k+1) = D1(k+1) + D(mm+1) * bkm_den(k+1,mm+1);

end

D1(k+1) = D1(k+1) / wmax^k;

end

B = real(D1(length(D1):-1:1))'

return





%------------------------------------------------------------------------
function [P,u,v,Tac] = forsythe(xi,q,m)

L = length(xi);
P = zeros(L,m+1);
u = zeros(m+1,1);
v = zeros(m+1,1);

%
% Define zero-th order poly.
%

P,1) = ones(L,1);

%
% Define first order poly.
%

d(1) = sum( P,1).^2 .* q.^2);
u(2) = sum( xi .* P,1).^2 .* q.^2) / d(1);

P,2) = (xi - u(2)) .* P,1);

%
% Define higher-order polys.
%

for k = 3m+1)

d(k-1) = sum( P,k-1).^2 .* q.^2);

v(k-1) = sum( xi .* P,k-1) .* P,k-2) .* q.^2) / d(k-2);

u(k) = sum( xi .* P,k-1).^2 .* q.^2) / d(k-1);

P,k) = (xi - u(k)) .* P,k-1) - v(k-1) .* P,k-2);

end

d(m+1) = sum( P,m+1).^2 .* q.^2);

return



function bkm = bkm_comp(mmax,u,v)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for k = 0:mmax,

for m = k:mmax

% First, call bkmcheck to see if recursion formula applies

bv = bkmcheck(k,m);

if bv == -1,

bv1 = bkmcheck(k-1,m-1);
bv2 = bkmcheck(k,m-1);
bv3 = bkmcheck(k,m-2);

if bv1 == -1,

bkm1 = bkm(k,m);

else

bkm1 = bv1;

end

if bv2 == -1,

bkm2 = bkm(k+1,m);

else

bkm2 = bv2;

end

if bv3 == -1,

bkm3 = bkm(k+1,m-1);

else

bkm3 = bv3;

end

bkm(k+1,m+1) = bkm1 - u(m+1)*bkm2 - v(m)*bkm3;

else

bkm(k+1,m+1) = bv;

end

end

end

return



function bv = bkmcheck(k,m)
%
% Name: bkmcheck
%
% Usage: bv = bkmcheck(k,m)
%
% Description: Subroutine of bkm_comp -- checks
% value of k and m to see which
% formula applies
%


if k < 0 | k > m

bv = 0;

elseif k == m

bv = 1;

else

bv = -1; % Need to use full formula (or do table lookup)

end


return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


And also I wrote the code for the half-range polynoms, but I didn't add it to the hole code because I don't know how to return to the primer polynoms. The reference I saw (URL

http://www.vibetech.com/papers/paper07.pdf)

doesn't explains how to do it, neither &quot;Theoretical and Experimental Modal Analysis&quot; (Maia, Silva); and, in the second, it is said that you can return to the original polynomial coordinates with a matrix product.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%------------------------------------------------------------------------
function [P,u,v,Tac] = forsythe(xi,q,m)

L = length(xi);
P = zeros(L,m+1);
u = zeros(m+1,1);
v = zeros(m+1,1);

%
% Define zero-th order poly.
%

P,1) = ones(L,1);

%
% Define first order poly.
%

d(1) = 2*sum( q.^2);
P,2) = (xi) .* P,1);

%
% Define higher-order polys.
%

for k = 3m+1)

d(k-1) = 2*sum( P,k-1).^2 .* q.^2);

v(k-1) = 2* sum( xi .* P,k-1) .* P,k-2) .* q.^2) / d(k-2);


P,k) = (xi) .* P,k-1) - v(k-1) .* P,k-2);

end

d(m+1) = sum( P,m+1).^2 .* q.^2);

% I think that here it should be a normalization of the polynoms.

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


I tested both source codes and I found that they only works properly for one FRF with one mode. I don't know why! I generated an IRF and then I aplied fft to get the FRF. I hope reply is usefull for you. Your source code has been very usefull for me because I didn't know how to return to the original coordinates. Now, I would want to find out why I can only use it with 1 dof and 1 frf. I don't know what I'm doing wrong.

 

[mkoijn]

The denominator orthogonal poly coeffs are not the same for each frf from the structure. You have to calculate for the normal denominator coeffs which are the same for each frf.


p p
SUM[Uk]^2{B} = SUM[Uk]{Vk}
k=1 k=1

where
[Uk] = [I-[Xk]'[Xk]][BMk]^-1
{Vk} = [Xk]'{Hk}
[Xk] = -re(2[Pk*]'[Tk])
{Hk] = re(2[Pk*]'[Wk])

[Pk] = numer. polys
[Tk] = denom. ploys
[Wk] = denom. polys and data
[BMk] = orthogonal to ordinary transform

Therefor you have to use the transform matrix[BMk] before doing the least squared solution. I also rewrote the code to calculate denom. coeffs and numer. coeffs seperately, The above is taken from the richardson formenti papers on global rational fraction polynomial curve fitting.
also you are solving for the {c}'s wrong. you have to create an orthogonal matrix[Z] using 1/ denom. poly coeffs{B} as a weighting function then solve
{C} = re([Z*]'{Y})
{Y} = measurement data
see the above papers for the formula and exact formula

Best regards

finnigan

 

[maruli]

I have already understood how to undo the change of coordinates; I hadn't realised that
the polynomial are ordered for the &quot;freqs&quot; function opposite to our formulation. I applied
normalisation to the polynomials, I employed q instead of q.^2 because that agrees
with the formulation I found and I changed what you said and I could only make it work with
1 dof and 1 frf:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [A,B] = rpcoeffs(h,wi,m,n)

% Usage: [A,B] = rpcoeffs(h,wi,m,n)


% List of variables:
% m = order of numerator polynomial
% n = order of denomenator polynomial
% L = number of frequency bins to curve fit
% h = FRF matrix values at selected frequencies
% wi = vector of frequencies to use (indexed 1 to L)
% Pnum, Pden = matrices of complex orthogonal Forsythe Polynomials (one-sided)
% row i = ith frequency
% column j = (j-1)th order polynomial
% qnum = weights for numerator Forsythe poly. (m+1 x 1)
% qden = weights for denomenator Forsythe poly. (n+1 x 1)
%

%
% Normalize frequency variables
% and make FRF curve 2-sided

j = sqrt(-1);
L1 = length(wi);

wmax = wi(L1);

wi = wi/wmax;
wi = [-wi(L1:-1:1);wi];
xi = j*wi;

h = [conj(h(L1:-1:1,);h];

L = length(wi);
nh = size(h,2); % Total number of FRFs to Curve-fit

%
% Generate Orthogonal Forsythe Polynomials
%

% Numerator Polynomials

qnum = ones(L,1); % numerator poly. weights

[Pnum,Tac] = forsythe(xi,qnum,m);


for hc=1:nh
% Denomenator Polynomials
qden=conj(h,hc)).*h,hc);
[Pden,Tbd,:,hc)]=forsythe(xi,qden,n);
R,hc)=Tbd(1:n,n+1,hc);

%Form matrices for least Squares Solution.
for s=1:n
T,s)=Pden,s).*h,hc);
end
X,:,hc)=-real(Pnum'*T);

W=h,hc).*Pden,n+1);
G,hc)=real(Pnum'*W);

Ug((hc-1)*+1:hc*,=(eye-X,:,hc)'*X,:,hc))/Tbd(1:n,1:n,hc);
Vg((hc-1)*+1:hc*,1)=real(Ug((hc-1)*+1:hc*,*R,hc))-X,:,hc)'*G,hc);
end

%Solve B.
B=real((Ug'*Ug)\(Ug'*Vg));

% Solve A.
for hc=1:nh
D=real(Tbd(1:n,1:n,hc)\(B-R,hc)));
C=real(G,hc)-X,:,hc)*D);
A,hc)=real(Tac*C);
end


% Return to the original frequencies.
for hc = 1:nh
% Numerator
for k=1:m+1
A(k, = A(k, ./ wmax.^(k-1);
end
A,hc) = A(length(A,hc)):-1:1,hc);
end

% Denomenator
B=real([B;sum(Tbd(n+1,:,2))]);
for k=1:n+1
B(k) = B(k) / wmax^(k-1);
end
B) = B(length(B)):-1:1);

return




% forsythe polynomials.

function [P,Tac] = forsythe(xi,q,m)

L = length(xi);
P = zeros(L,m+1);
Tac=zeros(m+1,m+1);
u = zeros(m+1,1);
v = zeros(m+1,1);

%
% Define zero-th order poly.
%

P,1) = ones(L,1);
Tac(1,1)=1;

%
% Define first order poly.
%

d(1) = sum( P,1).^2 .* q);
u(2) = sum( xi .* P,1).^2 .* q) / d(1);
P,2) = (xi - u(2)) .* P,1);
Tac(2,2)=1;Tac(1,2)=-u(2);

%
% Define higher-order polys.
%

for k = 3m+1)

d(k-1) = sum( P,k-1).^2 .* q);

v(k-1) = sum( xi .* P,k-1) .* P,k-2) .* q) / d(k-2);

u(k) = sum( xi .* P,k-1).^2 .* q) / d(k-1);

P,k) = (xi - u(k)) .* P,k-1) - v(k-1) .* P,k-2);

Tac,k)=[0;Tac(1:m,k-1)]-u(k)*Tac,k-1)-v(k-1)*Tac,k-2);

end

d(m+1) = sum( P,m+1).^2 .* q);

% Normalization of the polynoms.
for k=1m+1)
Tac,k)=Tac,k)/(d(k)^.5);
P,k)=P,k)/(d(k)^.5);
end

return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


Alberto, Coruña, civil engineering

 

[willowmy]

Dear finnigan:

I am trying to use DIAMOND Modal Toolbox in analysing bridge vibration, but it seems that some part of the tool box does not work
For example when i try to use the Complex Exponentail method to identify modes, i get the error message below:

??? Undefined function or variable 'onHndlVector'.

Error in ==> h:\diamond_distribution\diamond\cewind.m
On line 33 ==> onHndlVector=handact(onHndlVector,hndlMatrix,9);

??? Error while evaluating uimenu Callback.


I am using MATLAB 6.5. could the problems be due to different version of MATLAB!

Iwas wondering, did you use this toolbox youself and does it work? did you face problems like the ones I am facing?

Thank you

 

[mkoijn]

Hi,

No, I never did use the Diamond toolbox as I don't have Matlab! I use GNU octave which is a freeware matlab clone. Unfortunately Octave doesn't support 3d matrices which are required in the Diamond toolbox. What I basically did was rewrite the GRFP Mfiles to work in Octave, I found however due to the method of implementation that theres a fairly low limit to the number of frf's , frequency bins and modes one can process at anyone time. I rewrote it using the Richardson Formenti implementation.

I did look at the complex Exponential method from the Diamond toolbox but haven't got round to getting it working in Octave yet. However it does say in the Diamond documentation that only the Rational fraction method is somewhat reliable and that at least Matlab 5 is needed plus signal processing toolbox.

Best regards

finnigan

 

[willowmy]

Hi finnigan:

Thank for the info.
I am having troubles with the Complex exponential method as well as the MRPT method. Aslo I amuggling with the input files becuase the documentation that come with DIAMOND do not give much details.

I will try to use it with an older version of MATLAB, may be version 5 and see what happens.

Thank again