共享有理分式法识别模态参数的程序
看论坛中有些同志在寻找模态参数识别的算法，我就把资源拿出来共享一下，共同学习，共同提高。声明一下，版权不是我的。该算法实现了模态参数识别的有理多项式拟合算法，用有理多项式来拟合频响函数曲线，其中fitting函数是用来求有理分式分子分母多项式系数，result为求拟合后频响函数以及其极点、留数表示法。
fitting函数如下：
function =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: =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=; %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=;
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 =result(rec,omega,N)
% Script to compute Modal Parameter Estimation from Frequency Response
% Measurements Using Rational Fraction Polynomials Method.
%
% Syntax: =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
%**********************************************************************
=size(omega);
if r<c
omega=omega.'; %omega is now a column
end
=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
=fitting(rec,omega,1,m);
=fitting(rec,omega,2,n);
=size(phimatrix);
Phi=phimatrix(:,1:c); %phi matrix
=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} 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;
=size(A);
A=A(r:-1:1).'; %{A} standard polynomial coefficients
B=coeff_B*D;
=size(B);
B=B(r:-1:1).'; %{B} standard polynomial coefficients
%calculation of the poles and residues
=residue(A,B);
=size(R);
for n=1:(r/2),
Residuos(n,1)=R(2*n-1);
Polos(n,1)=P(2*n-1);
end
=size(Residuos);
Residuos=Residuos(r:-1:1)*nom_omega; %residues
Polos=Polos(r:-1:1)*nom_omega; %poles
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[ 本帖最后由 ddlluu 于 2007-5-12 10:16 编辑 ]
给出一个示例数据
FRF_NOISE.txt 为混有噪声的频响函数，为拟合效果图检验程序：
plot(w./(2*pi),20*real(log10(FRF_noise)),'b'), hold on,
omega=w(50:450,1);
rec=FRF_noise(50:450,1);
N=3; =rfp(rec,omega,N);
plot(omega./(2*pi),20*real(log10(alpha)),'r'), hold off, 你这个程序，有噪音的话误差很大吧？ 找了好久,十分感谢.正在研究中.有不懂的地方还要向各位请教. 程序的算法哪里有介绍比较详细的?看了很多,都叫粗糙.另外运行了程序,总是出错.ddlluu能否给出更详细的程序执行过程呢?我matlab也是学习的不深入,谢谢了. 我现在很需要这样的程序，但是我看不懂，希望楼主对程序多加一些注释 Fraction Polynomial Method. (This code must be used with rfp.m)
rpf.m 文件呢？？？楼主 哪位高手能讲解一下，上面的这段程序？？不怎么会用啊 全部程序和理论背景可从如下网址下载：
http://www.mathworks.com/matlabcentral/fileexchange/3805
回复 9楼 handful 的帖子
十分感激，不需要注册吧？ 谢谢楼主，但是从结果看，并不是很理想啊...... close all;clear all;clc;load FRF_noise.mat
plot(FRF_noise(:,1),20*log10(abs(FRF_noise(:,2))),'b'), hold on,
xlabel('Frequency (rad/sec)'),ylabel('Magnitude (dB)'),
omega=FRF_noise(50:450,1);
rec=FRF_noise(50:450,2);
N=3;
=rfp(rec,omega,N);
fn=par(:,1) %natural frequencies
xi=par(:,2) %damping ratios
C= %modal constant (magnitud,phase)
plot(omega,20*log10(abs(alpha)),'r'), hold off,
其中rec 和omega 是什么意思？ N值说是自由度，等于模态阶数不？ 大神，有了频响函数的数据，能否模态参数（模态刚度，质量，阻尼比）具体怎么实现，急求 高手，真心不会，学习 ddlluu 发表于 2007-5-12 10:37
FRF_NOISE.txt 为混有噪声的频响函数，未命名.jpg 为拟合效果图
检验程序：
plot(w./(2*pi),20*real(log1 ...
权限不够，看不到
页:
[1]