马上注册，结交更多好友，享用更多功能，让你轻松玩转社区。
您需要 登录 才可以下载或查看，没有账号？我要加入
x
看论坛中有些同志在寻找模态参数识别的算法，我就把资源拿出来共享一下，共同学习，共同提高。声明一下，版权不是我的。
该算法实现了模态参数识别的有理多项式拟合算法，用有理多项式来拟合频响函数曲线，其中fitting函数是用来求有理分式分子分母多项式系数，result为求拟合后频响函数以及其极点、留数表示法。
fitting函数如下：
function [P,coeff]=fitting(rec,omega,phitheta,kmax)
% Scripts to compute the orthogonal polynomials required for the Rational
% Fraction Polynomial Method. (This code must be used with rfp.m)
%
% Syntax: [P,coeff]=fitting(rec,omega,phitheta,kmax)
%
% rec = FRF measurement (receptance).
% omega = frequency range vector (rad/sec).
% phitheta = weighting function (must be 1 for phi matrix or 2 for theta matrix).
% kmax = degree of the polynomial.
% P = matrix of the orthogonal polynomials evaluated at the frequencies.
% coeff = matrix used to transform between the orthogonal polynomial coefficients
% and the standard polynomial.
%
%Reference: Mark H.Richardson & David L.Formenti "Parameter Estimation from
% Frequency Response Measurements Using Rational Fraction Polynomials",
% 1篒MAC Conference, Orlando, FL. November, 1982.
%**********************************************************************
%Chile, March 2002, Cristian Andres Gutierrez Acu馻, crguti@icqmail.com
%**********************************************************************
if phitheta==1
q=ones(size(omega)); %weighting function for phi matrix
elseif phitheta==2
q=(abs(rec)).^2; %weighting function for theta matrix
else
error('phitheta must be 1 or 2.')
end
R_minus1=zeros(size(omega));
R_0=1/sqrt(2*sum(q)).*ones(size(omega));
R=[R_minus1,R_0]; %polynomials -1 and 0.
coeff=zeros(kmax+1,kmax+2);
coeff(1,2)=1/sqrt(2*sum(q));
%generating the orthogonal polynomials matrix and transform matrix
for k=1:kmax,
Vkm1=2*sum(omega.*R(:,k+1).*R(:,k).*q);
Sk=omega.*R(:,k+1)-Vkm1*R(:,k);
Dk=sqrt(2*sum((Sk.^2).*q));
R=[R,(Sk/Dk)];
coeff(:,k+2)=-Vkm1*coeff(:,k);
coeff(2:k+1,k+2)=coeff(2:k+1,k+2)+coeff(1:k,k+1);
coeff(:,k+2)=coeff(:,k+2)/Dk;
end
R=R(:,2:kmax+2); %orthogonal polynomials matrix
coeff=coeff(:,2:kmax+2); %transform matrix
%make complex by multiplying by i^k
i=sqrt(-1);
for k=0:kmax,
P(:,k+1)=R(:,k+1)*i^k; %complex orthogonal polynomials matrix
jk(1,k+1)=i^k;
end
coeff=(jk'*jk).*coeff; %complex transform matrix
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
result函数如下：
function [alpha,Residuos,Polos]=result(rec,omega,N)
% Script to compute Modal Parameter Estimation from Frequency Response
% Measurements Using Rational Fraction Polynomials Method.
%
% Syntax: [alpha,Residuos,Polos]=result(rec,omega,N)
%
% rec = FRF measurement (receptance)
% omega = frequency range vector (rad/sec).
% N = degree of freedom.
% alpha = FRF regenerated (receptance).
%
%Reference: Mark H.Richardson & David L.Formenti "Parameter Estimation from
% Frequency Response Measurements Using Rational Fraction Polynomials",
% 1篒MAC Conference, Orlando, FL. November, 1982.
%**********************************************************************
%Chile, March 2002, Cristian Andres Gutierrez Acu馻, crguti@icqmail.com
%**********************************************************************
[r,c]=size(omega);
if r<c
omega=omega.'; %omega is now a column
end
[r,c]=size(rec);
if r<c
rec=rec.'; %rec is now a column
end
nom_omega=max(omega);
omega=omega./nom_omega; %omega normalization
m=2*N-1; %number of polynomial terms in numerator
n=2*N; %number of polynomial terms in denominator
%orthogonal function that calculates the orthogonal polynomials
[phimatrix,coeff_A]=fitting(rec,omega,1,m);
[thetamatrix,coeff_B]=fitting(rec,omega,2,n);
[r,c]=size(phimatrix);
Phi=phimatrix(:,1:c); %phi matrix
[r,c]=size(thetamatrix);
Theta=thetamatrix(:,1:c); %theta matrix
T=sparse(diag(rec))*thetamatrix(:,1:c-1);
W=rec.*thetamatrix(:,c);
X=-2*real(Phi'*T);
G=2*real(Phi'*W);
d=-inv(eye(size(X))-X.'*X)*X.'*G;
C=G-X*d; %{C} orthogonal polynomial coefficients
d2N=1;
D=[d;d2N]; %{D} orthogonal polynomial coefficients
%calculation of FRF (alpha) regenerated
for n=1:length(omega),
numer=sum(C.'.*Phi(n,:));
denom=sum(D.'.*Theta(n,:));
alpha(n)=numer/denom;
end
A=coeff_A*C;
[r,c]=size(A);
A=A(r:-1:1).'; %{A} standard polynomial coefficients
B=coeff_B*D;
[r,c]=size(B);
B=B(r:-1:1).'; %{B} standard polynomial coefficients
%calculation of the poles and residues
[R,P,K]=residue(A,B);
[r,c]=size(R);
for n=1:(r/2),
Residuos(n,1)=R(2*n-1);
Polos(n,1)=P(2*n-1);
end
[r,c]=size(Residuos);
Residuos=Residuos(r:-1:1)*nom_omega; %residues
Polos=Polos(r:-1:1)*nom_omega; %poles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ 本帖最后由 ddlluu 于 2007-5-12 10:16 编辑 ] |
发表于 2007-5-12 10:13
提升卡
变色卡
显身卡
楼主
