对振动信号的积分处理--程序运行问题

wm123

需要对一个txt文件内的路面激励振动信号（加速度）进行一次、二次积分，使用了王济老师的书上的程序如下，但是不能出结果，看到论坛的一些帖子，使用论坛的程序，但是还是运行出错，结果不对，请各位专家们指导一下：该对程序进行如何修改之后，再如何操作才能得到结果？因为刚刚接触这些，还望各位能多多指导，不胜感激！！！由于数据有570K上传不了，所以将其复制了前158k上传了，不知道能否使用。

书上的程序如下：
%频率积分
%%%%%%%
clear
clc
close all hidden
%%%%%%%%%
fni=input('频率积分-输入数据文件名;','s');
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);%采样频率值
fmin=fscanf(fid,'%f',1);
fmax=fscanf(fid,'%f',1);
c=fscanf(fid,'%f',1);
it=fscanf(fid,'%f',2);
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);
%取信号长度
n=length(x);
t=0:1/sf:(n-1)/sf
nfft=2^nextpow2(n);
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);
y=ifft(a,nfft);
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;
fid=fopen(fno,'w');
for k=1:n
   fprintf(fid,'%f%f\n',t(k),y(k));
end
status=fclose(fid);



在论坛上找到的程序如下：
%频域积分%%%%%%%%%%%%%%%%%%%%%%%%%%
clear;
clc;
close all hidden
%%%%%%%%%%%%%%%%%%
fni=input('频域积分-输入数据文件名：','s');
fid=fopen(fni,'r');
sf=fscanf(fid,'%f',1);%采样频率
fmin=fscanf(fid,'%f',1);%最小截止频率
fmax=fscanf(fid,'%f',1);%最大截止频率
c=fscanf(fid,'%f',1);%单位变换系数
it=fscanf(fid,'%f',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);
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+1:nfft-ni+1)=y(nfft-na+1:nfft-ni+1);
y=ifft(a,nfft); %IFFT变换%取逆变换的实部n个元素并乘以单位变换系数为积分结果
y=real(y(1:n))*c;subplot(2,1,1);
plot(t,x);
xlabel(sx);
ylabel(sy1);
grid on; %绘制积分前的时程曲线图形

%打开文件输出积分后的数据
fid=fopen(fno,'w');
for k=1:n, fprintf(fid,'%f \n',y(k));
end
status=fclose(fid);

wm123