]> Creatis software - CreaPhase.git/blobdiff - octave_packages/audio-1.1.4/auload.m
Add a useful package (from Source forge) for octave
[CreaPhase.git] / octave_packages / audio-1.1.4 / auload.m
diff --git a/octave_packages/audio-1.1.4/auload.m b/octave_packages/audio-1.1.4/auload.m
new file mode 100644 (file)
index 0000000..833e470
--- /dev/null
@@ -0,0 +1,383 @@
+## Copyright (C) 1999 Paul Kienzle
+##
+## This program is free software; you can redistribute it and/or modify
+## it under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2 of the License, or
+## (at your option) any later version.
+##
+## This program is distributed in the hope that it will be useful,
+## but WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+## GNU General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with this program; If not, see <http://www.gnu.org/licenses/>.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{x},@var{fs},@var{sampleformat}] =} auload (@var{filename})
+##
+## Reads an audio waveform from a file given by the string @var{filename}.  
+## Returns the audio samples in data, one column per channel, one row per 
+## time slice.  Also returns the sample rate and stored format (one of ulaw, 
+## alaw, char, int16, int24, int32, float, double). The sample value will be 
+## normalized to the range [-1,1] regardless of the stored format.
+##
+## @example
+##    [x, fs] = auload(file_in_loadpath("sample.wav"));
+##    auplot(x,fs);
+## @end example
+##
+## Note that translating the asymmetric range [-2^n,2^n-1] into the 
+## symmetric range [-1,1] requires a DC offset of 2/2^n. The inverse 
+## process used by ausave requires a DC offset of -2/2^n, so loading and 
+## saving a file will not change the contents.  Other applications may 
+## compensate for the asymmetry in a different way (including previous 
+## versions of auload/ausave) so you may find small differences in 
+## calculated DC offsets for the same file.
+## @end deftypefn
+
+## 2001-09-04 Paul Kienzle <pkienzle@users.sf.net>
+## * skip unknown blocks in WAVE format.
+## 2001-09-05 Paul Kienzle <pkienzle@users.sf.net>
+## * remove debugging stuff from AIFF format.
+## * use data length if it is given rather than reading to the end of file.
+## 2001-12-11 Paul Kienzle <pkienzle@users.sf.net>
+## * use closed interval [-1,1] rather than open interval [-1,1) internally
+
+function [data, rate, sampleformat] = auload(path)
+
+  if (nargin != 1)
+    usage("[x, fs, sampleformat] = auload('filename.ext')");
+  end
+  data = [];       # if error then read nothing
+  rate = 8000;
+  sampleformat = 'ulaw';
+  ext = rindex(path, '.');
+  if (ext == 0)
+    usage('x = auload(filename.ext)');
+  end
+  ext = tolower(substr(path, ext+1, length(path)-ext));
+
+  [file, msg] = fopen(path, 'rb');
+  if (file == -1)
+    error([ msg, ": ", path]);
+  end
+
+  msg = sprintf('Invalid audio header: %s', path);
+  ## Microsoft .wav format
+  if strcmp(ext,'wav') 
+
+    ## Header format obtained from sox/wav.c
+    ## April 15, 1992
+    ## Copyright 1992 Rick Richardson
+    ## Copyright 1991 Lance Norskog And Sundry Contributors
+    ## This source code is freely redistributable and may be used for
+    ## any purpose.  This copyright notice must be maintained. 
+    ## Lance Norskog And Sundry Contributors are not responsible for 
+    ## the consequences of using this software.
+
+    ## check the file magic header bytes
+    arch = 'ieee-le';
+    str = char(fread(file, 4, 'char')');
+    if !strcmp(str, 'RIFF')
+      error(msg);
+    end
+    len = fread(file, 1, 'int32', 0, arch);
+    str = char(fread(file, 4, 'char')');
+    if !strcmp(str, 'WAVE')
+      error(msg);
+    end
+
+    ## skip to the "fmt " section, ignoring everything else
+    while (1)
+      if feof(file)
+       error(msg);
+      end
+      str = char(fread(file, 4, 'char')');
+      len = fread(file, 1, 'int32', 0, arch);
+      if strcmp(str, 'fmt ')
+       break;
+      end
+      fseek(file, len, SEEK_CUR);
+    end
+
+    ## read the "fmt " section
+    formatid = fread(file, 1, 'int16', 0, arch);
+    channels = fread(file, 1, 'int16', 0, arch);
+    rate = fread(file, 1, 'int32', 0, arch);
+    fread(file, 1, 'int32', 0, arch);
+    fread(file, 1, 'int16', 0, arch);
+    bits = fread(file, 1, 'int16', 0, arch);
+    fseek(file, len-16, SEEK_CUR);
+
+    ## skip to the "data" section, ignoring everything else
+    while (1)
+      if feof(file)
+       error(msg);
+      end
+      str = char(fread(file, 4, 'char')');
+      len = fread(file, 1, 'int32', 0, arch);
+      if strcmp(str, 'data')
+       break;
+      end
+      fseek(file, len, SEEK_CUR);
+    end
+
+    if (formatid == 1)
+      if bits == 8
+       sampleformat = 'uchar';
+       precision = 'uchar';
+        samples = len;
+      elseif bits == 16
+       sampleformat = 'int16';
+       precision = 'int16';
+        samples = len/2;
+      elseif bits == 24
+       sampleformat = 'int24';
+       precision = 'int24';
+        samples = len/3;
+      elseif bits == 32
+       sampleformat = 'int32';
+       precision = 'int32';
+        samples = len/4;
+      else
+               error(msg);
+      endif
+    elseif (formatid == 3)
+      if bits == 32
+       sampleformat = 'float';
+       precision = 'float';
+       samples = len/4;
+      elseif bits == 64
+       sampleformat = 'double';
+       precision = 'double';
+       samples = len/8;
+      else
+       error(msg);
+      endif
+    elseif (formatid == 6 && bits == 8)
+      sampleformat = 'alaw';
+      precision = 'uchar';
+      samples = len;
+    elseif (formatid == 7 && bits == 8)
+      sampleformat = 'ulaw';
+      precision = 'uchar';
+      samples = len;
+    else
+      error(msg);
+      return;
+    endif
+
+  ## Sun .au format
+  elseif strcmp(ext, 'au')
+
+    ## Header format obtained from sox/au.c
+    ## September 25, 1991
+    ## Copyright 1991 Guido van Rossum And Sundry Contributors
+    ## This source code is freely redistributable and may be used for
+    ## any purpose.  This copyright notice must be maintained. 
+    ## Guido van Rossum And Sundry Contributors are not responsible for 
+    ## the consequences of using this software.
+
+    str = char(fread(file, 4, 'char')');
+    magic=' ds.';
+    invmagic='ds. ';
+    magic(1) = char(0);
+    invmagic(1) = char(0);
+    if strcmp(str, 'dns.') || strcmp(str, magic)
+      arch = 'ieee-le';
+    elseif strcmp(str, '.snd') || strcmp(str, invmagic)
+      arch = 'ieee-be';
+    else
+      error(msg);
+    end
+    header = fread(file, 1, 'int32', 0, 'ieee-be');
+    len = fread(file, 1, 'int32', 0, 'ieee-be');
+    formatid = fread(file, 1, 'int32', 0, 'ieee-be');
+    rate = fread(file, 1, 'int32', 0, 'ieee-be');
+    channels = fread(file, 1, 'int32', 0, 'ieee-be');
+    fseek(file, header-24, SEEK_CUR); % skip file comment
+
+    ## interpret the sample format
+    if formatid == 1
+      sampleformat = 'ulaw';
+      precision = 'uchar';
+      bits = 12;
+      samples = len;
+    elseif formatid == 2
+      sampleformat = 'uchar';
+      precision = 'uchar';
+      bits = 8;
+      samples = len;
+    elseif formatid == 3
+      sampleformat = 'int16';
+      precision = 'int16';
+      bits = 16;
+      samples = len/2;
+    elseif formatid == 5
+      sampleformat = 'int32';
+      precision = 'int32';
+      bits = 32;
+      samples = len/4;
+    elseif formatid == 6
+      sampleformat = 'float';
+      precision = 'float';
+      bits = 32;
+      samples = len/4;
+    elseif formatid == 7
+      sampleformat = 'double';
+      precision = 'double';
+      bits = 64;
+      samples = len/8;
+    else
+      error(msg);
+    end
+      
+  ## Apple/SGI .aiff format
+  elseif strcmp(ext,'aiff') || strcmp(ext,'aif')
+
+    ## Header format obtained from sox/aiff.c
+    ## September 25, 1991
+    ## Copyright 1991 Guido van Rossum And Sundry Contributors
+    ## This source code is freely redistributable and may be used for
+    ## any purpose.  This copyright notice must be maintained. 
+    ## Guido van Rossum And Sundry Contributors are not responsible for 
+    ## the consequences of using this software.
+    ##
+    ## IEEE 80-bit float I/O taken from
+    ##        ftp://ftp.mathworks.com/pub/contrib/signal/osprey.tar
+    ##        David K. Mellinger
+    ##        dave@mbari.org
+    ##        +1-831-775-1805
+    ##        fax       -1620
+    ##        Monterey Bay Aquarium Research Institute
+    ##        7700 Sandholdt Road
+    ## check the file magic header bytes
+    arch = 'ieee-be';
+    str = char(fread(file, 4, 'char')');
+    if !strcmp(str, 'FORM')
+      error(msg);
+    end
+    len = fread(file, 1, 'int32', 0, arch);
+    str = char(fread(file, 4, 'char')');
+    if !strcmp(str, 'AIFF')
+      error(msg);
+    end
+
+    ## skip to the "COMM" section, ignoring everything else
+    while (1)
+      if feof(file)
+       error(msg);
+      end
+      str = char(fread(file, 4, 'char')');
+      len = fread(file, 1, 'int32', 0, arch);
+      if strcmp(str, 'COMM')
+       break;
+      end
+      fseek(file, len, SEEK_CUR);
+    end
+
+    ## read the "COMM" section
+    channels = fread(file, 1, 'int16', 0, arch);
+    frames = fread(file, 1, 'int32', 0, arch);
+    bits = fread(file, 1, 'int16', 0, arch);
+    exp = fread(file, 1, 'uint16', 0, arch);    % read a 10-byte float
+    mant = fread(file, 2, 'uint32', 0, arch);
+    mant = mant(1) / 2^31 + mant(2) / 2^63;
+    if (exp >= 32768), mant = -mant; exp = exp - 32768; end
+    exp = exp - 16383;
+    rate = mant * 2^exp;
+    fseek(file, len-18, SEEK_CUR);
+
+    ## skip to the "SSND" section, ignoring everything else
+    while (1)
+      if feof(file)
+       error(msg);
+      end
+      str = char(fread(file, 4, 'char')');
+      len = fread(file, 1, 'int32', 0, arch);
+      if strcmp(str, 'SSND')
+       break;
+      end
+      fseek(file, len, SEEK_CUR);
+    end
+    offset = fread(file, 1, 'int32', 0, arch);
+    fread(file, 1, 'int32', 0, arch);
+    fseek(file, offset, SEEK_CUR);
+
+    if bits == 8
+      precision = 'uchar';
+      sampleformat = 'uchar';
+      samples = len - 8;
+    elseif bits == 16
+      precision = 'int16';
+      sampleformat = 'int16';
+      samples = (len - 8)/2;
+    elseif bits == 32
+      precision = 'int32';
+      sampleformat = 'int32';
+      samples = (len - 8)/4;
+    else
+      error(msg);
+    endif
+    
+  ## file extension unknown
+  else
+    error('auload(filename.ext) understands .wav .au and .aiff only');
+  end
+
+  ## suck in all the samples
+  if (samples <= 0) samples = Inf; end
+  if (precision == 'int24')
+    data = fread(file, 3*samples, 'uint8', 0, arch);
+    if (arch == 'ieee-le')
+      data = data(1:3:end) + data(2:3:end) * 2^8 + cast(typecast(cast(data(3:3:end), 'uint8'), 'int8'), 'double') * 2^16;
+    else
+      data = data(3:3:end) + data(2:3:end) * 2^8 + cast(typecast(cast(data(1:3:end), 'uint8'), 'int8'), 'double') * 2^16;
+    endif
+  else
+    data = fread(file, samples, precision, 0, arch);
+  endif
+  fclose(file);
+
+  ## convert samples into range [-1, 1)
+  if strcmp(sampleformat, 'alaw')
+    alaw = [ \
+            -5504,  -5248,  -6016,  -5760,  -4480,  -4224,  -4992,  -4736, \
+            -7552,  -7296,  -8064,  -7808,  -6528,  -6272,  -7040,  -6784, \
+            -2752,  -2624,  -3008,  -2880,  -2240,  -2112,  -2496,  -2368, \
+            -3776,  -3648,  -4032,  -3904,  -3264,  -3136,  -3520,  -3392, \
+           -22016, -20992, -24064, -23040, -17920, -16896, -19968, -18944, \
+           -30208, -29184, -32256, -31232, -26112, -25088, -28160, -27136, \
+           -11008, -10496, -12032, -11520,  -8960,  -8448,  -9984,  -9472, \
+           -15104, -14592, -16128, -15616, -13056, -12544, -14080, -13568, \
+             -344,   -328,   -376,   -360,   -280,   -264,   -312,   -296, \
+             -472,   -456,   -504,   -488,   -408,   -392,   -440,   -424, \
+              -88,    -72,   -120,   -104,    -24,     -8,    -56,    -40, \
+             -216,   -200,   -248,   -232,   -152,   -136,   -184,   -168, \
+            -1376,  -1312,  -1504,  -1440,  -1120,  -1056,  -1248,  -1184, \
+            -1888,  -1824,  -2016,  -1952,  -1632,  -1568,  -1760,  -1696, \
+             -688,   -656,   -752,   -720,   -560,   -528,   -624,   -592, \
+             -944,   -912,  -1008,   -976,   -816,   -784,   -880,   -848 ];
+    alaw = ([ alaw,-alaw]+0.5)/32767.5;
+    data = alaw(data+1);
+  elseif strcmp(sampleformat, 'ulaw')
+    data = mu2lin(data, 0);
+  elseif strcmp(sampleformat, 'uchar')
+    ## [ 0, 255 ] -> [ -1, 1 ]
+    data = data/127.5 - 1;
+  elseif strcmp(sampleformat, 'int16')
+    ## [ -32768, 32767 ] -> [ -1, 1 ]
+    data = (data+0.5)/32767.5;
+  elseif strcmp(sampleformat, 'int32')
+    ## [ -2^31, 2^31-1 ] -> [ -1, 1 ]
+    data = (data+0.5)/(2^31-0.5);
+  end
+  data = reshape(data, channels, length(data)/channels)';
+
+endfunction
+
+%!demo
+%! [x, fs] = auload(file_in_loadpath("sample.wav"));
+%! auplot(x,fs);