位移转加速度(谱转换法)

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列加速度
现在要求用频谱转换法,将加速度转为位移

(算法描述见附件)
下面是我自己用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 于 2008-1-27 00:05 编辑 ]

songzy41

可参看王济和胡晓编的 “MATLAB在振动信号处理中的应用”一书的笫5章中，就有在频率域把加速度积分转成位移，又有MATLAB程序。同时可参看
http://forum.vibunion.com/forum/thread-47704-1-1.html

wanyeqing2003

位移转加速度应该是两次微分。

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)

songzy41

1，在一楼中给出的数据1.txt，这次6楼读的数据是aa.txt，采样频率变了？一楼中给出的采样频率为20000，而这次T=1。因为要计算频率，采样频率是一个重要参数。
2，OMEGA的计算错了，OMEGA（1）=0，OMEGA（502）~OMEGA（1000）都是负频率。
3，建议看一下王济和胡晓编的 “MATLAB在振动信号处理中的应用”一书和帖子
http://forum.vibunion.com/forum/thread-47704-1-1.html