| % 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的博客，程序版权见程序注释 |
GMT+8, 2024-11-25 17:51 , Processed in 0.066894 second(s), 22 queries , Gzip On.
Powered by Discuz! X3.4
Copyright © 2001-2021, Tencent Cloud.
