传感器与测试技术-基于MATLAB的二阶系统固有频率和阻尼测量方法

完整版的程序下载：

0 系统采集背景

用力锤对时间进行敲击，产生一个宽频带的激励，它能在很宽频率范围内激励出各种模态。锤击力函数及频谱如图3所示。所以在锤击法中的脉冲，其应用频率的主瓣应尽量宽。因此要根据被激振结构选择不同的锤帽，以调节应用频率的范围。脉冲锤击法测量系统的示意图如图4所示。采用脉冲锤击法时，为了消除噪声干扰，必须采用多次测量实验，结果取平均值。

图3 锤击脉冲函数及其频谱

图4 脉冲锤击法的测量系统

由于采用快速正弦扫描激振法需要用到激振器，结构复杂，成本大，因此采用操作简单的锤击法对水泥混凝土结构件的阻尼特性进行测量。

1 问题描述

采用锤击法，获得了水泥试件的响应。

响应数据中：第一列为采样时间，第二列为响应（加速度信号）

2 数据源

系统数据如下： 

3 方法

FFT 求取固有频率（看前面的例子）

时域极值比或者半功率法求取阻尼

3.1 时域法

二阶系统的单位脉冲响应传递函数为：

          （4）

二阶系统的特征方程：

                   （5）

由此得两个极点：

                （6）

由于测量系统中水泥混凝土材料结构件的阻尼系数 ，二阶系统的极点是一对共轭复根。

                 

整理后得到：

                （9）

上式经过拉氏反变换，得：

                 （10）

图5为二阶系统的时域衰减振动波形图，分析一下振动波形峰值的变化规律可知，当 ，振动波形图出现第一个峰值 ：

               （12）

对上式两端取自然对数，整理得：

                   （14）

由式（14）可知，在系统的固有频率 确定的条件下，阻尼系数可通过时域波形的两个峰值计算获得。

图5 衰减振动的波形

3.2 半功率谱法

由上文已知阻尼系数测量系统传递函数（式3），通过FFT运算（快速傅里叶变换）可以得出：

                 （15）

（不能直接贴入公式，我截图）

 2.3 结果分析

 

 前期用LabVIEW写的，   空了补充 MATLAB 的算法。

3 MATLAB 分析

3.1 数据导入

  导入数据，导入为矩阵。

在workspace把变量名称改为Data。

3.2 数据分析

数据分析， 3次锤击结果。

% analyze the impuse data and find epsilon and wn

dt=(Data(5,1)-Data(2,1))/3 % sample time interval
fs=1/dt  %sample frequency, 5kHz

figure
plot(Data(:,1),Data(:,2))
xlabel('Time/s')
ylabel('Amplitude/V')

fs=5kHz。

锤击作用时间<0.2s.

取前0.15s数据进行分析。 标准的 单位脉冲响应

3.3 FFT

我选择前面的第一个锤击数据进行分析，如果是实际应用，需要将3次结果进行平均，自动搜索出锤击数据（找大于一定阈值作为起始点，截取0.15s长的数据），或者通过鼠标框来选择需要分析的区域。也就是程序中的samples范围。

%%FFT
% FFT 
samples=ceil(0.1/dt); % get 0.1 s samples
ydata_fft=fft(Data(1:samples,2));             % 傅里叶变换 ,
ydata_abs=2*abs(ydata_fft(1:samples/2))';% 取绝对值
f = fs*(0:(samples/2)-1)/samples;% 频率轴
 figure
plot(f,ydata_abs)
xlabel('Frequency/Hz');ylabel('Amplitude/mv');
title('Frequnecy domain')
set(gca,'xscale','log') % log in x

 

3.4 半功率法

按照定义来搜索 半功率的两个点。

以下程序作了简化：1 第一个点是搜> 0.707A; 2.第二个点是将最大值右边数组去除，搜≤0.707A。 

实际使用，还需要进行插值（根据搜寻到的频率及其幅值，在相邻点进行线性插值），以定位半功率点。因为搜索到的点不会正好是半功率点。

%% caculate the parameters
% find maximum Ampliude
[A,idx]=max(ydata_abs) %find A and its corresponding index
%idx is the index of maxium amplitude
f0=f(idx) % 共振频率


%find 1/sqrt(2) of A.
A3db=A/sqrt(2)

%search A3db's corresponding index

idx1=find(ydata_abs>=A3db) % find w1,

f1=f(idx1(1)) 

% search w2

idx2=find(ydata_abs(idx:end)<=A3db) % find w2 
f2=f(idx2(1)+idx)

%caculate episilon, damping ratio

epsilon=(f2-f1)/(2*f0)*100

% replot the frequency figure

figure
plot(f,ydata_abs)
xlabel('Frequency/Hz');ylabel('Amplitude/mv');
title('Frequnecy domain')
set(gca,'xscale','log') % log in x
hold on
 f=[f1,f0,f2];
 A2=[A3db,A,A3db];
 plot(f,A2,'o')
 hold off

 

最终结果！！

按照定义来理解！