|
楼主 |
发表于 2008-7-30 21:31
|
显示全部楼层
请问这个信号通过短时傅立叶变换后会有突变么?
非常感谢songzy41!!您回复的帖子我看了非常感激!!!!!!you are a good man!!!
非常抱歉,上面我发的帖子没有写完全,点错了,直接就发出去了,请大家原谅我的马虎失误,还请大家帮忙看看我的问题,多谢了!!
这个信号是G1,是由以下程序得到的:
M=[2573*10^3 0 0;0 2558*10^3 0;0 0 562*10^3];
>> k1=554*10^3;
>> k2=926*10^3;
>> k3=836*10^3;
>> K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
>> M1=inv(M);
>> Q=K*M1;
>> W=sqrt(eig(Q));
>> w1=W(1,1);w2=W(2,1);
>> xishu=2*0.05/(w1+w2)*[w1*w2;1];
>> a=xishu(1,1);b=xishu(2,1);
>> C=a*M+b*K;
>> U=[0;0;0];U1=[0;0;0];U2=[0;0;0];
t=0;
g=[0;0;0];
dt=0.02;sei=1.4;tt=sei*dt;
i=1;
G1=zeros(800,1); G2=zeros(800,1); G3=zeros(800,1);
G=zeros(801,1);
while t<8
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
i=i+1;
t=t+dt;
U2=-M1*(M*g+C*U1+K*U);
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
end
k1=398*10^3;
K=[k1+k2 -k2 0;-k2 k2+k3 -k3;0 -k3 k3];
Q=K*M1;
W=sqrt(eig(Q));
w1=W(1,1);w2=W(2,1);
xishu=2*0.05/(w1+w2)*[w1*w2;1];
a=xishu(1,1);b=xishu(2,1);
C=a*M+b*K;
while t<=16
dgtt=(sin(t+tt)-sin(t))*[1;1;1];
dptt=-M*dgtt;
dgdt=(sin(t+dt)-sin(t))*[1;1;1];
dpdt=-M*dgdt;
K1=K+6*M/(tt^2)+3*C/tt;
P1=dptt+M*(6*U1/tt+3*U2)+C*(3*U1+0.5*U2*tt);
dutt=inv(K1)*P1;
du2tt=6*(dutt-U1*tt-0.5*U2*tt^2)/tt^2;
du1tt=3*(dutt-U1*tt-(U2*tt^2)/6)/tt;
du2dt=du2tt/sei;
du1dt=U2*dt+0.5*du2dt*dt;
dudt=U1*dt+0.5*U2*dt^2+(du2dt*dt^2)/6;
U=U+dudt;
U1=U1+du1dt;
g=g+dgdt;
G(i,1)=t;
t=t+dt;
G1(i,1)=U2(1,1);
G2(i,1)=U2(2,1);
G3(i,1)=U2(3,1);
i=i+1;
U2=-M1*(M*g+C*U1+K*U);
end
plot(G,G1)就可以得到这个信号的图形,这里我没法把图给传过来,麻烦大侠们运行一下我的程序,就可以得到G1的图形了,现在我用短时傅立叶变换对这个信号分解,程序如下:
X=G1;
T=1:length(X);
N=length(X);
H=20;
TRACE=0;
[TFR,T,F]=TFRSTFT(X,T,N,H,TRACE);
mesh(T,F,TFR)
但看不到有任何突变,我的目的是证明这个信号在8秒时有明显突变,
是不是短时傅立叶变换做不到这个啊?还请您多费心帮忙看看好么?谢谢!! |
|