回复 #46 zhlong 的帖子
是采样频率的问题，我的0.005，你的0.001，我做出来了，但现在的问题还是，为什么在5秒出的损伤却出现在5秒之前的那个位置呢？[ 本帖最后由 zhlong 于 2007-7-10 13:35 编辑 ]
回复 #47 猴王 的帖子
换别的小波试试看 现在这个答案应该比较完整了，不用ode45，改用定步长的4阶龙格库塔法积分，即下面贴的ode4。function Y = ode4(odefun,tspan,y0,varargin)
%ODE4 Solve differential equations with a non-adaptive method of order 4.
% Y = ODE4(ODEFUN,TSPAN,Y0) with TSPAN = integrates
% the system of differential equations y' = f(t,y) by stepping from T0 to
% T1 to TN. Function ODEFUN(T,Y) must return f(t,y) in a column vector.
% The vector Y0 is the initial conditions at T0. Each row in the solution
% array Y corresponds to a time specified in TSPAN.
%
% Y = ODE4(ODEFUN,TSPAN,Y0,P1,P2...) passes the additional parameters
% P1,P2... to the derivative function as ODEFUN(T,Y,P1,P2...).
%
% This is a non-adaptive solver. The step sequence is determined by TSPAN
% but the derivative function ODEFUN is evaluated multiple times per step.
% The solver implements the classical Runge-Kutta method of order 4.
%
% Example
% tspan = 0:0.1:20;
% y = ode4(@vdp1,tspan,);
% plot(tspan,y(:,1));
% solves the system y' = vdp1(t,y) with a constant step size of 0.1,
% and plots the first component of the solution.
%
if ~isnumeric(tspan)
error('TSPAN should be a vector of integration steps.');
end
if ~isnumeric(y0)
error('Y0 should be a vector of initial conditions.');
end
h = diff(tspan);
if any(sign(h(1))*h <= 0)
error('Entries of TSPAN are not in order.')
end
try
f0 = feval(odefun,tspan(1),y0,varargin{:});
catch
msg = ['Unable to evaluate the ODEFUN at t0,y0. ',lasterr];
error(msg);
end
y0 = y0(:); % Make a column vector.
if ~isequal(size(y0),size(f0))
error('Inconsistent sizes of Y0 and f(t0,y0).');
end
neq = length(y0);
N = length(tspan);
Y = zeros(neq,N);
F = zeros(neq,4);
Y(:,1) = y0;
for i = 2:N
ti = tspan(i-1);
hi = h(i-1);
yi = Y(:,i-1);
F(:,1) = feval(odefun,ti,yi,varargin{:});
F(:,2) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,1),varargin{:});
F(:,3) = feval(odefun,ti+0.5*hi,yi+0.5*hi*F(:,2),varargin{:});
F(:,4) = feval(odefun,tspan(i),yi+hi*F(:,3),varargin{:});
Y(:,i) = yi + (hi/6)*(F(:,1) + 2*F(:,2) + 2*F(:,3) + F(:,4));
end
Y = Y.';
结构的动力学方程函数
function xdot=danziyoudukuangjia(t,x)
M=31.2;
if(t<5)
K=67.51*10^5;
else
K=67.51*10^5*0.95;
end
C=1.4*10^3;
xdot=;
信噪比计算
ss=sum(y1.^2);
nn=sum((y2).^2);
snr=10*log10((ss/nn))
主程序
t0=0.001;
te=10;
dt=0.001;%采样间隔
t=';
x0=;
=ode4('danziyoudukuangjia',t,x0);
d=x(:,1);
v=x(:,2);
s1=(diff(v(1:end-1))+diff(v(2:end)))/2;%中心差分法
s1=; %少的两个点用前向和后向差分法补齐
s1=s1/dt; %近似为导数
figure;set(gcf,'Color','W')
subplot(311);plot(t,d); title('位移')
subplot(312);plot(t,v); title('速度')
subplot(313);plot(t,s1); title('加速度')
y1=s1(150:end);%去掉信号前面一部分的瞬态成份
y2=0.005*randn(size(y1));%添加噪声
s=y1+y2;
%下面进行离散的单尺度小波变换并生成，各尺度上的信号
=wavedec(s,5,'bior6.8');%对第一信号进行3尺度一维离散小波分解，采用墨西哥小帽函数
%提取结构的低频和高频信号
ca3=appcoef(c,l,'bior6.8',5);%提取第三尺度系数的低频
=detcoef(c,l,);%提取第一、二、三尺度系数的高频
%重构信号的低频和高频部分
a3=wrcoef('a',c,l,'bior6.8',3);
d1=wrcoef('d',c,l,'bior6.8',1);
d2=wrcoef('d',c,l,'bior6.8',2);
d3=wrcoef('d',c,l,'bior6.8',3);
d4=wrcoef('d',c,l,'bior6.8',4);
d5=wrcoef('d',c,l,'bior6.8',5);
%显示多尺度一维信号的分解结果
figure;set(gcf,'Color','W')
subplot(611);plot(t(150:end),a3);title('信号的低频细节部分');%AXIS();
subplot(612);plot(t(150:end),d1);title('Detailsd1');%AXIS();
subplot(613);plot(t(150:end),d2);title('Detailsd2');%AXIS();
subplot(614);plot(t(150:end),d3);title('Detailsd3');%AXIS();
subplot(615);plot(t(150:end),d4);title('Detailsd4');%AXIS();
subplot(616);plot(t(150:end),d5);title('Detailsd5');%AXIS();
这是得到的位移、速度、加速度：
这是噪声为0.005*randn(size(y1))，信噪比约13db时的小波分解图和EMD分解结果：
这是完全无噪声的分解情况：
[ 本帖最后由 zhlong 于 2007-7-10 14:57 编辑 ]
回复 #50 猴王 的帖子
小波函数没改，只是将变步长积分改为了定步长积分。不用客气，欢迎常来论坛。
[ 本帖最后由 zhlong 于 2007-7-10 16:12 编辑 ] 这个积分程序不错，收藏了。
回复 #51 zhlong 的帖子
请问，为什么我把损伤时刻改到30s的时候，损伤就是别出不来了呢，我只是把损伤的时刻，和采样频率改变了，损伤就看不出来了，而且还出现在5s时刻，这是不是也是因为m,k,c这几个系数阿function xdot=danziyoudukuangjia(t,x)
F=20*sin(25*t+14);
M=31.2;
if(t<30)
K=67.51*10^5;
else
K=67.51*10^5*0.95;
end
C=1.4*10^3;
xdot=;
close all;
t0=0.001;
te=60;
dt=0.005;%采样间隔
t=';
x0=;
=ode4('danziyoudukuangjia',t,x0);
d=x(:,1);
v=x(:,2);
s1=(diff(v(1:end-1))+diff(v(2:end)))/2;%中心差分法
s1=; %少的两个点用前向和后向差分法补齐
s1=s1/dt; %近似为导数
figure;set(gcf,'Color','W')
subplot(311);plot(t,d); title('位移')
subplot(312);plot(t,v); title('速度')
subplot(313);plot(t,s1); title('加速度')
y1=s1(150:end);%去掉信号前面一部分的瞬态成份
y2=0.005*randn(size(y1));%添加噪声
s=y1+y2;
ss=sum(y1.^2);
nn=sum((y2).^2);
snr=10*log10((ss/nn))
%下面进行离散的单尺度小波变换并生成，各尺度上的信号
=wavedec(s,5,'bior6.8');%对第一信号进行3尺度一维离散小波分解，采用墨西哥小帽函数
%提取结构的低频和高频信号
ca3=appcoef(c,l,'bior6.8',5);%提取第三尺度系数的低频
=detcoef(c,l,);%提取第一、二、三尺度系数的高频
%重构信号的低频和高频部分
a3=wrcoef('a',c,l,'bior6.8',3);
d1=wrcoef('d',c,l,'bior6.8',1);
d2=wrcoef('d',c,l,'bior6.8',2);
d3=wrcoef('d',c,l,'bior6.8',3);
d4=wrcoef('d',c,l,'bior6.8',4);
d5=wrcoef('d',c,l,'bior6.8',5);
%显示多尺度一维信号的分解结果
figure;set(gcf,'Color','W')
subplot(611);plot(t(150:end),a3);title('信号的低频细节部分');%AXIS();
subplot(612);plot(t(150:end),d1);title('Detailsd1');%AXIS();
subplot(613);plot(t(150:end),d2);title('Detailsd2');%AXIS();
subplot(614);plot(t(150:end),d3);title('Detailsd3');%AXIS();
subplot(615);plot(t(150:end),d4);title('Detailsd4');%AXIS();
subplot(616);plot(t(150:end),d5);title('Detailsd5');%AXIS();