加速度转换成位移的matlab代码及说明

加速度转换成位移的matlab代码及说明
由测量的加速度离散数据数据转化成位移数据一般不直接在时域进行积分处理，而是由时域转换成频域在频域中进行二次积分再转化到时域中得到位移结果。

相关matlab处理程序方法参考王济老师的matlab在振动信号处理中的应用中的程序如下：
%频域积分
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc; close all hidden
%%%%%%%%%%%%%%%%%%
fni=input('频域积分-输入数据文件名：','s');
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);%采样频率
fmin=fscanf(fid,'%f',1);%最小截止频率
fmax=fscanf(fid,'%f',1);%最大截止频率
c=fscanf(fid,'%f',1);%单位变换系数
it=fscanf(fid,'%f',1);%积分次数
sx=fscanf(fid,'%s',1);%横向坐标轴的标注
sy1=fscanf(fid,'%s',1);%纵向坐标轴输入单位的标注
sy2=fscanf(fid,'%s',1);%纵向坐标轴输出单位的标注
fno=fscanf(fid,'%s',1);%输出数据文件名
x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
status=fclose(fid);
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*(0.5*sf-df);
%建立负的离散圆频率向量
w2=2*pi*(0.5*sf-df):-dw:0;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
%以积分次数为指数，建立圆频率变量向量
w=w.^it;
%进行积分的频域变换
a=zeros(1,nfft); a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);
if it == 2
y=-a; %进行二次积分的相位变换
else
a1=imag(a); a2=real(a); y=a1-a2*i; %进行一次积分的相位变换
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变换
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
subplot(2,1,1); plot(t,x); xlabel(sx); ylabel(sy1); grid on; %绘制几分钱的时程曲线图形
subplot(2,1,2); plot(t,y); xlabel(sx); ylabel(sy2); grid on; %绘制积分后的时程曲线图形
%打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n, fprintf(fid,'%f \n',y(k)); end
status=fclose(fid);
程序使用说明：运行该程序会显示频域积分-输入数据文件名：在之后写出输入数据文件的完整文件名（包括扩展名，例如data.txt）要求数据文件必须在matlab工作空间路径下，该数据文件必须满足如下格式要求：
采样频率
下限截止频率
上限截止频率
单位变换系数
积分阶数（加速度变位移为2，变速度为1）
时间（s）
加速度（g）
位移（mm）
out.txt
具体加速度数据（回车或空格隔开）
下图为书上给出的例子：
满足这样的格式的数据文件才能正常处理，注意单位变换系数和坐标有关，如果按书上的例子加速度以g为单位位移以mm为单位，单位变换系数应该为9.8*1000=9800；而用米每二次方秒为加速度单位以mm为位移单位的话就是1*1000=1000；如果m为位移单位的话那单位变换系数就应该为1了。

要注意单位变换系数的选取和自己设定的坐标单位有关。

而加速度数据的问题，如果数据为一维数据即全为振动加速度的值可以直接使用该程序，但如果像我一样是测量的原始结果，为振动加速度对应着时间需要对该程序进行修改，在读取数据的代码上修改为读取的x数组分离出加速度数据和时间数据。

分别代替原来的加速度x 和时间t。

我修改后的程序如下，之前除了数据部分直接写到程序里，不在需要输入了，而测量数据为时间对应着相应的加速度，m为位移单位m/s^2为加速度单位，单位变换系数为1，程序如下：
注意：如果直接复制程序的话注意中英文输入法切换导致程序字符有问题，自己在matlab 中改一下单引号为英文输入状态下的就好。

%频域积分
%%%%%%%%%%%%%%%%%%%%%%%%%%
clear; clc; close all hidden
%%%%%%%%%%%%%%%%%%
fni=input('频域积分-输入数据文件名：','s');
fid=fopen(fni,'r');
sf=12000%采样频率
fmin=0.1%最小截止频率
fmax=6000%最大截止频率
c=1%单位变换系数
it=2%积分次数
sx=’时间（s）’;%横向坐标轴的标注
sy1=’加速度（m/s^2）’;%纵向坐标轴输入单位的标注
sy2=’位移（m）’;%纵向坐标轴输出单位的标注
out.txt;%输出数据文件名
x=fscanf(fid,'%f',[1,inf]);%输入数据存成行向量
%加速度时间数据分离
for i=1:1:(length(x)/2)
%时间数据
t(i) = x(2*i-1);
%加速度数据
xx(i) = x(2*i);
end
status=fclose(fid);
n=length(xx);
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(xx,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*(0.5*sf-df);
%建立负的离散圆频率向量
w2=2*pi*(0.5*sf-df):-dw:0;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
%以积分次数为指数，建立圆频率变量向量
w=w.^it;
%进行积分的频域变换
a=zeros(1,nfft); a(2:nfft-1) =y(2:nfft-1)./w(2:nfft-1);
if it == 2
y=-a; %进行二次积分的相位变换
else
a1=imag(a); a2=real(a); y=a1-a2*i; %进行一次积分的相位变换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变换
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
subplot(2,1,1); plot(t,xx); xlabel(sx); ylabel(sy1); grid on; %绘制几分钱的时程曲线图形subplot(2,1,2); plot(t,y); xlabel(sx); ylabel(sy2); grid on; %绘制积分后的时程曲线图形%打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n, fprintf(fid,'%f \n',y(k)); end
status=fclose(fid);。