振动台试验数据处理程序，内含超详细的注释
% 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的博客，程序版权见程序注释
声明：本微信转载文章出于非商业性的教育和科研目的，并不意味着赞同其观点或证实其内容的真实性。版权归原作者所有，如转载稿涉及版权等问题，请立即联系我们，我们会予以更改或删除相关文章，保证您的权利！
