24小时热门版块排行榜

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

baobiao007

木虫 (职业作家)

中国特色

应助: 201 (大学生)
金币: 6482.7
散金: 557
红花: 40
帖子: 3050
在线: 1009.9小时
虫号: 505962
注册: 2008-02-18
专业: 应用地球物理学

[交流] 分享用matlab显示地震记录的波形变面积图已有9人参与

wigb函数是网上找的，自己添了个主程序方便使用。

主程序
function PlotRecord_Wigb(filename,nx,nt,dx,dt,zy)
%波形变面积显示地震记录
% clear; clf;clc;
% nx=200; %道数
% nt=2001; %采样点数
% dx=10.;    %道间距
% dt=0.001; %时间采样间隔
% zy=1.0;    %增益
% filename='r200-2001-1ms.dat';
clf;
fid=fopen(filename,'r');
[A,COUNT]=fread(fid,[nt,nx],'float32');
fclose(fid);

x=0:dx: (nx-1)*dx;
z=0:dt: (nt-1)*dt;
wigb(A,zy,x,z);

end

子函数
function wigb (a,scal,x,z,amx)
%WIGB: Plot seismic data using wiggles.
%
%  WIGB(a,scal,x,z,amx)
%
%  IN a:    seismic data
%       scale: multiple data by scale
%       x:    x-axis;
%       z:    vertical axis (time or depth)
%
%    If only 'a' is enter, 'scal,x,z,amn,amx' are decided automatically;
%    otherwise, 'scal' is a scalar; 'x, z' are vectors for annotation in
%    offset and time, amx are the amplitude range.
%
%
%  Author(s): Xingong Li (Exxon-Mobil)
%  Copyright 1998-2003 Xingong
%  Revision: 1.2  Date: Dec/2002
%
%  Signal Analysis and Imaging Group (SAIG)
%  Department of Physics, UofA
%

if nargin == 0, nx=10;nz=10; a = rand(nz,nx)-0.5; end;

[nz,nx]=size(a);

trmx= max(abs(a));
if (nargin <= 4); amx=mean(trmx);  end;
if (nargin <= 2); x=[1:nx]; z=[1:nz]; end;
if (nargin <= 1); scal =1; end;

if nx <= 1 ; disp(' ERR: PlotWig: nx has to be more than 1');return;end;

% take the average as dx
dx1 = abs(x(2:nx)-x(1:nx-1));
dx = median(dx1);

dz=z(2)-z(1);
xmx=max(max(a)); xmn=min(min(a));

if scal == 0; scal=1; end;
a = a * dx /amx;
a = a * scal;

fprintf(' PlotWig: data range [%f, %f], plotted max %f \n',xmn,xmx,amx);

% set display range
x1=min(x)-2.0*dx; x2=max(x)+2.0*dx;
z1=min(z)-dz; z2=max(z)+dz;

set(gca,'NextPlot','add','Box','on', ...
  'XLim', [x1 x2], ...
  'YDir','reverse', ...
  'YLim',[z1 z2]);

fillcolor = [0 0 0];
linecolor = [0 0 0];
linewidth = 0.1;

z=z'; % input as row vector
zstart=z(1);
zend  =z(nz);

for i=1:nx,

  if trmx(i) ~= 0; % skip the zero traces
tr=a(:,i); % --- one scale for all section
    s = sign(tr) ;
    i1= find( s(1:nz-1) ~= s(2:nz) ); % zero crossing points
npos = length(i1);

%12/7/97
zadd = i1 + tr(i1) ./ (tr(i1) - tr(i1+1)); %locations with 0 amplitudes
aadd = zeros(size(zadd));

[zpos,vpos] = find(tr >0);
[zz,iz] = sort([zpos; zadd]); % indices of zero point plus positives
aa = [tr(zpos); aadd];
aa = aa(iz);

% be careful at the ends
if tr(1)>0, a0=0; z0=1.00;
else, a0=0; z0=zadd(1);
end;
if tr(nz)>0, a1=0; z1=nz;
else, a1=0; z1=max(zadd);
end;

zz = [z0; zz; z1; z0];
aa = [a0; aa; a1; a0];

zzz = zstart + zz*dz -dz;

patch( aa+x(i) , zzz,  fillcolor);

line( 'Color',[1 1 1],'EraseMode','background',  ...
      'Xdata', x(i)+[0 0], 'Ydata',[zstart zend]); % remove zero line

%'LineWidth',linewidth, ...
%12/7/97 'Xdata', x(i)+[0 0], 'Ydata',[z0 z1]*dz); % remove zero line

line( 'Color',linecolor,'EraseMode','background',  ...
   'LineWidth',linewidth, ...
   'Xdata', tr+x(i), 'Ydata',z); % negatives line

else % zeros trace
line( 'Color',linecolor,'EraseMode','background',  ...
   'LineWidth',linewidth, ...
      'Xdata', [x(i) x(i)], 'Ydata',[zstart zend]);
end;
end;

[ Last edited by baobiao007 on 2012-4-9 at 10:41 ]

回复此楼