examp.m

来自「这是在网上下的一个东东」· M 代码 · 共 176 行

clear;
noise=0.05;
N=211;
M=211;
dt=.5;
dF=1/dt;
t=linspace(-5,100,N);
F=linspace(0,dF/2,N-1);
%instrument response is a critically-damped pulse
sigi=10;
for i=1:N-1
if (t(i)<0)
g(i) = 0;
else
g(i) = t(i)*exp(-t(i)/sigi);
end
end
gmax=max(g);
g=g'/gmax;

%true signal is two pulses of sig standard deviation
sig=2;
mtrue = (exp(-((t(1:N-1)-8).^2/(sig^2*2)))+0.5*exp(-((t(1:N-1)-25).^2/(sig^2*2))))';
mtrue=mtrue/max(mtrue);
d=conv(g,mtrue);
nn=noise*randn(length(d),1);
dn=d+nn;

figure(1)
bookfonts
plot(t(1:N-1),mtrue,'k')
xlabel('Time (s)')
ylabel('Acceleration (m/s^2)')
axis tight
%print -deps truemodel.eps

figure(2)
bookfonts
plot(t(1:N-1),g,'k')
xlabel('Time (s)')
ylabel('V')
axis tight
%print -deps impresp.eps

figure(3)
bookfonts
plot(t(1:N-1),d(1:N-1),'k')
xlabel('Time (s)')
ylabel('V')
axis tight
%print -deps dsynfig.eps

%generate spectra
%%pad all time series to 2*(N-1) points to avoid wrap-around effects
dspec=fft([d; 0]);
dnspec=fft([dn; 0]);
gspec=fft([g; zeros(length(g),1)]);
mperf=real(ifft(dspec./gspec));
mn=real(ifft(dnspec./gspec));

figure(4)
bookfonts
plot(t(1:N-1),mperf(1:N-1),'k')
xlabel('Time (s)')
ylabel('Acceleration (m/s^2)')
axis tight
%print -deps perfectsol.eps

figure(5)
bookfonts
plot(t(1:N-1),dn(1:N-1),'k')
xlabel('Time (s)')
ylabel('V')
axis tight
%print -deps dnoisefig.eps

figure(6)
bookfonts
plot(t(1:N-1),mn(1:N-1),'k')
xlabel('Time (s)')
ylabel('Acceleration (m/s^2)')
axis tight
%print -deps m2noise.eps

gs=abs(gspec(1:length(gspec)/2));
ds=abs(dspec(1:length(dspec)/2));
dns=abs(dnspec(1:length(dspec)/2));

figure(7)
bookfonts
loglog(F,gs,'k');
xlabel('f (Hz)')
ylabel('Spectral Amplitude')
ylim([1e-3 1e2]);
xlim([1/N 1]);
hold on
loglog(F,ds,'k-.')
loglog(F,dns,'k-')
hold off
%print -deps spectra1.eps

figure(8)
bookfonts
loglog(F,ds./gs,'k')
loglog(F,dns./gs,'k-')
xlabel('f (Hz)')
ylabel('Spectral Amplitude')
axis tight
%print -deps spectra2.eps

wcount=0;
for wlevexp=-1:.1:2
  wcount=wcount+1;
  wlev=10^wlevexp
%
%apply water level damping
%
  for i=1:length(gspec)
    if abs(gspec(i)) < wlev
      gwspec(i,1)=wlev*(gspec(i)./abs(gspec(i)));
    else
      gwspec(i,1)=gspec(i);
    end
  end

  mw=real(ifft(dnspec./gwspec));
  wval(wcount)=wlev;
  mnorm(wcount)=norm(mw);
  resid(wcount)=norm((dnspec./gwspec).*gspec-dnspec);
  mstore(:,wcount)=mw;
%
%animation of the evolving solution (not saved)
%
  figure(9)
  bookfonts
  plot(t(1:N-1),mw(1:N-1),'k')
  xlabel('Time (s)')
  ylabel('Acceleration (m/s^2)')
  drawnow;
  pause(0.1);
end

figure(10)
bookfonts
plot(resid,mnorm,'ko')
xlabel('Misfit 2-Norm');
ylabel('Model 2-Norm');
for i=1:10:wcount
  H=text(resid(i),mnorm(i)+.2,['  w=',num2str(wval(i))]);
  set(H,'Fontsize',18);
end
%print -deps mwtradeo.eps

figure(11)
bookfonts
for i=1:31
hold on
plot(t(1:N-1),mstore(1:N-1,i)+10*log10(wval(i)),'k')
plot(t(1:N-1),mtrue+10*log10(wval(i)),'-.r')
end
hold off
xlabel('Time (s)')
ylabel('10 log_{10} (w)')
axis tight
%print -deps mwrange.eps

figure(12)
bookfonts
plot(t(1:N-1),mstore(1:N-1,16),'k')
hold on
plot(t(1:N-1),mtrue(1:N-1),'-.r')
xlabel('Time (s)')
ylabel('Acceleration (m/s^2)')
axis tight
%print -deps mwbest.eps