extract.m

来自「对sgy格式的地震数据」· M 代码 · 共 84 行

function [value]=extract(filename,hw)
%EXTRACT: Extract a header word from a segy file 
%
%   This function return the numerical values
%   of a SEGY header word. hw is the desired header 
%
%   example:    cdp = extract('data','cdp') will provide
%               the cdp numbers of each trace in 'data' 
%               See documentation to see the rest of the header
%               words.
%
%   example:    of = extract('data','offset') will extract offset
%               of each trace.
%
%   example:    D = extract('data','trace') will extract the traces
%               into a matrix called D.
%
%
%  M.D.Sacchi, July 1997, Dept. of Physics, UofA.
%        
%  [e-mail]
%
%  Modified by geoer. successful 2007-7-11-20-12.
 

   FID = fopen(filename,'r','n');      % n if for native,n的含义:Numeric format of the machine on which MATLAB is running (the default)

   segy=segy_struct;                   % load the definitions of
                                       % the header words

   count=count_struct;                 % load the position of   
                                       % each word in the header (in bytes)

   status = fseek(FID,count.ns,'bof'); % go to begging of file,to read the ns

   ns = fread(FID,1,segy.ns);          % read ns from first trace

   if strcmp(hw,'ns') == 0;            % get only ns and return,ns较特殊，直接给出并返回
    else
     value = ns;
      return;
       end

   total = 60+ns;                      % total num of 4bytes words 
                                       % in the header

   if strcmp(hw,'trace') == 0;                    
    elements = 1;                     %如果道头字不是trace,那么到时候读取的就是一个元素
     word=zeros(1,1);
      else
      elements = ns;                  %如果道头字是trace，则读取整道，
       word=zeros(ns,1);
        end

   hc=eval(strcat('count.',hw));      %得到要取的道头字的位置和数据类型
   hp=eval(strcat('segy.',hw));


   max_traces = 999999;               %最大读取道数

   for j=1:max_traces;

     position = total*(j-1)*4+hc;     %将要读取的道头字的位置
      status = fseek(FID,position,'bof');  %这里有一个严重的bug,就是tracl读完文件后，再读下一道的时候，fseek跳到eof,依然是合法的，
                                           %其他道头字在读取的时候，则直接失败，跳出.所以,必须进行文件末
                                           %尾的检测.
       if status == 0                      % check status and read,0成功，-1失败     
        word = fread(FID,elements,hp);
       if isempty(word)                    %modified by geoer
           return;
       end
        value(:,j)=word(:,1);
        j = j + 1; 
       else
       return;                       % 读取失败的话，就返回
      end 
   end 

fclose(FID);