Reading mini-Seed

by

 

m-files for reading mini-Seed files with STEIN-1 or STEIN-2 compression protocol.

minSeedST2()
function signal = minSeedST2()
%-------------------------------------------------------------------------
% Reading mini-Seed data
% signal      :    seismic trace stored in mini-Seed file
%
% Compresion algorithm is STEIM-2 meaning that "Encoding_format=11"
%
% For experience Matlab users, program is self explanatory.
% Good Lack.
%
% Ljupco Jordanovski
% Seismological Observatory at University Scc. Kiril and Methodius
% Skopje,  R. Macedonia.
% ljordanovski@gmail.com
% ------------------------------------------------------------------------
% ------------------SEED data encoding types STEIM-2 ----------------------
%--------------------------------------------------------------------------%--------------------------------------------

END_DATA        = 0;
BgEndian    = false;
LtEndian    = false;

[FileName,PathName] = uigetfile('*.*','Select the Seed ST2-file');
inpFile= fullfile(PathName,FileName);

fid = fopen(inpFile,'rb','ieee-be');
%---------------------   Figure out what is it: Big or Little endian
out                         = fread(fid,20,'*char');
Start_time_Year             = fread(fid,1,'uint16');
if Start_time_Year < 2056
    BgEndian = true;
    fclose(fid);
    fid = fopen(inpFile,'rb','ieee-be');
else
    LtEndian = true;
    fclose(fid);
    fid = fopen(inpFile,'rb','ieee-le');
end
%------------------------------------------------------------
stream_counter  =1;
signal=[];
sequence_counter=-1;
while 1
if END_DATA
        disp('End of reading mSeed file');
        plot(signal);
        fclose(fid);
        return
end

%__________________ exmple:     |000001D | 
    Start_of_block = ftell(fid);
    first_data_line = 0;
       
    Sequence_number     = fread(fid,6,'*char')';
    sequence_counter=sequence_counter+1;
    if isempty(Sequence_number)
        disp('End of reading mSeed file');
        plot(signal);
        return
    end
    %----------------------------------------------------------------------
    Header_type         = fread(fid,1,'*char')';    
    % V - Volume Index Control Header
    % A - Abbreviation Control Header
    % S - Station Control Header
    % T - Time span Control Header
    % D,R or Q - Data Records
    %--------------------------------------------------------------------------
    Continuation_code   = fread(fid,1,'*char')';
 
    % Next, depending of Header_type we read additional information
    switch Header_type
        
        case {'D','R','Q'}
            %disp('Data records')
            stream = Read_Data_record_info();
        case 'V'
            disp('Volume Index Control Header');
        case 'A'
            disp('Abbreviation Control Header');
        case 'S'
            disp('Station Control Header');
        case 'T'
            disp('Time Span Control Headers');
        otherwise
            disp('Unknown Control Header');
            fclose(fid);
            return
    end
end
fclose(fid);
%-----------------------------------------------------------------------
function stream = Read_Data_record_info()

    Station_Code                = fread(fid,5,'*char')';
    Location_id                 = fread(fid,2,'*char')';
    Chann_id                    = fread(fid,3,'*char')';
    Network_code                = fread(fid,2,'*char')';
    Start_time_Year             = fread(fid,1,'uint16');
    Start_time_Day_of_Year      = fread(fid,1,'uint16');
    Start_time_Hours_of_day     = fread(fid,1,'uint8');
    Start_time_Minutes_of_day   = fread(fid,1,'uint8');
    Start_time_Seconds_of_day   = fread(fid,1,'uint8');
    Dummy                       = fread(fid,1,'uint8');
    Start_time_0001_seconds     = fread(fid,1,'uint16');
    Numb_samples                = fread(fid,1,'uint16');
    Sample_rate_factor          = fread(fid,1,'int16');
    Sample_rate_mult            = fread(fid,1,'int16');
    Activity_flag               = fread(fid,1,'uint8');
    I_O_Flag                    = fread(fid,1,'uint8');
    Data_quality_flag           = fread(fid,1,'uint8');
    Num_blockett_follow         = fread(fid,1,'uint8');
    Time_correction             = fread(fid,1,'float32');
    Offset_begin_data           = fread(fid,1,'uint16');
    Offset_first_blockette      = fread(fid,1,'uint16');
    Offset_next_blockette       = Offset_first_blockette;
    seconds = Start_time_Seconds_of_day+0.0001*Start_time_0001_seconds;
    while Offset_next_blockette
        position                = Offset_next_blockette+Start_of_block;
        fseek(fid,position, 'bof');
        Blockette_type          = fread(fid,1,'uint16');
        switch Blockette_type
            case 1000
                Offset_next_blockette   = fread(fid,1,'uint16');
                Encoding_format         = fread(fid,1,'uint8');
                STEIM                   = Encoding_format;
                Word_order              = fread(fid,1,'uint8');
                Data_record_length      = fread(fid,1,'uint8');
                record_size             = 2^Data_record_length;
                Reserved                = fread(fid,1,'uint8');
            case 1001
                Offset_next_blockette   = fread(fid,1,'uint16');
                Timing_quality          = fread(fid,1,'uint8');
                Micro_sec               = fread(fid,1,'uint8');
                Reserved                = fread(fid,1,'uint8');
                Frame_count             = fread(fid,1,'uint8');
            case 100
                Offset_next_blockette   = fread(fid,1,'uint16');
                Actual_Sample_rate      = fread(fid,1,'float32');
                Flags = fread(fid,1,'uint8');
                Reserve_byte            = fread(fid,1,'uint8');
            otherwise
            disp('Unknown Blockette');
            fclose(fid);
            return
        end
    end
        pos = Offset_begin_data+Start_of_block;
        fseek(fid,pos,'bof');
        chk = exist('record_size','var');
        if chk==0
            record_size=4096;
        end
    tot_smpl = 0;
    Number_blocs = (record_size-Offset_begin_data)/64;
    stream = Read_first_data_frame();
    signal = [signal stream];
    tot_smpl=tot_smpl+length(stream);
        if END_DATA 
           return
        end
    for rc = 2:Number_blocs
        stream = Read_rest_data_frame();
        signal = [signal stream];
        tot_smpl=tot_smpl+length(stream);
        if END_DATA && (rc < Number_blocs)
           disp(num2str(tot_smpl))
           return
        end        
    end
