[转]振动试验数据整理和处理程序

unknowno

转自：http://blog.sina.com.cn/u/462c63480100097s


%   STDATA.m
%   振动台试验数据整理和处理软件
%   Copyright：Concrete
%        江苏省土木工程与防灾减灾重点实验室
%本程序部分代码使用了王济、胡晓编著的《MATLAB在振动信号处理中的应用》一书中的
%代码,在此表示感谢；
%注意：
%  1. 因编程时间仓促，本软件为非通用软件，内部参数与具体试验相关;
%  2. 使用者在使用时应确保明了程序含义，并根据试验内容作相应修改;
%  3. 本软件为自由软件，使用者在注明出处后可以自由使用；
%  4. 对使用本程序造成的任何后果,本人不承担任何责任；
%  5. 欢迎指出程序错误之处，联系方法：concretelu@gmail.com。
%
%  数据文件示例：
%     2  ％需要处理的文件的数量
%     AE1a.C02  ％以下为需要处理的文件名
%     AE3a.C02
%     128  ％采样频率
%     0  ％是否需要进行峰值调整；0:不进行；
%     -1   ％是否需要进行基线调整：－1:进行但不输出文件
%     3  ％基线调整拟合多项式阶数
%     0  ％是否需要进行滤波0:不进行；
%     1  ％是否需要进行积分：1:进行且输出文件；
%     0.5  ％最小截止频率
%     40   ％最大截止频率
%     9800  ％单位转换系数
%     2  ％积分次数(加速度求位移）
%     1  ％是否需要进行通道间计算
%     0  ％是否需要进行绘图

clear                %清除内存中所有变量和函数
clc                 %清除工作窗口中所显示的内容
close all hidden          %关闭所有隐藏的窗口
echo off;              %关闭回显
fnc=input('数据处理控制文件名:','s');
fid=fopen(fnc,'r');
nf=fscanf(fid,'%f',1);       %需处理的文件数量
for i=1:nf
  fname(i,:)=fscanf(fid,'%s',1); %读入各数据文件名
end
sf=fscanf(fid,'%f',1);       %读入采样频率
flag_adj=fscanf(fid,'%d',1);
if flag_adj~=0           
  %是否需要进行峰值调整；0:不进行；1:进行且输出文件；－1:进行但不输出文件
  adj=fscanf(fid,'%f',[1,nf]);  %读入各文件调整系数
end
flag_bas=fscanf(fid,'%d',1);
if flag_bas~=0           
  %是否需要进行基线调整0:不进行；1:进行且输出文件；－1:进行但不输出文件
  basm=fscanf(fid,'%d',1);    %读入拟合多项式阶数
end
flag_fil=fscanf(fid,'%d',1);
if flag_fil~=0           
  %是否需要进行滤波0:不进行；1:进行且输出文件；－1:进行但不输出文件
  fmin=fscanf(fid,'%f',1);    %最小截止频率
  fmax=fscanf(fid,'%f',1);    %最大截止频率
end
flag_int=fscanf(fid,'%d',1);
if flag_int~=0           
  %是否需要进行积分：0:不进行；1:进行且输出文件；－1:进行但不输出文件
  fmin=fscanf(fid,'%f',1);    %最小截止频率
  fmax=fscanf(fid,'%f',1);    %最大截止频率
  c=fscanf(fid,'%f',1);      %单位转换系数
  it=fscanf(fid,'%f',1);     %积分次数
end
flag_cmp=fscanf(fid,'%d',1);    %是否需要进行通道间计算
flag_fig=fscanf(fid,'%d',1);    %是否需要进行绘图

status=fclose(fid);         %关闭控制文件
delete('Acc.xls');         %删除存在的结果文件
delete('Disp.xls');
delete('Drift.xls');
delete('Model.xls');
delete('TotalAcc.xls');
delete('TotalDisp.xls');
delete('TotalDrift.xls');
Acc_string={'项目','台面加速度','台面位移','底层','二层','三层',...
  '四层','五层','六层','七层','八层','九层','屋面','一层南','一层北',...
  '六层南','六层北','屋面南','屋面北'};        %建立加速度表头字符串
Disp_string={'项目',,'底层','二层','三层','四层','五层','六层','七层','八层',...
  '九层','屋面','一层南','一层北','六层南','六层北','屋面南','屋面北'};
                           %建立位移表头字符串
Drift_string={'项目','底层','二层','三层','四层','五层','六层','七层',...
  '八层','九层','一层南北','六层南北','屋面南北'};  %建立层间位移表头字符串
for fi=1:nf
  x=dlmread(fname(fi,:),'\t',1,1); %从数据文件读入数据,去除标题行和时间列
  [n,ch]=size(x);          %确定数据长度n和通道数ch
  t=0:1/sf:(n-1)/sf;        %建立信号离散时间向量t
  TableDisp=x(:,2);         %读取第二通道：台面位移
  if flag_adj~=0          %确定是否需要进行加速度峰值调整
    x=adj(fi).*x;         %进行峰值调整
    if flag_adj==1
      filename=strcat('Adj_',fname(fi,:));
      dlmwrite(filename,x);   %输出峰值调整后的数据
    end
   
    figure(1);           %绘制调整后的加速度信号
    subplot(gl,4,[1 4]);      %采用两栏显示台面加速度   
    plot(t,x(:,1));
    xlabel('时间 (s)');
    ylabel('加速度 (g)');
    title('台面加速度信号');
    grid on;
    for i=3:ch;           %绘制各通道信号调整峰值后随时间变化的图形
      subplot(gl,4,i+2);   
      plot(t,x(:,i));
      grid on;
    end
  end
  
  case_name=fname(fi,:);       %建立工况名称向量
  x2=x.^2;              %计算加速度有效值
  xsum=cumsum(x2);
  rmsx=sqrt(xsum(n,:)/n);
  xlswrite('Acc.xls',Acc_string,case_name,'A1'); %输出加速度时程及统计值
  xlswrite('Acc.xls',{'Max Acc.';'Min Acc.';'MaxAbs Acc.';'Val Acc.'},...
    case_name,'A2');
  acc_out=[max(x);min(x);max(-min(x),max(x));rmsx];
  xlswrite('Acc.xls',acc_out,case_name,'B2');
  xlswrite('Acc.xls',[t',x],case_name,'A6');
  
  MaxAcc=[MaxAcc;max(x)];         %建立每个工况的统计值,后期统一输出
  MinAcc=[MinAcc;min(x)];
  AbsMaxAcc=[AbsMaxAcc;max(-min(x),max(x))];
  ValAcc=[ValAcc;rmsx];
  
  if flag_bas~=0            %确定是否需要进行基线调整
    for i=1:ch
      tt=(0:1/sf:(n-1)/sf)';    %建立信号离散时间列向量t
      a=polyfit(tt,x(:,i),basm);  %计算趋势项的多项式待定系数向量a
      x(:,i)=x(:,i)-polyval(a,tt); %形成消除趋势项的数据
    end
    if flag_bas==1
      filename=['Bas_',fname(fi,:)];
      dlmwrite(filename,x);     %输出基线调整后的数据
    end
  end
  
  if flag_fil~=0           %确定是否需要进行单独滤波
    fmin=fscanf(fid,'%f',1);    %最小截止频率
    fmax=fscanf(fid,'%f',1);    %最大截止频率
  end
  
  if flag_int~=0           %确定是否需要进行积分
    tt=(0:1/sf:(n-1)/sf)';     %建立信号离散时间列向量t
    nfft=2^nextpow2(n);       %取大于并接近n的2的幂次方为FFT长度
    df=sf/nfft;           %计算频率间隔(Hz/s)
    ni=round(fmin/df+1);      %计算指定频带对应数组的下标
    na=round(fmax/df+1);
    dw=2*pi*df;           %计算圆频率间隔(rad/s)
    w1=0:dw:2*pi*(0.5*sf-df);    %建立正的离散圆频率向量
    w2=2*pi*(0.5*sf-df):-dw:0;   %建立负的离散圆频率向量
    w=[w1,w2];           %将正负圆频率向量合成单一向量
    w=w.^it;            %以积分次数为指数,建立圆频率变量向量
   
    disp=zeros(n,ch);        %为避免不同向量长度影响,变量清零
    adisp=zeros(n,ch);
    Drift=zeros(n,ch-6);
   
    for i=1:ch
      tx=x(:,i)';
      y=fft(tx,nfft);       %FFT变换
      a=zeros(1,nfft);      %进行积分的频域变换
      a(2:nfft-1)=y(2:nfft-1)./w(2:nfft-1);
      if it==2          %进行二次积分的相位变换
        y=-a;
      else
        real(y)=imag(a);    %进行一次积分的相位变换
        imag(y)=-real(a);
      end
      a=zeros(1,nfft);
      a(ni:na)=y(ni:na);     %消除指定正频带外的频率成分
      a(nfft-na+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
                    %消除指定负频带外的频率成分
      y=ifft(a,nfft);       %IFFT变换
      y=real(y(1:n))*c;      
            %取逆变换的实部n个元素并乘以单位变换系数作为积分结果
      ab=polyfit(tt,y',2);        %计算趋势项的多项式待定系数向量a
      adisp(:,i)=y'-polyval(ab,tt);    %形成消除趋势项的数据(绝对位移)
      disp(:,i)=adisp(:,i)-TableDisp;   %计算各点的相对位移
      ac=polyfit(tt,disp(:,i),2);     %计算趋势项的多项式待定系数向量ac
      disp(:,i)=disp(:,i)-polyval(ac,tt); %形成消除趋势项的数据(相对位移)
    end
  end
  
  if flag_int==1             %如果需要输出位移文件
    y1=min(disp(:,3:ch));        %搜寻每层相对位移最小值（负向）
    y2=max(disp(:,3:ch));        %搜寻每层相对位移最大值（正向）
    d2=disp(:,3:ch).^2;
    dsum=cumsum(d2);
    rmsd=sqrt(dsum(n,:)/n);       %搜寻每层相对位移有效值
    xlswrite('Disp.xls',Disp_string,case_name,'A1'); %输出位移时程及统计值
    xlswrite('Disp.xls',{'Max Disp';'Min Disp';'MaxAbs Disp';'Val Disp'},...
      case_name,'A2');
    dispout=[y2;y1;max(-y1,y2);rmsd];
    xlswrite('Disp.xls',dispout,case_name,'B2');
    xlswrite('Disp.xls',[t' disp(:,3:ch)],case_name,'A6');
                  %输出基线调整后的相对位移数据
    MaxDisp=[MaxDisp;y2];
    MinDisp=[MinDisp;y1];
    AbsMaxDisp=[AbsMaxDisp;max(-y1,y2)];
    ValDisp=[ValDisp;rmsd];
  end
  
  if flag_cmp~=0           %确定是否需要进行通道间计算
    for i=1:9
      Drift(:,i)=disp(:,i+3)-disp(:,i+2); %计算层间位移
    end
    Drift(:,10)=disp(:,14)-disp(:,13);   %计算两角点位移差
    Drift(:,11)=disp(:,16)-disp(:,15);
    Drift(:,12)=disp(:,18)-disp(:,17);
  end
  
  if flag_cmp==1
    y3=min(Drift);         %搜寻各层间位移最小值
    y4=max(Drift);         %搜寻各层间位移最大值
    dr2=Drift.^2;
    drsum=cumsum(dr2);
    rmsdr=sqrt(drsum(n,:)/n);    %计算层间位移有效值
    xlswrite('Drift.xls',Drift_string,case_name,'A1');
    xlswrite('Drift.xls',{'Max Drift';'Min Drift';'MaxAbs Drift';'Val Drift'},...
      case_name,'A2');
    driftout=[y4;y3;max(-y3,y4);rmsdr];
    xlswrite('Drift.xls',driftout,case_name,'B2');
    xlswrite('Drift.xls',[t' Drift],case_name,'A6');
                    %输出层间位移数据
    MaxDrift=[MaxDrift;y4];
    MinDrift=[MinDrift;y3];
    AbsMaxDrift=[AbsMaxDrift;max(-y3,y4)];
    ValDrift=[ValDrift;rmsdr];
  end
  [Txy,F]=tfestimate(x(:,3),x(:,12),nfft/2,[],[],40);
                    %求台面信号和顶层信号间的传递函数
  %aTxy=abs(Txy);            %输出模态幅值
  %xlswrite('Model.xls',{'频率','幅值','相位角'},case_name,'A1');
  %xlswrite('Model.xls',[F aTxy angle(Txy)],case_name,'A2');

  if flag_fig~=0           %确定是否需要进行绘图
    gl=ceil((ch-2)/4)+1;      %计算图形显示所需要的行数，列数为4
    figure(2);
    subplot(gl,4,[1 4]);      %采用两栏显示台面加速度   
    plot(t,adisp(:,1),t,TableDisp,'r');
    xlabel('时间 (s)');
    ylabel('位移 (mm)');
    title([case_name,' 积分位移信号']);
    grid on;
    ymin=floor(min(y1'));      %搜寻全部相对位移最小值
    ymax=ceil(max(y2'));      %搜寻全部相对位移最大值
    xmax=ceil((n-1)/sf);
    for i=3:ch;           %绘制各通道相对位移随时间变化的图形
      subplot(gl,4,i+2);   
      plot(t,disp(:,i));
      axis([0  xmax ymin ymax]);
      %xlabel('时间  (s)');
      %ylabel('加速度  (g)');
      %title('通道'+int2str(ch));
      grid on;
    end
    figure(3);
    ymin=floor(min(y3'));      %搜寻全部相对位移最小值
    ymax=ceil(max(y4'));      %搜寻全部相对位移最大值
    xmax=ceil((n-1)/sf);
    for i=1:12;           %绘制各通道信号x随时间变化的图形
      subplot(4,3,i);     
      plot(t,Drift(:,i));
      axis([0  xmax ymin ymax])
      grid on;
    end
    figure(4);
    subplot(2,1,1);
    plot(F,aTxy);
    grid on;
    subplot(2,1,2);
    plot(F,angle(Txy)*180/pi());         %输出相位角
    grid on;
  end
end
xlswrite('TotalAcc.xls',Acc_string,'MaxAcc','A1');   %建立数据表名称
xlswrite('TotalAcc.xls',Acc_string,'MinAcc','A1');
xlswrite('TotalAcc.xls',Acc_string,'AbsMaxAcc','A1');
xlswrite('TotalAcc.xls',Acc_string,'ValAcc','A1');
xlswrite('TotalDisp.xls',Disp_string,'MaxDisp','A1');
xlswrite('TotalDisp.xls',Disp_string,'MinDisp','A1');
xlswrite('TotalDisp.xls',Disp_string,'AbsMaxDisp','A1');
xlswrite('TotalDisp.xls',Disp_string,'ValDisp','A1');
xlswrite('TotalDrift.xls',Drift_string,'MaxDrift','A1');
xlswrite('TotalDrift.xls',Drift_string,'MinDrift','A1');
xlswrite('TotalDrift.xls',Drift_string,'AbsMaxDrift','A1');
xlswrite('TotalDrift.xls',Drift_string,'ValDrift','A1');
xlswrite('TotalAcc.xls',MaxAcc,'MaxAcc','B2');   %输出数据
xlswrite('TotalAcc.xls',MinAcc,'MinAcc','B2');
xlswrite('TotalAcc.xls',AbsMaxAcc,'AbsMaxAcc','B2');
xlswrite('TotalAcc.xls',ValAcc,'ValAcc','B2');
xlswrite('TotalDisp.xls',MaxDisp,'MaxDisp','B2');
xlswrite('TotalDisp.xls',MinDisp,'MinDisp','B2');
xlswrite('TotalDisp.xls',AbsMaxDisp,'AbsMaxDisp','B2');
xlswrite('TotalDisp.xls',ValDisp,'ValDisp','B2');
xlswrite('TotalDrift.xls',MaxDrift,'MaxDrift','B2');
xlswrite('TotalDrift.xls',MinDrift,'MinDrift','B2');
xlswrite('TotalDrift.xls',AbsMaxDrift,'AbsMaxDrift','B2');
xlswrite('TotalDrift.xls',ValDrift,'ValDrift','B2');
close all;

[ 本帖最后由 unknowno 于 2007-7-6 20:56 编辑 ]

unknowno

可以对振动试验数据进行批量处理，如滤波、消除趋势项、积分等！
非常方便！

vib

看着挺复杂的，但是还看不懂

linlyjiang

还是不会用！

baobao1982

好多东西已经过时了
它们处理一般的问题还行
换成实际的东西就不管用了
就拿消除趋势线过来说吧
实际中的心混入的趋势项是线性的或是多项式的情况有多大？
呵呵，没有多大
大多情况下我们根被就不知道趋势项的具体形式
那又怎么用最小二乘法去除呢？
个人愚见
仅供参考

vibration2008

原帖由 baobao1982 于 2009-3-1 14:43 发表

                               
登录/注册后可看大图

好多东西已经过时了
它们处理一般的问题还行
换成实际的东西就不管用了
就拿消除趋势线过来说吧
实际中的心混入的趋势项是线性的或是多项式的情况有多大？
呵呵，没有多大
大多情况下我们根被就不知道趋势项的 ...

你有更好的程序或者方法吗？

zw393778182

老师给了我三组台架振动数据 要我用MATLAB进行预处理 时域分析 频域分析 以及可视化输出 按照MATLAB在振动信号处理中的应用一书编出的程序运行不了

dreamstone

二次积分个人觉得程序无法解决，最近偶在对加速度进行二次积分得位移，用matlab指令，自己通过梯形、5点公式、差分法，都试过了，差异还是比较大，个人愚见。
楼主给扫盲，或者给解决下，谢谢

wanyeqing2003

楼主的分析有针对性。别人使用，还要熟悉半天。

zhoulid2010

学学，挺好

0810064

对楼主的辛勤工作表示感谢，好贴挺起来

yyy101102

对楼主的辛勤工作表示感谢