由加速度时程到位移时程

wangbohust

本人已有编好的考虑地震动空间变异性的地震动加速度时程，但是要将其转为位移时程，这样通用程序使用才方便点，下面是变换思路：
1、对加速度时程进行傅立叶变换，得出其加速度的傅氏谱：
2、对要求的位移时程进行傅立叶变换，得出其位移的傅氏谱；
3、对位移傅氏谱进行二阶微分，得到加速度的傅氏谱值，利用已知的加速度谱和频率分量可以求出位移傅氏谱
4、对位移傅氏谱进行逆傅立叶变换求出位移时程。
具体推导方程和程序见附件，希望大家帮我解决一下这个问题，也希望大家讨论讨论该方法的可行性！谢谢！


   （附件已删除）

[ 本帖最后由 xinyuxf 于 2006-12-26 16:32 编辑 ]

toes

我只用时域方法做过，频域的方法没有试过，但是据说幅值需要概念清楚才能矫正。

其他的没有仔细考虑，请高手继续。

[ 本帖最后由 toes 于 2006-10-14 14:19 编辑 ]

wangbohust

toes，你好，时域的我也过了，但是积分出来的位移离零线偏离较远，资料上说是由于地震动加速度记录的误差分量造成的，即所谓的零点漂移现象，但我的波是自编的人工波，需要做滤波处理吗？你是怎么处理的，希望能和你交流，我的QQ:50235129

wangbohust

附上时程曲线，请高人分析原因，注：加速度为第一个附件中的加速度时程。

无水1324

建议楼上的看一下《MATLAB在振动信号处理中的应用》这本书
上面就有这个原例子
我没有电子版的，你可以去图书馆借或者买一本，才32块！

无水1324

[时间序列的积分变换]
实验室没有扫描的设备
我只好一个个地敲上去了，累！
要用的自己拿哈

clear;clc;close all;
%%%%%%%%%%%%%%%%%%%%%%%
fni=input('频域积分-输入数据文件名:','s');
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);           % 采样频率
fmin=fscanf(fid,'%f',1);         % 最小截止频率
fmax=fscanf(fid,'%d',1);         % 最大截止频率
c=fscanf(fid,'%d',1);            % 单位变换系数
it=fscanf(fid,'%s',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);
% 计算输入数据的长度
nt=length(x);
% 建立时间向量
t=0:1/sf:(n-1)/sf;
% 取大于n且与其最接近的2的整数次方为FFT长度
nf=2^nextpow2(nt);
% FFT 变换
y=fft(x,nfft);
% 计算频率间隔
df=sf/nfft;
% 计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
% 计算原频率间隔
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
    % 进行一次积分的相位变换
    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);

% ifft变换
y=ifft(a,nfft);
% 积分结果
y=real(y(1:n))*c;
subplot(211)
plot(t,x);
xlabel(sx);ylabel(sy1);
grid on;
% 积分后曲线
subplot(212)
plot(t,y);
xlabel(sx);ylabel(sy2);
grid on;
% 打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n
    fprint(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid)

toes

to wangbohust:
本来时域的程序是给别人做的，今天问别人要程序，结果他说我的是错的，只是一直没告诉我，晕死。

toes

无水1324 ：
不对啊，粗看了一下你贴的程序，至少这个地方不对：

    % 进行一次积分的相位变换
    real(y)=imag(a);
    imag(y)=-real(a);

这个赋值语句怎么对呢？

无水1324

原帖由 toes 于 2006-10-14 20:49 发表
无水1324 ：
不对啊，粗看了一下你贴的程序，至少这个地方不对：

    % 进行一次积分的相位变换
    real(y)=imag(a);
    imag(y)=-real(a);

这个赋值语句怎么对呢？

是对的，
保存的数据注意格式

toes

我编写的时域方法程序，方法很粗糙，只是尝试一下
信号太简单，可能需要加窗预处理一下；
积分应该用精度高一点的方法，比如辛普森积分；
去除直流的方法可以再考虑；

1. 
   
