【求助】位移转加速度(谱转换法)

[nwnmrj]

%修改后的算法
%读取数据
a='d:\aa.txt';
data=load(a);
[m,n]=size(data);

%提取加速度,
acce=data(:,4);
XK=fft(acce);

re=real(XK);
im=imag(XK);
%AK=zeros(m,1);

T=1;
for i=1:m
%计算幅值
AK_A(i,1)=sqrt(re(i,1)^2+im(i,1)^2);
%计算初相
PHY_A(i,1)=atan(im(i,1)/re(i,1));
%计算角频率
OMEGA(i,1)=2*pi*i/T;
end

for i=1:m
%计算幅值
AK_D(i,1)=AK_A(i,1)/(OMEGA(i,1)^2);
%计算初相
PHY_D(i,1)=PHY_A(i,1)-pi;
end

%计算位移
%%%%%%%%%%%%%%%%%%
t=1/1000;
step=t;
%%%%%%%%%%%%%%%%%%
D=zeros(m,1);
i=1;
for tt=t:step:1
for k=1:m
D(i,1)=D(i,1)+AK_D(k,1)*cos(OMEGA(k,1)*tt+PHY_D(k,1));
end
i=i+1;
end
plot(D)