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]

《新程序員》：云原生和全面數字化實踐50位技術專家共同創作，文字、視頻、音頻交互閱讀

總結

以上是生活随笔為你收集整理的Matlab画地球剖面图,分享用matlab显示地震记录的波形变面积图的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。