zhengjirui
2012-07-11, 16:48
sf=5000;
fno='4.doc';
x=load('2.txt','%d');
%x=y;
% function h=f(x,y);
% h=x;
%status=fclose(fid);
syms s;
f=s;
m=s^3/6;
N=length(x);
z=zeros(N,1);
t=(0:1/sf:(N-1)/sf)';
% x(1)=x(0);
y(1)=0;
for n=1:N-1;
% x(n+1)=x(n)+1/sf;
% z0=y(n)+(1/sf)*feval(h,x(n),y(n));
% y(n+1)=y(n)+1/sf/2*(feval(h,x(n),y(n))+feval(fun,x(n+1),z0));
% z0=y(n)+(1/sf)*subs(sym(f),findsym(sym(f)),x(n));
y(n+1)=y(n)+(1/sf)/2*(subs(sym(f),findsym(sym(f)),x(n))+subs(sym(f),findsym(sym(f)),x(n+1)));
z(n+1)=z(n)+(1/sf)/2*(subs(sym(m),findsym(sym(m)),x(n))+subs(sym(m),findsym(sym(m)),x(n+1)));
end
subplot(3,1,1);
plot(t,x);
grid on;
subplot(3,1,2);
plot(t,y);
subplot(3,1,3);
plot(t,z);
2.txt 就是加速度的信号 现在用的是sin（pi*t）产生的数据，但是结果不对
fno='4.doc';
x=load('2.txt','%d');
%x=y;
% function h=f(x,y);
% h=x;
%status=fclose(fid);
syms s;
f=s;
m=s^3/6;
N=length(x);
z=zeros(N,1);
t=(0:1/sf:(N-1)/sf)';
% x(1)=x(0);
y(1)=0;
for n=1:N-1;
% x(n+1)=x(n)+1/sf;
% z0=y(n)+(1/sf)*feval(h,x(n),y(n));
% y(n+1)=y(n)+1/sf/2*(feval(h,x(n),y(n))+feval(fun,x(n+1),z0));
% z0=y(n)+(1/sf)*subs(sym(f),findsym(sym(f)),x(n));
y(n+1)=y(n)+(1/sf)/2*(subs(sym(f),findsym(sym(f)),x(n))+subs(sym(f),findsym(sym(f)),x(n+1)));
z(n+1)=z(n)+(1/sf)/2*(subs(sym(m),findsym(sym(m)),x(n))+subs(sym(m),findsym(sym(m)),x(n+1)));
end
subplot(3,1,1);
plot(t,x);
grid on;
subplot(3,1,2);
plot(t,y);
subplot(3,1,3);
plot(t,z);
2.txt 就是加速度的信号 现在用的是sin（pi*t）产生的数据，但是结果不对