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

weixin

　　% 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的博客，程序版权见程序注释