加速度信号积分求速度信号

zhendongxue

加速度信号积分求速度信号
看了论坛里很多帖子，拿来用了但效果不好。z是我的数据，下面是我的程序，请大神指教。
%频域积分
%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc
fni=z;
sf =20000;     %采样频率
fmin =1;   %最小截止频率
fmax =19900;   %最大截止频率
c =9800;      %单位变换系数
it =1;     %积分次数
x=fni;%输入数据存成行向量  
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;                     
%进行积分的频域变换
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+2:nfft-ni+2)=y(nfft-na+2:nfft-ni+2);                                 
%IFFT变换                                         
y=ifft(a,nfft);               
%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;            
%绘制积分前的时程曲线图形                        
subplot(2,1,1);                                                                 
plot(t,x);      
%添加横向坐标轴的标注
xlabel('sx');           
%添加纵向坐标轴的标注   
ylabel('sy1');               
grid on;                     
%绘制积分后的时程曲线图形              
subplot(2,1,2);               
plot(t,y);
%添加横向坐标轴的标注
xlabel('sx');                  
%添加纵向坐标轴的标注   
ylabel('sy2');                  
grid on;