%--------------------------------------------------------------------------

    function stream = Read_first_data_frame()
        stream=[];
        stream_counter=0;
        cc          = zeros(1,16,'uint8');
        ww          = zeros(16,4,'uint8');
        pp          = uint32(0);
        for k=1:16
            ww(k,:)   = fread(fid,4,'uint8');
        end
        pp         = Get_cc(ww(1,1:4));
        for k=16:-1:1
                cc(k)           = bitand(pp,3);
                pp              = bitshift(pp,-2);
        end
        first_point                    = Byte_convert_to_32(ww(2,:),4,'L');
        last_point                     = Byte_convert_to_32(ww(3,:),4,'L');
        if ~isempty(signal)
            last                           = signal(end);
        else
            last                           = first_point;
        end
            for k = 4:16
                switch cc(k)
                    case 1
                        data                        = Byte_convert_to_32(ww(k,:),1,'L');
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last                    = stream(stream_counter);
                        end
                    case 2
                        data                        = stein2(ww(k,:),2);
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last = stream(stream_counter);
                        end
                    case 3
                        data                        = stein2(ww(k,:),3);
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last = stream(stream_counter);
                        end
                    otherwise
                        END_DATA = 1;
                end
            end
    end

%--------------------------------------------------------------------------

    function stream = Read_rest_data_frame()
        stream_counter          = 0;
        stream=[];
        if ~isempty(signal)
            last                           = signal(end);
        else
            last                           = first_point;
        end
        cc                  = zeros(1,16,'uint8');
        ww                  = zeros(16,4,'uint8');
        pp                  = uint32(0);
        for k=1:16
            ww(k,:)         = fread(fid,4,'uint8');
        end
        pp                  = Get_cc(ww(1,1:4));
        s=pp;
        index               = 0;
        for k=16:-1:1
            cc(k)           = bitand(pp,3);
            pp              = bitshift(pp,-2);
        end
        if cc(1)>0
            return
        end
            for k = 2:16
                switch cc(k)
                    case 1
                        data                        = Byte_convert_to_32(ww(k,:),1,'L');
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last = stream(stream_counter);
                        end
                    case 2
                        data                        = stein2(ww(k,:),2);
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last = stream(stream_counter);
                        end
                    case 3
                        data                        = stein2(ww(k,:),3);
                        for nn = 1:length(data)
                            stream_counter          = stream_counter+1;
                            stream(stream_counter)  = last+ data(nn);
                            last = stream(stream_counter);
                        end
                    otherwise
                        END_DATA = 1;
                end
            end
    end

end
end
%-------------------------------------------------------------
   function out = Byte_convert_to_32(number8,n,type)
 
   %------------------------------------------  One byte Data    START
   if n == 1
           out= int32(typecast(number8,'int8'));
           return
   end
   %-------------------------------------------  One byte Data     END
   
        switch type 
            case 'B'    % Big Endian
                switch n
                    case 2
                        numberArray = number8;
                        tmp = typecast(numberArray, 'int16');
                    case 4
                        numberArray = number8;
                        tmp = typecast(numberArray, 'int32');
                    otherwise
                        disp(' unknown packing number ')
                end

            case 'L'    % Little Endian
                switch n
                    case 2
                        numberArray = number8(end:-1:1);
                        tmp = typecast(numberArray, 'int16');
                    case 4
                        numberArray = number8(end:-1:1);
                        tmp = typecast(numberArray, 'int32');
                    otherwise
                        disp(' unknown packing number ')
                end
            otherwise
                disp(' Unkown Type Byte ordering ')
        end
   out = int32(tmp);

   end
%-------------------------------------------------------------
   function tmp = Get_cc(Array)
   numberArray = Array(end:-1:1);
   tmp = typecast(numberArray, 'uint32');
   end
   %-----------------------------------------------------------------------
   function out = stein2(nmb8,ck)
   Array = nmb8(end:-1:1);
   n32 = typecast(Array,'uint32');
   c0 = uint8(bitget(n32,32:-1:31));
   kode = c0(1)*2+c0(2);
   switch kode
       case 2
           mask = uint32(15);
           start =7;
           bshift = -4;
           bget = 4;
           if ck == 2
               mask = uint32(32767);
               start = 2;
               bshift = -15;
               bget = 15;
           end
       case 1
           mask = uint32(31);
           start = 6;
           bshift = -5;
           bget = 5;
        
           if ck == 2
               mask = uint32(1073741823);
               start = 1;
               bshift = -30;
               bget = 30;
           end
       case 0
           mask = uint32(63);
           start = 5;
           bshift = -6;
           bget = 6;
       case 3
           mask = uint32(1023);
           start = 3;
           bshift = -10;
           bget = 10;
       otherwise
           disp(' unknown packing number ')
   end
   for k = start:-1:1
       data = bitand(n32,mask);
       sgn  = bitget(data,bget:bget);
       if sgn
           data = bitcmp(data);
           data = bitand(data,mask)+1;
           out(k)=-int32(data);
       else
           out(k) = int32(data);
       end
       n32 = bitshift(n32,bshift);
   end
   end

Contact us