diff options
Diffstat (limited to 'matfile.scm')
| -rw-r--r-- | matfile.scm | 187 |
1 files changed, 187 insertions, 0 deletions
diff --git a/matfile.scm b/matfile.scm new file mode 100644 index 0000000..2e3ff15 --- /dev/null +++ b/matfile.scm @@ -0,0 +1,187 @@ +; "matfile.scm", Read MAT-File Format version 4 (MATLAB) +; Copyright (c) 2001, 2002, 2003 Aubrey Jaffer +; +;Permission to copy this software, to modify it, to redistribute it, +;to distribute modified versions, and to use it for any purpose is +;granted, subject to the following restrictions and understandings. +; +;1. Any copy made of this software must include this copyright notice +;in full. +; +;2. I have made no warranty or representation that the operation of +;this software will be error-free, and I am under no obligation to +;provide any services, by way of maintenance, update, or otherwise. +; +;3. In conjunction with products arising from the use of this +;material, there shall be no use of my name in any advertising, +;promotional, or sales literature without prior written consent in +;each case. + +(require 'array) +(require 'byte) +(require 'byte-number) +(require-if 'compiling 'string-case) ; string-ci->symbol used by matfile:load + +;;@code{(require 'matfile)} +;;@ftindex matfile +;;@ftindex matlab +;; +;;@uref{http://www.mathworks.com/access/helpdesk/help/pdf_doc/matlab/matfile_format.pdf} +;; +;;@noindent +;;This package reads MAT-File Format version 4 (MATLAB) binary data +;;files. MAT-files written from big-endian or little-endian computers +;;having IEEE format numbers are currently supported. Support for files +;;written from VAX or Cray machines could also be added. +;; +;;@noindent +;;The numeric and text matrix types handled; support for @dfn{sparse} +;;matrices awaits a sample file. + +(define (bytes->long lst) + (bytes->integer lst -4)) +(define (bytes->short lst) + (bytes->integer lst -2)) +(define (bytes->ushort lst) + (bytes->integer lst 2)) + +;;Version 4 MAT-file endianness cannot be detected solely from the +;;first word; it is ambiguous when 0. +(define (matfile:read-matrix port) + (define null (integer->char 0)) + (define (read1 endian type mrows ncols imagf namlen) + (set! type (bytes->long type)) + (set! mrows (bytes->long mrows)) + (set! ncols (bytes->long ncols)) + (set! imagf (bytes->long imagf)) + (set! namlen (+ -1 (bytes->long namlen))) + (let ((d-prot (modulo (quotient type 10) 10)) + (d-endn (case (quotient type 1000) + ((0 ;ieee-little-endian + 2 ;vax-d-float + 3) endian) ;vag-g-float + ((1 ;ieee-big-endian + 4) (- endian)) ;cray + (else #f))) + (m-type (case (modulo type 10) + ((0) 'numeric) + ((1) 'text) + ((2) 'sparse) + (else #f)))) + (define d-leng (case d-prot + ((0) 8) + ((1 2) 4) + ((3 4) 2) + ((5) 1) + (else #f))) + (define d-conv (case d-prot + ((0) (case (quotient type 1000) + ((0 1) bytes->ieee-double) + ((2) bytes->vax-d-double) + ((3) bytes->vax-g-double) + ((4) bytes->cray-double))) + ((1) (case (quotient type 1000) + ((0 1) bytes->ieee-float) + ((2) bytes->vax-d-float) + ((3) bytes->vax-g-float) + ((4) bytes->cray-float))) + ((2) bytes->long) + ((3) bytes->short) + ((4) bytes->ushort) + ((5) (if (eqv? 'text m-type) + (lambda (lst) (integer->char (byte-ref lst 0))) + (lambda (lst) (byte-ref lst 0)))) + (else #f))) + ;;(@print d-leng d-endn m-type type mrows ncols imagf namlen d-conv) + (cond ((and (= 0 (modulo (quotient type 100) 10) (quotient type 65536)) + d-leng d-endn m-type + (<= 0 imagf 1) + (< 0 mrows #xFFFFFF) + (< 0 ncols #xFFFFFF) + (< 0 namlen #xFFFF)) + (set! imagf (case imagf ((0) #f) ((1) #t))) + (let ((namstr (make-string namlen)) + (mat (case m-type + ((numeric) (create-array + (case d-prot + ((0) ((if imagf Ac64 Ar64))) + ((1) ((if imagf Ac32 Ar32))) + ((2) (As32)) + ((3) (As16)) + ((4) (Au16)) + ((5) (Au8)) + (else (slib:error 'p 'type d-prot))) + mrows ncols)) + ((text) (create-array "." mrows ncols)) + ((sparse) (slib:error 'sparse '?))))) + (do ((idx 0 (+ 1 idx))) + ((>= idx namlen)) + (string-set! namstr idx (read-char port))) + ;;(@print namstr) + (if (not (eqv? null (read-char port))) + (slib:error 'matfile 'string 'missing null)) + (do ((jdx 0 (+ 1 jdx))) + ((>= jdx ncols)) + (do ((idx 0 (+ 1 idx))) + ((>= idx mrows)) + (array-set! mat (d-conv (read-bytes (* d-endn d-leng) port)) + idx jdx))) + (if imagf + (do ((jdx 0 (+ 1 jdx))) + ((>= jdx ncols)) + (do ((idx 0 (+ 1 idx))) + ((>= idx mrows)) + (array-set! mat + (+ (* (d-conv (read-bytes (* d-endn d-leng) port)) + +i) + (array-ref mat idx jdx)) + idx jdx)))) + (list namstr mat))) + (else #f)))) + ;;(trace read1) + (let* ((type (read-bytes 4 port)) + (mrows (read-bytes 4 port)) + (ncols (read-bytes 4 port)) + (imagf (read-bytes 4 port)) + (namlen (read-bytes 4 port))) + ;;Try it with either endianness: + (or (read1 1 type mrows ncols imagf namlen) + (read1 -1 + (bytes-reverse type) + (bytes-reverse mrows) + (bytes-reverse ncols) + (bytes-reverse imagf) + (bytes-reverse namlen))))) + +;;@body @1 should be a string naming an existing file containing a +;;MATLAB Version 4 MAT-File. The @0 procedure reads matrices from the +;;file and returns a list of the results; a list of the name string and +;;array for each matrix. +(define (matfile:read filename) + (call-with-open-ports + (open-file filename 'rb) + (lambda (port) + (do ((mat (matfile:read-matrix port) (matfile:read-matrix port)) + (mats '() (cons mat mats))) + ((or (not mat) (eof-object? (peek-char port))) + (if (and (null? mats) (not mat)) + '() + (reverse (cons mat mats)))))))) + +;;@body @1 should be a string naming an existing file containing a +;;MATLAB Version 4 MAT-File. The @0 procedure reads matrices from the +;;file and defines the @code{string-ci->symbol} for each matrix to its +;;corresponding array. @0 returns a list of the symbols defined. +(define (matfile:load filename) + (require 'string-case) + (let ((mats (matfile:read filename))) + (for-each (lambda (nam-mat) + (and nam-mat + (slib:eval + (list 'define + (string-ci->symbol (car nam-mat)) + (list 'quote (cadr nam-mat)))))) + mats) + (map string-ci->symbol (map car mats)))) + +;;(trace-all "/home/jaffer/slib/matfile.scm") |