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

[nwnmrj]

那啥,小弟初学matlab,比较菜

有这样一组数据(样本见附件)
0.3987 1.3757 -5.0000 -0.8385
0.3984 1.3759 -5.0000 -0.8220
0.3989 1.3765 -5.0000 -0.7935
0.3986 1.3760 -5.0000 -0.7515
0.3984 1.3759 -5.0000 -0.7213
0.3976 1.3756 -5.0000 -0.7024
0.3973 1.3753 -5.0000 -0.6691
0.3978 1.3756 -5.0000 -0.6296
0.3969 1.3754 -5.0000 -0.5919
... .... ... ...
说明:
1.采样速率20K，k=1000
2.文件长度22K
3.第1，2 列激光位移，第3列是干扰,(忽略)第4列加速度
现在要求用频谱转换法,将加速度转为位移

算法是这个样子的:
如果在时间T内采集N个数据,

x(n)经离散傅里叶变换后得到的X(k)是一个长度为N的复数序列(频谱),

X(k)=DFT[x(n)]=[(a0，jb0)，(a1，jb1 )，...,(an-1 ，jbn-1 )]



x(n)中各谐波分量的幅值、圆频率及初相角可由式(1)~(3)求出

Ak=sqrt(ak^2+bk^2) (1)

ωk=2*pi*k/T (2)

φ=atan(bk/ak) (3)




根据信号叠加原理，任何周期信号是若干简谐信号的叠加。加速度信号可用式(4)表示，其相应的位移用式(5)表示，两者幅值和相位之间的关系如式(6)

a=Aa0cos(ω0t+φa0t)+Aa0cos(ω0t+φa0t)+...+Aan-1cos(ωn-1t+φan-1t) (4)

d=Ad0cos(ω0t+φd0t)+Aa0cos(ω0t+φd0t)+...+Aan-1cos(ωn-1t+φdn-1t) (5)

Adi=Aai/ωi^2 , φdi=φai-pi

下面是我自己用matlab写该算法的实现代码,比较臭
%读取数据
a='d:\1.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/20000;
%%%%%%%%%%%%%%%%%%
D=zeros(m,1);
D(i,1)=AK_D(i,1)*cos(OMEGA(i,1)*t+PHY_D(i,1));
for i=2:m
D(i,1)=D(i-1,1)+AK_D(i,1)*cos(OMEGA(i,1)*t+PHY_D(i,1));
end
plot(D)

现在的问题是:
1 通过该算法得到的位移D与数据里给出的两列位移都不匹配 (说实话我不太清楚两列位移各表示什么,但我算出来的D却和这两列数据毫无关系)
2.我的算法实现是不是存在问题,如果是,那我应该如何得到位移呢?

希望达人能帮帮忙,谢谢啦

[nwnmrj]

我突然发现我的帖子标题写错了,应该是加速度转位移
两次积分的方法我可以理解,不过老师要求的是谱转换法

[nwnmrj]

我发现我的算法写的有问题,重新修改以后
计算的结果仍然不正确
按照我学过的知识,最后计算出的位移波形应该是在采集时间T内
不通频率的波形叠加而成的,可我查看了一下各次谐波的幅值AK_D都十分小,叠加之后很不明显
希望高人能指点一下
我写的谱转换算法是否有问题

[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)

[linweixu168]

作者: 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...

留下联系方式 我也遇到了此问题 商量商量

[ljq897]

我的QQ:850171172,我也遇到同样问题，