ITD模态参数识别matlab修改版

Matlab中的系统辨识与参数估计技术Matlab（Matrix Laboratory）是一款强大的数学软件，被广泛应用于科学计算、数据处理和工程设计等领域。

在实际工程项目中，经常需要通过已有的数据来推断系统的行为模型，这就涉及到系统辨识与参数估计技术。

本文将介绍在Matlab中使用系统辨识与参数估计技术的方法和步骤。

一、系统辨识与参数估计的概念系统辨识和参数估计是在给定输入输出数据的前提下，通过数学或统计方法来推断系统的动态模型和参数值的过程。

系统辨识旨在从实验数据中提取出模型的结构信息，而参数估计则是为了获得模型的具体参数值。

二、离散时间系统的辨识与参数估计对于离散时间系统，常用的辨识方法有ARX、ARMA和ARMAX等。

以ARX 模型为例，其数学表达式为：y(t) = -a(1)y(t-1) - a(2)y(t-2) - … - a(na)y(t-na) + b(1)u(t-1) + b(2)u(t-2) + … +b(nb)u(t-nb)其中，y(t)表示系统的输出，u(t)表示系统的输入，a和b分别是系统的参数。

在Matlab中，可以使用System Identification Toolbox来进行辨识和参数估计。

首先，需要将实验数据导入到Matlab中，然后根据数据的性质选择合适的辨识方法和模型结构。

接下来，使用辨识工具箱提供的函数，通过最小二乘法或最大似然估计等算法来得到系统的参数估计值。

三、连续时间系统的辨识与参数估计对于连续时间系统，常用的辨识方法有传递函数模型、状态空间模型和灰色系统模型等。

以传递函数模型为例，其数学表达式为：G(s) = num(s)/den(s)其中，num(s)和den(s)分别是系统的分子和分母多项式。

在Matlab中，可以使用System Identification Toolbox或Control System Toolbox 来进行连续时间系统的辨识和参数估计。

ITD（Intervaled Time Domain）法是一种在时域内识别模态参数的方法，它通过对信号的时域特征进行分析，得到系统的模态参数。

在MATLAB环境下，可以使用以下步骤实现ITD法识别模态参数的编程：
1. 读入信号数据
首先需要读入需要进行模态参数识别的信号数据。

可以使用MATLAB中的`load`函数或直接从文件中读取数据。

2. 对信号进行预处理
在进行模态参数识别之前，需要对信号进行预处理，以消除噪声等干扰因素。

可以使用MATLAB中的各种信号处理函数进行滤波、去噪等操作。

3. 计算信号的时域特征
ITD法需要计算信号的时域特征，包括信号的均值、方差、峰值等。

可以使用MATLAB中的`mean`、`var`、`peak`等函数计算这些特征值。

4. 确定模态参数
根据ITD法的原理，通过对信号的时域特征进行分析，可以得到系统的模态参数。

具体来说，可以根据信号的峰值分布情况确定系统的模态频率和阻尼比等参数。

可以使用MATLAB中的`histogram`函数绘制信号的峰值分布图，并手动或自动确定模态参数。

5. 实现ITD法识别程序
将以上步骤编写成程序，即可实现ITD法识别模态参数的功能。

可以使用MATLAB中的脚本文件或函数文件进行编写。

需要注意的是，ITD法是一种比较简单的模态参数识别方法，适用于某些特定的应用场景。

在实际应用中，需要根据具体情况选择合适的模态参数识别方法。

最近在做一个项目的方案设计，应各位老总的要求，只有系统框图和器件选型可不行，为了凸显方案设计的高大上，必须上理论分析，炫一下“技术富”，至于具体有多大实际指导意义，那就不得而知了！本人也是网上一顿百度，再加几日探索，现在对用matlab 实现系统辨识有了一些初步的浅薄的经验，在此略做一小节。

