测得的加速度信号求速度信号时出错!请各位大侠帮忙

hgdzxz

以下是时域积分法程序！
sf=10240;  %采样频率
n=length(x);  %信号长度
t=0:1/sf:(n-1)/sf;  %采用时间
t1=1/sf;
%采用辛布森时域积分法求速度
yvs(1)=t1*(x(1)+x(2))/2;
for k=2:n-1
yvs(k)=yvs(k-1)+t1*(x(k-1)+4*x(k)+x(k+1))/6;
end
yvs(n)=yvs(n-1);
%画图
plot(t,yvs);


以下是时域积分法程序
sf=10240;%采样频率
fmin=10;%最小截止频率
fmax=1500;%最大截止频率
it=1;%积分次数
n=length(x);
%建立时间向量
t=0:1/sf:(n-1)/sf;
%大于并最接近n的2的幂次方为FFT长度
nfft=2^nextpow2(n);
%FFT变换
y=fft(x,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*(0.5*sf);
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%将正负圆频率向量组合成一个向量
w=[w1,w2];
%以积分次数为指数，建立圆频率变量向量
w=w.^it;
ww=w';
%进行积分的频域变换
a=zeros(1,nfft);
a(2:nfft-1) =y(2:nfft-1)./ww(2:nfft-1);
%进行一次积分的相位变换
a1=imag(a);
a2=real(a);
y=a1-a2*i;
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);
%画图
plot(t,y);

问题：为什么时域积分法和频域积分法求出来的差别这么大呢？
时域积分法求出的速度随时间向下倾斜？为什么？
频域积分法求出的速度似乎是两个信号，是与我测量的数据有关吗？还是其他原因？

clownzhi

我自己的理解
积分上有时域频域上
一般时域上直接的累加会存在漂移，也就是累计误差
频域上积分的话，由于信号不是严格的整周期采样，积分变化中会存在频率分散和超低频部分
建议是做积分前先做高通滤波处理——积分
一般一次积分还是比较准确的
二次的话重复

hyl2323

1）时域方法，你没有在积分前去零均值化。
2）频域方法没你写的这么复杂，怀疑你写错了。

songzy41

时域积分正如教研室主任指出的没有消除直流分量。而频域积分中LZ是引用王济一书中的程序，在本论坛中已指出该程序有错：
http://forum.vibunion.com/thread-47704-1-1.html
针对LZ的数据修改程序如下，数据在x中：
x0=mean(x)
x=(x-x0)';
sf=10240;  %采样频率
t1=1/sf;
%%%%%辛普森（simpson）算法时域积分求速度
yvs(1)=t1*(x(1)+x(2))/2;
n=length(x);
t=(1:n)/sf;
for k=2:n-1
yvs(k)=yvs(k-1)+t1*(x(k-1)+4*x(k)+x(k+1))/6;
end
yvs(n)=yvs(n-1);
%%%%%辛普森（simpson）算法时域积分求速度
figure(1)
subplot 211
title('实际信号输入信号')
plot(t,x); grid;
subplot 212
plot(t,yvs,'r'); grid;
title('时域积分信号')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfft=n;
c=1;
fmin=10;%最小截止频率
fmax=1500;%最大截止频率
it=1;%积分次数
%FFT变换
y=fft(x,nfft);
%计算频率间隔（Hz/s）
df=sf/nfft;
%计算指定频带对应频率数组的下标
ni=round(fmin/df+1);
na=round(fmax/df+1);
%计算圆频率间隔（rad/s）
dw=2*pi*df;
%建立正的离散圆频率向量
w1=0:dw:2*pi*0.5*sf;
%建立负的离散圆频率向量
w2=-2*pi*(0.5*sf-df):dw:-dw;
%将正负圆频率向量组合成一个向量
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
%进行一次积分的相位变换
a1=imag(a);
a2=real(a);
y=a1-a2*i;
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);
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;
%绘制积分前的时程曲线图形
figure
subplot 211;
plot(t,x);
title('实际时程信号');
grid on;
%绘制积分后的时程曲线图形
subplot 212;
plot(t,y);
hold on
plot(t,yvs,'r')
legend('频域积分','时域积分',2);
title('频域积分、时域积分信号')
grid on;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
得图如下，频域积分和时域积分比较接近。

李元

x0=mean(x)
x=(x-x0)';
sf=10240;  %采样频率
t1=1/sf;
%%%%%辛普森（simpson）算法时域积分求速度
yvs(1)=t1*(x(1)+x(2))/2;
n=length(x);
t=(1:n)/sf;
for k=2:n-1
yvs(k)=yvs(k-1)+t1*(x(k-1)+4*x(k)+x(k+1))/6;
end
yvs(n)=yvs(n-1);
%%%%%辛普森（simpson）算法时域积分求速度
figure(1)
subplot 211
title('实际信号输入信号')
plot(t,x); grid;
subplot 212
plot(t,yvs,'r'); grid;
title('时域积分信号')
%这里需要加这个吗？
meanYvs=mean(yvs);
yvs=yvs-meanYvs;
%%%%%辛普森（simpson）算法时域积分求位移
yds(1)=t1*(yvs(1)+yvs(2))/2;
for k=2:n-1
yds(k)=yds(k-1)+t1*(yvs(k-1)+4*yvs(k)+yvs(k+1))/6;
end
yds(n)=yds(n-1);
figure(2)
plot(t,yds);grid;
title('位移信号')

songzy41大侠，以上是我的二次积分程序，二次积分后的图为
file:///C:/DOCUME~1/ADMINI~1/LOCALS~1/Temp/NQ%7BAY4W1L1KU4)OT9$A%7DD(1.jpg

请问如何矫正这个图呢？请指明方向，谢谢！ 等待…

李元

老师您好，我想问您savitzky-golay七点三次平滑器的代码该是如何的？
h2=[-2,3,6,7,6,3,-2]';
%以下三句如何修改？y原始信号需要进行平滑处理，y2是经过平滑处理后的信号
y2(1)=y(1);
y2(2)=y(2);
y2(3)=y(3);

for i=4:(n-3)
yDTemp=[y(i-3),y(i-2),y(i-1),y(i),y(i+1),y(i+2),y(i+3)];
trendTemp=(yDTemp*h2)/21;
y2(i)=y(i)-trendTemp;
end;
%同上，以下三句如何修改？
y2(n-2)=y(n-2)-trendTemp;
y2(n-1)=y(n-1)-trendTemp;
y2(n)=y(n)-trendTemp;
谢谢您，不甚感激！

李元

经过学习，在MATLAB中的sgolay库函数正是该平滑器的系数生成器；help sgolay，可以知道该矩阵前几行和后几行分别对应平滑数据的开头和结尾的数据平滑，中间的一行对应中间数据的平滑。希望对向我这样的初学者有用。

李元

大侠，请问有没有比这个辛普森积分法更加精确的方法来进行时域积分呢？有没有改进型的积分方法呢？请指点。

165390292

很有用，学习一下。