modal parameter code, run out of memory
modal parameter code, run out of memory
(OP)
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 + 1:(m+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 = 3:(m+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
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 + 1:(m+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 = 3:(m+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
RE: modal parameter code, run out of memory
M
RE: modal parameter code, run out of memory
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;
RE: modal parameter code, run out of memory
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
RE: modal parameter code, run out of memory
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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)*(n)+1:hc*(n),1)=G(:,hc);
Ug((hc-1)*(n)+1:hc*(n),:)=-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 = 3:(m+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 "Theoretical and Experimental Modal Analysis" (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 = 3:(m+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.
RE: modal parameter code, run out of memory
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
RE: modal parameter code, run out of memory
the polynomial are ordered for the "freqs" 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)*(n)+1:hc*(n),:)=(eye(n)-X(:,:,hc)'*X(:,:,hc))/Tbd(1:n,1:n,hc);
Vg((hc-1)*(n)+1:hc*(n),1)=real(Ug((hc-1)*(n)+1:hc*(n),:)*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 = 3:(m+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=1:(m+1)
Tac(:,k)=Tac(:,k)/(d(k)^.5);
P(:,k)=P(:,k)/(d(k)^.5);
end
return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Alberto, Coruña, civil engineering
RE: modal parameter code, run out of memory
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
RE: modal parameter code, run out of memory
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
RE: modal parameter code, run out of memory
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