声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1395|回复: 1

振动台试验数据处理程序,内含超详细的注释

[复制链接]
发表于 2018-4-26 09:42 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
  % 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;

  本文转载自新浪Glulam的博客,程序版权见程序注释

回复
分享到:

使用道具 举报

您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-11-14 12:40 , Processed in 0.066569 second(s), 20 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表