由加速度时程到位移时程
本人已有编好的考虑地震动空间变异性的地震动加速度时程，但是要将其转为位移时程，这样通用程序使用才方便点，下面是变换思路：1、对加速度时程进行傅立叶变换，得出其加速度的傅氏谱：
2、对要求的位移时程进行傅立叶变换，得出其位移的傅氏谱；
3、对位移傅氏谱进行二阶微分，得到加速度的傅氏谱值，利用已知的加速度谱和频率分量可以求出位移傅氏谱
4、对位移傅氏谱进行逆傅立叶变换求出位移时程。
具体推导方程和程序见附件，希望大家帮我解决一下这个问题，也希望大家讨论讨论该方法的可行性！谢谢！
（附件已删除）
[ 本帖最后由 xinyuxf 于 2006-12-26 16:32 编辑 ] 我只用时域方法做过，频域的方法没有试过，但是据说幅值需要概念清楚才能矫正。
其他的没有仔细考虑，请高手继续。
[ 本帖最后由 toes 于 2006-10-14 14:19 编辑 ] toes，你好，时域的我也过了，但是积分出来的位移离零线偏离较远，资料上说是由于地震动加速度记录的误差分量造成的，即所谓的零点漂移现象，但我的波是自编的人工波，需要做滤波处理吗？你是怎么处理的，希望能和你交流，我的QQ:50235129 附上时程曲线，请高人分析原因，注：加速度为第一个附件中的加速度时程。 建议楼上的看一下《MATLAB在振动信号处理中的应用》这本书
上面就有这个原例子
我没有电子版的，你可以去图书馆借或者买一本，才32块！ [时间序列的积分变换]
实验室没有扫描的设备
我只好一个个地敲上去了，累！
要用的自己拿哈
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',); % 按行读入原始信号数据
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=;
%以积分次数为指数，建立元频率向量
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) to wangbohust:
本来时域的程序是给别人做的，今天问别人要程序，结果他说我的是错的，只是一直没告诉我，晕死。 无水1324 ：
不对啊，粗看了一下你贴的程序，至少这个地方不对：
% 进行一次积分的相位变换
real(y)=imag(a);
imag(y)=-real(a);
这个赋值语句怎么对呢？ 原帖由 toes 于 2006-10-14 20:49 发表
无水1324 ：
不对啊，粗看了一下你贴的程序，至少这个地方不对：
% 进行一次积分的相位变换
real(y)=imag(a);
imag(y)=-real(a);
这个赋值语句怎么对呢？
是对的，
保存的数据注意格式 我编写的时域方法程序，方法很粗糙，只是尝试一下
信号太简单，可能需要加窗预处理一下；
积分应该用精度高一点的方法，比如辛普森积分；
去除直流的方法可以再考虑；
function Acc2Dis
% 加速度信号转换成位移信号
% 用时域法。
clear
clc
fs=1024; %采样频率
N=20000; %采样点数
t=(0:N-1)/fs; %采样时间序列s
A = 2; % 幅值
w = 4; % 信号频率
a = A*sin(w*t); % 加速度信号
subplot(4,1,1);
plot(t,a);
title('加速度');
subplot(4,1,2);
% 积分。用了矩形面积计算积分，应该用精确一点的积分方法。
a_Mut_dt = a./fs;
for ii = 1:N-1
s(ii) = sum( a_Mut_dt(1:ii+1) );
end
% 去除直流分量
s_fft = fft(s,N-1);
s_fft( 1:2 ) = 0;
s_fft( end-2: end) = 0;
s_ifft = real( ifft(s_fft,N-1) );
plot(t(2:end),s_ifft);
title('速度');
ylim( [-A/w A/w])
subplot(4,1,3);
% 积分。用了矩形面积计算积分，应该用精确一点的积分方法。
s_Mut_dt = s_ifft./fs;
for ii = 1:N-2
d(ii) = sum( s_Mut_dt(1:ii+1) );
end
% 去除直流分量
d_fft = fft(d,N-2);
d_fft( 1:2 ) = 0;
d_fft( end-2: end) = 0;
d_ifft = real( ifft(d_fft,N-2) );
plot(t(3:end),d_ifft);
title('位移');
ylim( [-A/w^2 A/w^2])
[ 本帖最后由 toes 于 2006-10-16 14:17 编辑 ] 不错了
书上的的好啰嗦的
简单点，适用点，精度差不多就可以了 感谢以上两位高人的指点，让我受益非潜啊，那本书在我这个城市居然没有得卖，我只得网上定购了！
另外问toes两个问题：
"信号太简单，可能需要加窗预处理一下；
去除直流的方法可以再考虑；"
请问这两句话何解？我的程序何你差不多，只是没有你的去除直流方法这一点，能否详细解释一下，多谢。 “信号太简单，可能需要加窗预处理一下；”
我的意思是我的这个例子分析的信号太简单，处理复杂信号时效果可能不好；你看，在位移信号的后面一部分幅值已经偏差比较大了。
采集的信号一般都要预处理一下的。我只是惯性思维，觉得加窗处理一下可能有助于减小截断误差，在这个程序中有没有必要我不知道，需要进一步尝试。
“去除直流的方法可以再考虑”
去除直流的方法比较多，可能用数字滤波器处理会更好；
你的程序没有去除直流那就会有信号的“漂移”，这是积分造成的，做出来当然不对。 MATLAB在振动信号处理中的应用 我用无水兄给的那个程序做过，原始信号是实验拿到的，定频激励下非线性减振器的响应，每组数据都有4000个周期以上，做出来的结果并不好，截断频率稍微做一点改动，结果就会有很大的变化，非常敏感。
后来就用时域方法，用了各种数值积分方法，结果都很差。在较少的周期数之内，出来的结果很好，但周期一多就不行了，得到的位移时程飘得一塌糊涂，不只是零偏和趋势项，高阶的大幅振荡也很明显。
为此很是头疼了一阵子，后来看到有信号放大器带硬件积分的，得到的位移时程还可以，目前还没试。
不知各位达人有什么高招？