必须要指出的是，本文研究对象是经典控制论理最简单最常用的线性时不变的siso 系统，而且是2阶的哦，也就是具有如下形式的传递函数：121)(22++=Ts s T s G ξ 本文要做的就是，对于有这样传递函数的一个系统，要辨识得到其中的未知数T ， ξ!!这可是控制系统设计分析的基础哦，没有系统模型，啥理论、算法都是白扯，在实际工程中非常重要哦！ 经过总结研究，在得到系统阶跃响应实验数据之后（当然如果是其他响应，也有办法可以辨识，在此还是只讨论最简单的阶跃响应实验曲线，谁让你我是菜鸟呢），利用matlab 至少可以有两种方法实现实现（目前我只会两种，呵呵）！一、函数法二、GUI 系统辨识工具箱下面分别作详细介绍！一、 函数法看官别着急，先来做一段分析（请看下面两排红*之间部分），这段分析是网上找来的，看看活跃一下脑细胞吧，如果不研读一下，对于下面matlab 程序，恐怕真的就是一头雾水咯！*******************************************************************************G(s)可以分解为：))((1)(212ωω++=s s T s G其中， [][]11112221--=-+=ξξωξξωTT1ω、2ω都是实数且均大于零。

则有：211ωω=T ，21212ωωωωξ+=传递函数进一步化为：))(()(2121ωωωω++=s s s G 因此，辨识传递函数就转化为求解1ω、2ω。

当输入为单位阶跃函数时，对上式进行拉普拉斯反变换，得系统时域下的单位阶跃响应为：tteet y 212111221)(ωωωωωωωω---+--=即 tteet y 21211122)(1ωωωωωωωω-----=-令1ω=2ωk )1(>k，得tk t ek e k k t y 22111)(1ωω-----=-⎥⎦⎤⎢⎣⎡--=---t k t e k e k k 2)1(2111ωω 对上式两边取以e 为底的对数得[]⎥⎦⎤⎢⎣⎡-+--=---t k e k t k k t y 2)1(211ln 1ln )(1ln ωω 当∞→t 时，⎥⎦⎤⎢⎣⎡---t k e k 2)1(11ln ω0→，则上式化简为[]t k k t y 21ln )(1ln ω--=-该式的形式满足直线方程b at t y +=)(*其中，)(*t y =[])(1ln t y -，1ln ,2-=-=k kb a ω)1(>k通过最小二乘算法实现直线的拟合，得到a ，b 的值，即可得到1ω、2ω的值，进而可得系统的传递函数。

基于速度响应信号的模态参数识别ITD算法研究
武晓飞;郑海起;栾军英;赵亚星
【期刊名称】《机械制造》
【年(卷),期】2009(047)009
【摘要】以环境随机激励作用下的大型塔架结构为研究背景,提出了一种基于速度响应信号的ITD方法.首先应用随机减量技术从速度随机响应信号中提取速度自由响应信号,然后推导出基于速度信号的ITD法并应用于识别系统的模态参数.该方法不用将速度信号积分成位移信号,而是直接对速度信号进行模态参数识别,避免了积分过程中引入误差.将试验结果与传统的基于位移响应的ITD法和峰值谱法进行对比,证明该方法具有较高的计算精度.
【总页数】4页(P29-32)
【作者】武晓飞;郑海起;栾军英;赵亚星
【作者单位】军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003;军械工程学院火炮工程系,石家庄,050003
【正文语种】中文
【中图分类】TH113
【相关文献】
1.基于响应协方差小波变换和SVD的结构工作模态参数识别 [J], 徐晓霞;任伟新;韩建刚
2.基于传递函数有理分式的模态参数识别算法研究 [J], 于亚婷;李敏;姜敏;邹宇
3.基于响应功率谱传递比的应变模态参数识别方法研究 [J], 曹林波; 颜王吉; 任伟新
4.基于脉冲响应数据的ARMA法建模以及模态参数识别 [J], 郭永刚;许亮华;水小平
5.基于速度信号的大型结构模态参数识别ITD方法 [J], 栾军英;张晓涛;武晓飞因版权原因，仅展示原文概要，查看原文内容请购买。

模态参数是指结构动力特性的基本参数，是描述结构动力特性的基本概念，包括固有频率、阻尼比、振型等。

结构模态参数的准确识别，是进行结构健康监测及故障诊断的重要基础，直接关系到结构安全，因此，开展结构模态参数识别技术研究具有重要的理论意义与工程实用价值。

近年来，利用环境激励已大量应用于土木工程的结构动力特性测试中。

环境激励测试能够在结构的实际工作状态下进行，更真实地了解结构的动力特性和结构性能。

本文将对各种模态识别方法进行分类汇总、论述，并对环境激励下模态参数识别算法有待进一步研究的问题进行了展望。

1频域识别算法1.1峰值拾取法基于结构的频响函数在其固有频率位置处会出现峰值的特征，可以实现对结构的模态参数识别。

由于环境激励下无法得到结构的频响函数，用功率谱密度函数代替结构的频响函数实现模态参数的识别，功率谱由实测的随机振动信号快速傅立叶变化转化得到。

姜蕾蕾[1]将幂指数窗应用于多种结构中，并与其他五种窗函数对比研究，确定能够有效改善傅立叶变换后频谱的质量,从而提高峰值拾取法的频率和阻尼比识别精度,拓宽峰值拾取法对阻尼比的适用范围。

陈涛[2]将测点传递率函数矩阵的第2阶奇异值倒数的均值为模态指示函数，建立基于多参考测点平均的峰值拾取法，准确识别系统的模态频率及振型。

在实际应用中，该方法只需计算少量的局部极值点，识别速度快，适用性广泛，被大量使用在实测实验中。

但由于峰值拾取法对峰值的选择较为敏感，对于峰值存在干扰或者峰值较小的信号，可能导致参数提取不准确，并且输出结果可能受到峰值选择的主观性影响，存在一定的不确定性。

因此，在使用时需要综合考虑实际需求和信号特征，选择合适的峰值。

1.2频域分解法频域分解法是峰值拾取法的优化算法，基本原理是根据振动响应构建谱函数矩阵，通过奇异值分解，将多自由度系统转换为单自由度体系，依靠峰值法选取特征频率，进而对系统进行识别。

频域分解法在20世纪80年代由Prevosto[3]所提出。

下面我来介绍关于系统如何构建。

对于一般的Simulink建模方法可以分为两类。

第一类是首要原则法，就是根据系统的物理意义构造微分方程，得到状态方程后，利用Simulink里自带的乘法和加法器进行模块构建，或者直接用状态方程模块分别写入A，B，C，D四个矩阵得到系统。

还有一种方法是利用Simulink里面自带的Simscape物理建模平台，里面有基本的建立机械系统，动力传动系统和机械和电子系统的模块搭建系统。

最后一种就是搭建CAD模型，利用SimMechanics把CAD模型转换成MATLAB模块进行系统的搭建。

所以这种方法需要玩家的数学功底比较强，知识面广阔。

第二类是数据驱动法，利用MAtLAB系统辨识工具箱，根据实测数据反推系统的方程，从而达到建模目的。

这种方法的首要前提是要有搭建好的物理模型进行测试。

还有一个工具叫SimulinkDesign Optimization，根据实测数据调节首要原则法建立的模型中不准确的参数。

从这里我们很快就能发现LEGO最大的优点就是很快的建立好物理模型，利用蓝牙传送被测的物理量。

所以利用系统辨识工具箱可以很快建立出模型。

下面我对系统辨识工具箱的应用进行讲解。

在讲解之前，首先我们对这次建模利用的Simulink模块进行介绍。

首先是SignalBuilder模块，如图1所示。

打开模块的设置界面，可以方便地通过GUI（用户交互式）界面画出任意需要的波形，如图2所示。

波形的输入可以通过手动的修改波形，也可以通过Excel导入的形式输入到模块内。

这个模块的意义就是写入测试数据，这个就是数据驱动法的首要条件。

当然在以后的设计中，当我们得到了系统，建立模型进行测试的时候不是通过上NXT进行测试，而是经过Simulink仿真，所以我们经常要写一些测试数据，所以这个模块的意义也变得十分重要。

在MATLAB帮助文件里面也经常会看见这个模块的出现。

然后是LEGO MINDSTORMSNXT模块，如图3所示。

实用标准文案4. 设某物理量Y 与X 满足关系式Y=aX 2+bX+c ，实验获得一批数据如下表，试辨识模型参数a ，b 和c 。

（50分）报告要求：要有问题描述、参数估计原理、程序流程图、程序清单，最后给出结果及分析。

(1)问题描述：由题意知，这是一个已知模型为Y=aX 2+bX+c ，给出了10组实验输入输出数据，要求对模型参数a ，b ，c 进行辨识。

这里对该模型参数辨识采用递推最小二乘法。

(2)参数估计原理对该模型参数辨识采用递推最小二乘法，即RLS （ recurisive least square ），它是一种能够对模型参数进行在线实时估计的辨识方法。

其基本思想可以概括为：新的估计值)(ˆk θ=旧的估计值)1(ˆ-k θ+修正项 下面将批处理最小二乘法改写为递推形式即递推最小二乘参数估计的计算方法。

批处理最小二乘估计θˆ为Y T TΦΦΦ=-1)(ˆθ，设k 时刻的批处理最小二乘估计为：k T k k T k Y ΦΦΦ=-1)(ˆθ令111)]1()()1([)()(----+-=ΦΦ=k k k P k P T kT k ϕϕ K 时刻的最小二乘估计可以表示为k T k Y k P k Φ=)()(ˆθ=)]()()[(11k y k Y k P k T k ϕ+Φ-- =)]1(ˆ)()()[()1(ˆ--+-k k k y k K k Tθϕθ；式中)()()(k k P k K ϕ=，因为要推导出P(k)和K(k)的递推方程，因此这里介绍一下矩阵求逆引理：设A 、（A+BC ）和（I +B CA 1-）均为非奇异方阵，则111111)()(------+-=+CA B CA I B A A BC A 通过运用矩（3)程序流程图(如右图1所示）递推最小二乘法（RLS）步骤如下：已知：n、b n和d。

aStep 1 ：设置初值)0(ˆθ和P(0)，输入初始数据；Step2 ：采样当前输出y(k)、和输入u(k)Step3 ：利用上面式①②③计算)(k K、)(ˆkθ和)(k P；Step4 ：k→k+1，返回step2，继续循环。

%ITDxx识别模态参数clearclcclose all hiddenformat long%% txt文件下输入fni=input('ITD法模态参数识别-输入数据文件名：','s');fid=fopen(fni,'r');mn=fscanf(fid,'%d',1);%模态阶数%定义输入实测数据类型%ig=1时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果%ig=2频域数据如频响函数实部和虚部数据ig=fscanf(fid,'%f',1);%ig=1时，f为采样频率sf，ig=2时，f为频率间隔dff=fscanf(fid,'%f',1);fno=fscanf(fid,'%s',1);%输出数据文件名b=fscanf(fid,'%f',[ig,inf]);%实测时域或频域数据status=fclose(fid);%%clc;clear all;format long[FileName,PathName] = uigetfile('*.mat', 'Select the Mat-files of time signal'); %窗口读文件,并获取包含路径的文件名if isequal(FileName,0)disp('User cancel the selection');%如果取消选择则显示提示return;elseFULLFILE=fullfile(PathName,FileName);Signal_str= sprintf('User selected signal file:%s',FULLFILE);disp(Signal_str);Struct=load(FULLFILE);endc=fieldnames(Struct);%得到一个元胞数组，包含Struct中各个域名（倘若有多个的话）b=getfield(Struct,c{1}); %获取c{1}对应的域中的内容b=b(3601:9600);%%%ig=1时域数据如冲击响应、自由振动、互相关函数、随机减量法处理结果%ig=2频域数据如频响函数实部和虚部数据ig=input('数据类型ig=');f=input('采样频率f=');%指定采样频率mn=input('计算模态阶数mn=');%指定计算模态阶数%建立特征方程矩阵的阶数（为模态阶数的2倍）nm=2*mn;%组织识别计算多用的时域数据及参数if ig==1%实测时域数据sf=f;%采样频率n=fix(length(b)/2);%向0靠拢取整，取时域数据的长度h=b(1,1:2*n)';%将输入时域数据赋值给列向量hdt=1/sf;%时间间隔t=0:dt:(2*n-1)*dt;%建立离散时间向量else %实测频域数据df=f;%取频率间隔n=length(b(1,:));f=0:df:(n-1)*df;%建立离散频率向量H=b(1,:)'+b(2,:)'*i;%建立对应正负频率的实测频响函数向量H(n+1)=real(H(n));H(n+2:2*n)=conj(H(n:-1:2));%conj求负数的共轭值h=real(ifft(H));%频响函数经IFFT并取实部变换成脉冲响应函数t=linspace(0,1/df,2*n);%建立离散时间向量dt=t(2)-t(1);%计算时间间隔end%计算自由振动响应矩阵L=length(h);M=L/2;for k=1:nmx1(k,:)=h(k:L-(nm-k+1))';x2(k,:)=h(k+1:L-(nm-k))';end%用最小二乘法求解特征方程矩阵B=x1\x2;%B=x2*x1'*inv(x1*x1');[A,V]=eig(B);%计算特征值及特征向量（特征值V，特征向量A）%变换特征值对角阵为一向量for k=1:nmU(k)=V(k,k);endF1=abs(log(U'))./(2*pi*dt);%计算模态频率向量D1=sqrt(1./(((imag(log(U'))./real(log(U'))).^2)+1));%计算阻尼比向量%计算振型系数向量l=1;for k=0:(2*n-1)Va(k+1,:)=[conj(U).^k];endS1=(inv(conj(Va')*Va)*conj(Va')*h);%inv矩阵求逆h1=real(Va*S1);%计算生成的脉冲响应函数%绘制脉冲响应函数拟合曲线图figure(1);plot(t,h,':',t,h1);xlabel('时间（s）');ylabel('幅值');legend('实测','拟合');grid on;if ig>1H1=fft(Va*S1);%计算生成的频响函数%绘制频响函数实部拟合曲线图figure(2);nn=1:n;subplot(2,1,1);plot(f,real(H(nn)),':',f,real(H1(nn)),'r'); xlabel('频率(Hz)');ylabel('实部');legend('实测','拟合');grid on;%绘制频响函数虚部拟合曲线图subplot(2,1,2);plot(f,b(2,:),':',f,imag(H1(nn)),'r');xlabel('频率(Hz)');ylabel('虚部');legend('实测','拟合');grid on;end[F2,I]=sort(F1);%将自振频率从小到大排列%剔除方程中的非模态项（非共轭根）和共轭项（重复项）m=0;for k=1:1:nm-1if F2(k)~=F2(k+1)continue;endm=m+1;l=I(k);F(m)=F1(l);%自振频率D(m)=D1(l);%阻尼比S(m)=S1(l);%振型系数end%打开文件输出识别的模态参数数据fno='out.txt';fid=fopen(fno,'w');fprintf(fid,'频率（Hz）阻尼比（%%）振型系数\n');for k=1:mfprintf(fid,'%10.4f %10.4f%10.6f\n',F(k),D(k)*100.0,imag(S(k))); endstatus=fclose(fid);。