Thread: response spectra

Started: 2006-09-20 15:49:39
Last activity: 2006-09-23 02:29:03
Topics: SAC Help
feigelso@rohan.sdsu.edu
2006-09-20 15:49:39
Sac-help,
Can SAC make response spectra? If not, how can sac files be converted for
MATLAB to read?
Thanks for the help,
Leah



  • Frederik Tilmann
    2006-09-23 02:29:03
    Alternative matlab routines for reading sac files for Sun binary order.



    http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=5546&objectType=FILE

    http://www.mathworks.com/matlabcentral/fileexchange/loadFile.do?objectId=3859&objectType=FILE


    Attached file is the same as first link but for PC byte order.

    Frederik


    feigelso<at>rohan.sdsu.edu wrote:
    Sac-help,
    Can SAC make response spectra? If not, how can sac files be converted for
    MATLAB to read?
    Thanks for the help,
    Leah


    _______________________________________________
    sac-help mailing list
    sac-help<at>iris.washington.edu
    http://www.iris.washington.edu/mailman/listinfo/sac-help


    --
    Frederik Tilmann
    Bullard Laboratories Tel. +44 1223 765545
    Department of Earth Sciences Fax. +44 1223 360779
    University of Cambridge email: tilmann<at>esc.cam.ac.uk
    Madingley Road
    Cambridge CB3 0EZ
    UK



    function [SACdata,SeisData,sacfiles]=sacpc2mat(varargin)
    % [SACdata,SeisData,filenames] = SACPCMAT('file1','file2',..., 'filen' )
    %
    % reads n SAC files file1, file2, filen (SAC files are assumed to have
    % PC byte order) and converts them to matlab
    % format. The filenames can contain globbing characters (e.g. * and ?).
    % These are expanded and all matching files loaded.
    %
    % SACPCMAT( cellarray ) where cellarray={'file1','file2',...,'filen'}
    % is equivalent to the standard form.
    %
    % SACdata is an n x 1 struct array containing the header variables
    % in the same format as is obtained by using MAT function
    % of SAC2000.
    % SACdata(i).trcLen contains the number of samples.
    %
    % SeisData is an m x n array (where m=max(npts1, npts2, ...) )
    % containing the actual data.
    %
    % filenames is a n x 1 string cell array with the filenames actually read.
    %
    % Note that writing
    %
    % [SACdata,SeisData] = sacsun2mat('file1','file2',..., 'filen' )
    %
    % is equivalent to the following sequence
    %
    % sac2000
    % READ file1 file2 .. filen
    % MAT
    %
    % (in fact the failure of above sequence to work properly on my
    % system motivated this script).
    %
    %
    % SACPC2MAT was written by F Tilmann (tilmann<at>esc.cam.ac.uk)
    % based on sac_sun2pc_mat by C. D. Saragiotis (I copied the
    % routines doing the actual work from this code but
    % used a different header structure and made the routine
    % flexible).
    % It was tested on MATLAB5 on a PC but
    % should work on newer versions, too.
    %
    % (C) 2004
    %

    F = 4-1; % float byte-size minus 1;
    K = 8-1; % alphanumeric byte-size minus 1
    L = 4-1; % long integer byte-size minus 1;

    fnames={};
    for i=1:nargin
    if ischar(varargin{i})
    fnames=cell([fnames; cellstr(varargin{i})]);
    elseif iscellstr(varargin{i}) & size(varargin{i},1)==1
    fnames=cell([fnames; varargin{i}']);
    elseif iscellstr(varargin{i}) & size(varargin{i},2)==1
    fnames=cell([fnames; varargin{i}]);
    end
    end
    % expand globs
    sacfiles={};k=1;
    for i=1:length(fnames)
    dirlist=dir(fnames{i});
    for j=1:length(dirlist)
    if ~dirlist(j).isdir
    sacfiles{k,1}=dirlist(j).name;
    k=k+1;
    end
    end
    end

    maxnpts=0;
    for i=1:length(sacfiles)
    fid=fopen(sacfiles{i},'rb');
    if fid==-1
    error(sprintf('Could not open SAC file %s',fnames{i}))
    end
    SACdata(i,1)=readSacHeader(fid,F,K,L);
    npts=SACdata(i).trcLen;
    if npts>maxnpts
    maxnpts=npts;
    end
    fprintf('Processing file %d: %s\n',i,sacfiles{i});
    SeisData(npts,i)=0; % Magnify seis matrix if necessary
    SeisData(:,i)=[ readSacData(fid,npts,F+1); zeros(maxnpts-npts,1)]; % Have to pad with zeros if new data have less data points than some previously read file
    end

    function hdr = readSacHeader(fileId,F,K,L)
    % hdr = readSacHeader(FID)
    % sacReadAlphaNum reads the SAC-format header fields and returns most of them.
    %
    % The output variable, 'hdr' is a structure, whose fields are
    % the fields as in the SACdata structure generated by SAC's
    % matlab command MAT
    %
    % Created by C. D. Saragiotis, August 5th, 2003, modified by F Tilmann
    headerBytes = 632;
    chrctr = fread(fileId,headerBytes,'uchar');
    chrctr = chrctr(:)';
    % Read floats
    hdr.times.delta = sacReadFloat(chrctr(1:1+F)); % increment between evenly spaced samples
    %hdr.DEPMIN = sacReadFloat(chrctr(5:5+F)); % MINimum value of DEPendent variable
    %hdr.DEPMAX = sacReadFloat(chrctr(9:9+F)); % MAXimum value of DEPendent variable
    % (not currently used) SCALE = sacReadFloat(chrctr(13:13+F)); % Mult SCALE factor for dependent variable
    %hdr.ODELTA = sacReadFloat(chrctr(17:17+F)); % Observd increment if different than DELTA
    hdr.times.b = sacReadFloat(chrctr(21:21+F)); % Begining value of the independent variable
    hdr.times.e = sacReadFloat(chrctr(25:25+F)); % Ending value of the independent variable
    hdr.times.o = sacReadFloat(chrctr(29:29+F)); % event Origin time
    hdr.times.a = sacReadFloat(chrctr(33:33+F)); % first Arrival time
    hdr.times.t0 = sacReadFloat(chrctr(41:41+F));
    hdr.times.t1 = sacReadFloat(chrctr(45:45+F));
    hdr.times.t2 = sacReadFloat(chrctr(49:49+F));
    hdr.times.t3 = sacReadFloat(chrctr(53:53+F));
    hdr.times.t4 = sacReadFloat(chrctr(57:57+F));
    hdr.times.t5 = sacReadFloat(chrctr(61:61+F));
    hdr.times.t6 = sacReadFloat(chrctr(65:65+F));
    hdr.times.t7 = sacReadFloat(chrctr(69:69+F));
    hdr.times.t8 = sacReadFloat(chrctr(73:73+F));
    hdr.times.t9 = sacReadFloat(chrctr(77:77+F));
    hdr.times.f = sacReadFloat(chrctr(81:81+F)); % Fini of event time


    hdr.response = sacReadFloat(reshape(chrctr(85:85+10*(F+1)-1),F+1,10)');

    hdr.station.stla = sacReadFloat(chrctr(125:125+F)); % STation LAttitude
    hdr.station.stlo = sacReadFloat(chrctr(129:129+F)); % STation LOngtitude
    hdr.station.stel = sacReadFloat(chrctr(133:133+F)); % STation ELevation
    hdr.station.stdp = sacReadFloat(chrctr(137:137+F)); % STation DePth below surface

    hdr.event.ecla = sacReadFloat(chrctr(141:141+F)); % EVent LAttitude
    hdr.event.evlo = sacReadFloat(chrctr(145:145+F)); % EVent LOngtitude
    hdr.event.evel = sacReadFloat(chrctr(149:149+F)); % EVent ELevation
    hdr.event.evdp = sacReadFloat(chrctr(153:153+F)); % EVent DePth below surface
    hdr.event.mag = sacReadFloat(chrctr(157:157+F)); % EVent DePth below surface

    userdata=sacReadFloat(reshape(chrctr(161:161+10*(F+1)-1),F+1,10)');

    hdr.evsta.dist = sacReadFloat(chrctr(201:201+F)); % station to event DISTance (km)
    hdr.evsta.az = sacReadFloat(chrctr(205:205+F)); % event to station AZimuth (degrees)
    hdr.evsta.baz = sacReadFloat(chrctr(209:209+F)); % station to event AZimuth (degrees)
    hdr.evsta.gcarc = sacReadFloat(chrctr(213:213+F)); % station to event Great Circle ARC length (degrees)

    %hdr.DEPMEN = sacReadFloat(chrctr(225:225+F)); % MEaN value of DEPendent variable

    hdr.station.cmpaz = sacReadFloat(chrctr(229:229+F)); % CoMPonent AZimuth (degrees clockwise from north)
    hdr.station.cmpinc = sacReadFloat(chrctr(233:233+F)); % CoMPonent INCident angle (degrees from vertical)

    hdr.llnl.xminimum = sacReadFloat(chrctr(237:237+L));
    hdr.llnl.xmaximum = sacReadFloat(chrctr(241:241+L));
    hdr.llnl.yminimum = sacReadFloat(chrctr(245:245+L));
    hdr.llnl.ymaximum = sacReadFloat(chrctr(249:249+L));

    % Read long integers
    hdr.event.nzyear = sacReadLong(chrctr(281:281+L)); % GMT YEAR
    hdr.event.nzjday = sacReadLong(chrctr(285:285+L)); % GMT julian DAY
    hdr.event.nzhour = sacReadLong(chrctr(289:289+L)); % GMT HOUR
    hdr.event.nzmin = sacReadLong(chrctr(293:293+L)); % GMT MINute
    hdr.event.nzsec = sacReadLong(chrctr(297:297+L)); % GMT SECond
    hdr.event.nzmsec = sacReadLong(chrctr(301:301+L)); % GMT MilliSECond

    hdr.llnl.norid = sacReadLong(chrctr(309:309+L));
    hdr.llnl.nevid = sacReadLong(chrctr(313:313+L));

    hdr.trcLen = sacReadLong(chrctr(317:317+L)); % Number of PoinTS per data component

    hdr.llnl.nwfid = sacReadLong(chrctr(325:325+L));
    hdr.llnl.nxsize = sacReadLong(chrctr(329:329+L));
    hdr.llnl.nysize = sacReadLong(chrctr(333:333+L));

    hdr.descrip.iftype = sacReadLong(chrctr(341:341+L)); % File TYPE
    hdr.descrip.idep = sacReadLong(chrctr(345:345+L)); % type of DEPendent variable
    hdr.descrip.iztype = sacReadLong(chrctr(349:349+L)); % reference time equivalence

    hdr.descrip.iinst = sacReadLong(chrctr(357:357+L)); % type of recording INSTrument

    % Before there were floats read here?
    hdr.descrip.istreg = sacReadLong(chrctr(361:361+F)); % STation geographic REGion
    hdr.descrip.ievreg = sacReadLong(chrctr(365:365+F)); % EVent geographic REGion
    hdr.descrip.ievtyp = sacReadLong(chrctr(369:369+F)); % EVent geographic REGion
    hdr.descrip.iqual = sacReadLong(chrctr(373:373+F)); % QUALity of data
    hdr.descrip.isynth = sacReadLong(chrctr(377:377+F)); % SYNTHetic data flag

    hdr.event.imagtyp = sacReadLong(chrctr(381:381+F));
    hdr.event.imagsrc = sacReadLong(chrctr(385:385+F));

    % no logical SAC header variables in matlab SACdata structure!
    % previous version set these defaults
    % $$$ % SAC defaults
    % $$$ hdr.IEVTYPE = 'IUNKN';
    % $$$ %IEVTYPE= sacReadFloat(chrctr(369:369+F)); % EVent TYPE
    % $$$ hdr.LEVEN = 'TRUE'; % true, if data are EVENly spaced (required)
    % $$$ %LEVEN= sacReadFloat(chrctr(421:521+F)); % true, if data are EVENly spaced (required)
    % $$$ hdr.LPSPOL = 'FALSE'; % true, if station components have a PoSitive POLarity
    % $$$ %LPSPOL= sacReadFloat(chrctr(425:425+F)); % true, if station components have a PoSitive POLarity
    % $$$ hdr.LOVROK = 'FALSE'; % true, if it is OK to OVeRwrite this file in disk
    % $$$ %LOVROK= sacReadFloat(chrctr(429:429+F)); % true, if it is OK to OVeRwrite this file in disk
    % $$$ hdr.LCALDA = 'TRUE'; % true, if DIST, AZ, BAZ and GCARC are to be calculated from station and event coordinates
    % $$$ %LCALDA= sacReadFloat(chrctr(433:433+F)); % true, if DIST, AZ, BAZ and GCARC are to be calculated from station and event coordinates

    % Read alphanumeric data
    hdr.station.kstnm = sacReadAlphaNum(chrctr(441:441+K)); % STation NaMe
    hdr.event.kevnm = sacReadAlphaNum(chrctr(449:449+2*(K+1)-1)); % EVent NaMe
    %hdr.KHOLE = sacReadAlphaNum(chrctr(465:465+K)); % HOLE identification, if nuclear event
    hdr.times.ko = sacReadAlphaNum(chrctr(473:473+K)); % event Origin time identification
    hdr.times.ka = sacReadAlphaNum(chrctr(481:481+K)); % first Arrival time identification

    hdr.times.kt0 = sacReadAlphaNum(chrctr(489:489+K));
    hdr.times.kt1 = sacReadAlphaNum(chrctr(497:497+K));
    hdr.times.kt2 = sacReadAlphaNum(chrctr(505:505+K));
    hdr.times.kt3 = sacReadAlphaNum(chrctr(513:513+K));
    hdr.times.kt4 = sacReadAlphaNum(chrctr(521:521+K));
    hdr.times.kt5 = sacReadAlphaNum(chrctr(529:529+K));
    hdr.times.kt6 = sacReadAlphaNum(chrctr(537:537+K));
    hdr.times.kt7 = sacReadAlphaNum(chrctr(545:545+K));
    hdr.times.kt8 = sacReadAlphaNum(chrctr(553:553+K));
    hdr.times.kt9 = sacReadAlphaNum(chrctr(561:561+K));


    hdr.times.kf = sacReadAlphaNum(chrctr(569:569+K)); % Fini identification
    kuser0 = sacReadAlphaNum(chrctr(577:577+K)); % USER-defined variable storage area
    kuser1 = sacReadAlphaNum(chrctr(585:585+K)); % USER-defined variable storage area
    kuser2 = sacReadAlphaNum(chrctr(593:593+K)); % USER-defined variable storage area
    hdr.station.kcmpnm = sacReadAlphaNum(chrctr(601:601+K)); % CoMPonent NaMe
    hdr.station.knetwk = sacReadAlphaNum(chrctr(609:609+K)); % name of seismic NETWorK
    %hdr.KDATRD = sacReadAlphaNum(chrctr(617:617+K)); % DATa Recording Date onto the computer
    %hdr.KINST = sacReadAlphaNum(chrctr(625:625+K)); % generic name of recording INSTrument

    usercell=num2cell(userdata);
    [usercell{find(userdata==-12345)}]=deal([]);
    [hdr.user(1:10).data]=deal(usercell{:});
    [hdr.user(1:10).label]=deal(kuser0, kuser1,kuser2,[], [], [], [], [], [], []);


    function X = readSacData(fid,N,F)
    % function data = readSacData('filename',NPTS,floatByteSize)
    chrctr = fread(fid,N*F,'uchar');
    x=reshape(chrctr,F,N)';
    %x
    X = sacReadFloat(x);


    function lNumber = sacReadLong(cb)
    % reads long integers (4 bytes long)
    % cb is the character buffer
    cb = cb(:);
    % SUN sac format use: lNumber = (256.^(3:-1:0))*cb;
    lNumber = (256.^(0:3))*cb; % changed line
    if lNumber == -12345, lNumber = []; end

    function fNumber = sacReadFloat(cb)
    % reads floats (4 bytes long)
    % cb is the character buffer
    C = size(cb,1);
    % SUN sac format use: stringOfBitSequence = [dec2bin(cb(:,1),8) dec2bin(cb(:,2),8) dec2bin(cb(:,3),8) dec2bin(cb(:,4),8)];
    stringOfBitSequence = [dec2bin(cb(:,4),8) dec2bin(cb(:,3),8) ...
    dec2bin(cb(:,2),8) dec2bin(cb(:,1),8)]; %changed line
    bitSequence = stringOfBitSequence=='1';
    fSign = -2*bitSequence(:,1)+1;
    fExp = bitSequence(:,2:9)*(2.^(7:-1:0)') - 127;
    fMantissa = [ones(C,1) bitSequence(:,10:32)]*(2.^(0:-1:-23)');
    fNumber = fSign.*fMantissa.*(2.^fExp);
    isZeroCheck = sum(bitSequence')'==0;
    fNumber = fNumber.*(~isZeroCheck);
    if C==1 & fNumber == -12345, fNumber = []; end


    function alphaNum = sacReadAlphaNum(cb)
    % reads alphanumeric data (8 or 16 bytes long). If it cb is empty, it returns a ' '
    % cb is the character buffer
    K = max(size(cb));
    alphaNumTemp = char(cb);
    if K == 8
    if alphaNumTemp == '-12345 '
    alphaNum = [];
    else
    alphaNum = alphaNumTemp;
    end
    else
    if K == 16
    if alphaNumTemp == '-12345 -12345 ' | alphaNumTemp == '-12345 '
    alphaNum = [];
    else
    alphaNum = alphaNumTemp;
    end
    end
    end

05:28:12 v.01697673