2. function Acc2Dis
   
3. % 加速度信号转换成位移信号
   
4. % 用时域法。
   
5. 
   
6. clear
   
7. clc
   
8. 
   
9. fs=1024;          %采样频率
   
10. N=20000;           %采样点数
    
11. t=(0:N-1)/fs;     %采样时间序列s
    
12. A = 2; % 幅值
    
13. 
    
14. w = 4; % 信号频率
    
15. a = A*sin(w*t); % 加速度信号
    
16. 
    
17. subplot(4,1,1);
    
18. plot(t,a);
    
19. title('加速度');
    
20. 
    
21. subplot(4,1,2);
    
22. 
    
23. % 积分。用了矩形面积计算积分，应该用精确一点的积分方法。
    
24. a_Mut_dt = a./fs;
    
25. for ii = 1:N-1
    
26.    s(ii) = sum( a_Mut_dt(1:ii+1) );
    
27. end
    
28. 
    
29. % 去除直流分量
    
30. s_fft = fft(s,N-1);
    
31. s_fft( 1:2 ) = 0;
    
32. s_fft( end-2: end) = 0;
    
33. s_ifft = real( ifft(s_fft,N-1) );
    
34. 
    
35. plot(t(2:end),s_ifft);
    
36. title('速度');
    
37. ylim( [-A/w A/w])
    
38. 
    
39. subplot(4,1,3);
    
40. 
    
41. % 积分。用了矩形面积计算积分，应该用精确一点的积分方法。
    
42. s_Mut_dt = s_ifft./fs;
    
43. for ii = 1:N-2
    
44.    d(ii) = sum( s_Mut_dt(1:ii+1) );
    
45. end
    
46. 
    
47. % 去除直流分量
    
48. d_fft = fft(d,N-2);
    
49. d_fft( 1:2 ) = 0;
    
50. d_fft( end-2: end) = 0;
    
51. d_ifft = real( ifft(d_fft,N-2) );
    
52. 
    
53. plot(t(3:end),d_ifft);
    
54. title('位移');
    
55. ylim( [-A/w^2 A/w^2])
    
56. 
    

复制代码

[ 本帖最后由 toes 于 2006-10-16 14:17 编辑 ]

无水1324

不错了
书上的的好啰嗦的
简单点，适用点，精度差不多就可以了

wangbohust

感谢以上两位高人的指点，让我受益非潜啊，那本书在我这个城市居然没有得卖，我只得网上定购了！
另外问toes两个问题：
"信号太简单，可能需要加窗预处理一下；
去除直流的方法可以再考虑；"

请问这两句话何解？我的程序何你差不多，只是没有你的去除直流方法这一点，能否详细解释一下，多谢。

toes

“信号太简单，可能需要加窗预处理一下；”
我的意思是我的这个例子分析的信号太简单，处理复杂信号时效果可能不好；你看，在位移信号的后面一部分幅值已经偏差比较大了。
采集的信号一般都要预处理一下的。我只是惯性思维，觉得加窗处理一下可能有助于减小截断误差，在这个程序中有没有必要我不知道，需要进一步尝试。

“去除直流的方法可以再考虑”
去除直流的方法比较多，可能用数字滤波器处理会更好；

你的程序没有去除直流那就会有信号的“漂移”，这是积分造成的，做出来当然不对。

sunworking

MATLAB在振动信号处理中的应用

nishishei

我用无水兄给的那个程序做过，原始信号是实验拿到的，定频激励下非线性减振器的响应，每组数据都有4000个周期以上，做出来的结果并不好，截断频率稍微做一点改动，结果就会有很大的变化，非常敏感。
后来就用时域方法，用了各种数值积分方法，结果都很差。在较少的周期数之内，出来的结果很好，但周期一多就不行了，得到的位移时程飘得一塌糊涂，不只是零偏和趋势项，高阶的大幅振荡也很明显。
为此很是头疼了一阵子，后来看到有信号放大器带硬件积分的，得到的位移时程还可以，目前还没试。
不知各位达人有什么高招？