octave_packages/audio-1.1.4/auload.m

   1 ## Copyright (C) 1999 Paul Kienzle
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify
   4 ## it under the terms of the GNU General Public License as published by
   5 ## the Free Software Foundation; either version 2 of the License, or
   6 ## (at your option) any later version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful,
   9 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  11 ## GNU General Public License for more details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License
  14 ## along with this program; If not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {[@var{x},@var{fs},@var{sampleformat}] =} auload (@var{filename})
  18 ##
  19 ## Reads an audio waveform from a file given by the string @var{filename}.
  20 ## Returns the audio samples in data, one column per channel, one row per
  21 ## time slice.  Also returns the sample rate and stored format (one of ulaw,
  22 ## alaw, char, int16, int24, int32, float, double). The sample value will be
  23 ## normalized to the range [-1,1] regardless of the stored format.
  24 ##
  25 ## @example
  26 ##    [x, fs] = auload(file_in_loadpath("sample.wav"));
  27 ##    auplot(x,fs);
  28 ## @end example
  29 ##
  30 ## Note that translating the asymmetric range [-2^n,2^n-1] into the
  31 ## symmetric range [-1,1] requires a DC offset of 2/2^n. The inverse
  32 ## process used by ausave requires a DC offset of -2/2^n, so loading and
  33 ## saving a file will not change the contents.  Other applications may
  34 ## compensate for the asymmetry in a different way (including previous
  35 ## versions of auload/ausave) so you may find small differences in
  36 ## calculated DC offsets for the same file.
  37 ## @end deftypefn
  38
  39 ## 2001-09-04 Paul Kienzle <pkienzle@users.sf.net>
  40 ## * skip unknown blocks in WAVE format.
  41 ## 2001-09-05 Paul Kienzle <pkienzle@users.sf.net>
  42 ## * remove debugging stuff from AIFF format.
  43 ## * use data length if it is given rather than reading to the end of file.
  44 ## 2001-12-11 Paul Kienzle <pkienzle@users.sf.net>
  45 ## * use closed interval [-1,1] rather than open interval [-1,1) internally
  46
  47 function [data, rate, sampleformat] = auload(path)
  48
  49   if (nargin != 1)
  50     usage("[x, fs, sampleformat] = auload('filename.ext')");
  51   end
  52   data = [];       # if error then read nothing
  53   rate = 8000;
  54   sampleformat = 'ulaw';
  55   ext = rindex(path, '.');
  56   if (ext == 0)
  57     usage('x = auload(filename.ext)');
  58   end
  59   ext = tolower(substr(path, ext+1, length(path)-ext));
  60
  61   [file, msg] = fopen(path, 'rb');
  62   if (file == -1)
  63     error([ msg, ": ", path]);
  64   end
  65
  66   msg = sprintf('Invalid audio header: %s', path);
  67   ## Microsoft .wav format
  68   if strcmp(ext,'wav')
  69
  70     ## Header format obtained from sox/wav.c
  71     ## April 15, 1992
  72     ## Copyright 1992 Rick Richardson
  73     ## Copyright 1991 Lance Norskog And Sundry Contributors
  74     ## This source code is freely redistributable and may be used for
  75     ## any purpose.  This copyright notice must be maintained.
  76     ## Lance Norskog And Sundry Contributors are not responsible for
  77     ## the consequences of using this software.
  78
  79     ## check the file magic header bytes
  80     arch = 'ieee-le';
  81     str = char(fread(file, 4, 'char')');
  82     if !strcmp(str, 'RIFF')
  83       error(msg);
  84     end
  85     len = fread(file, 1, 'int32', 0, arch);
  86     str = char(fread(file, 4, 'char')');
  87     if !strcmp(str, 'WAVE')
  88       error(msg);
  89     end
  90
  91     ## skip to the "fmt " section, ignoring everything else
  92     while (1)
  93       if feof(file)
  94         error(msg);
  95       end
  96       str = char(fread(file, 4, 'char')');
  97       len = fread(file, 1, 'int32', 0, arch);
  98       if strcmp(str, 'fmt ')
  99         break;
 100       end
 101       fseek(file, len, SEEK_CUR);
 102     end
 103
 104     ## read the "fmt " section
 105     formatid = fread(file, 1, 'int16', 0, arch);
 106     channels = fread(file, 1, 'int16', 0, arch);
 107     rate = fread(file, 1, 'int32', 0, arch);
 108     fread(file, 1, 'int32', 0, arch);
 109     fread(file, 1, 'int16', 0, arch);
 110     bits = fread(file, 1, 'int16', 0, arch);
 111     fseek(file, len-16, SEEK_CUR);
 112
 113     ## skip to the "data" section, ignoring everything else
 114     while (1)
 115       if feof(file)
 116         error(msg);
 117       end
 118       str = char(fread(file, 4, 'char')');
 119       len = fread(file, 1, 'int32', 0, arch);
 120       if strcmp(str, 'data')
 121         break;
 122       end
 123       fseek(file, len, SEEK_CUR);
 124     end
 125
 126     if (formatid == 1)
 127       if bits == 8
 128         sampleformat = 'uchar';
 129         precision = 'uchar';
 130         samples = len;
 131       elseif bits == 16
 132         sampleformat = 'int16';
 133         precision = 'int16';
 134         samples = len/2;
 135       elseif bits == 24
 136         sampleformat = 'int24';
 137         precision = 'int24';
 138         samples = len/3;
 139       elseif bits == 32
 140         sampleformat = 'int32';
 141         precision = 'int32';
 142         samples = len/4;
 143       else
 144         error(msg);
 145       endif
 146     elseif (formatid == 3)
 147       if bits == 32
 148         sampleformat = 'float';
 149         precision = 'float';
 150         samples = len/4;
 151       elseif bits == 64
 152         sampleformat = 'double';
 153         precision = 'double';
 154         samples = len/8;
 155       else
 156         error(msg);
 157       endif
 158     elseif (formatid == 6 && bits == 8)
 159       sampleformat = 'alaw';
 160       precision = 'uchar';
 161       samples = len;
 162     elseif (formatid == 7 && bits == 8)
 163       sampleformat = 'ulaw';
 164       precision = 'uchar';
 165       samples = len;
 166     else
 167       error(msg);
 168       return;
 169     endif
 170
 171   ## Sun .au format
 172   elseif strcmp(ext, 'au')
 173
 174     ## Header format obtained from sox/au.c
 175     ## September 25, 1991
 176     ## Copyright 1991 Guido van Rossum And Sundry Contributors
 177     ## This source code is freely redistributable and may be used for
 178     ## any purpose.  This copyright notice must be maintained.
 179     ## Guido van Rossum And Sundry Contributors are not responsible for
 180     ## the consequences of using this software.
 181
 182     str = char(fread(file, 4, 'char')');
 183     magic=' ds.';
 184     invmagic='ds. ';
 185     magic(1) = char(0);
 186     invmagic(1) = char(0);
 187     if strcmp(str, 'dns.') || strcmp(str, magic)
 188       arch = 'ieee-le';
 189     elseif strcmp(str, '.snd') || strcmp(str, invmagic)
 190       arch = 'ieee-be';
 191     else
 192       error(msg);
 193     end
 194     header = fread(file, 1, 'int32', 0, 'ieee-be');
 195     len = fread(file, 1, 'int32', 0, 'ieee-be');
 196     formatid = fread(file, 1, 'int32', 0, 'ieee-be');
 197     rate = fread(file, 1, 'int32', 0, 'ieee-be');
 198     channels = fread(file, 1, 'int32', 0, 'ieee-be');
 199     fseek(file, header-24, SEEK_CUR); % skip file comment
 200
 201     ## interpret the sample format
 202     if formatid == 1
 203       sampleformat = 'ulaw';
 204       precision = 'uchar';
 205       bits = 12;
 206       samples = len;
 207     elseif formatid == 2
 208       sampleformat = 'uchar';
 209       precision = 'uchar';
 210       bits = 8;
 211       samples = len;
 212     elseif formatid == 3
 213       sampleformat = 'int16';
 214       precision = 'int16';
 215       bits = 16;
 216       samples = len/2;
 217     elseif formatid == 5
 218       sampleformat = 'int32';
 219       precision = 'int32';
 220       bits = 32;
 221       samples = len/4;
 222     elseif formatid == 6
 223       sampleformat = 'float';
 224       precision = 'float';
 225       bits = 32;
 226       samples = len/4;
 227     elseif formatid == 7
 228       sampleformat = 'double';
 229       precision = 'double';
 230       bits = 64;
 231       samples = len/8;
 232     else
 233       error(msg);
 234     end
 235
 236   ## Apple/SGI .aiff format
 237   elseif strcmp(ext,'aiff') || strcmp(ext,'aif')
 238
 239     ## Header format obtained from sox/aiff.c
 240     ## September 25, 1991
 241     ## Copyright 1991 Guido van Rossum And Sundry Contributors
 242     ## This source code is freely redistributable and may be used for
 243     ## any purpose.  This copyright notice must be maintained.
 244     ## Guido van Rossum And Sundry Contributors are not responsible for
 245     ## the consequences of using this software.
 246     ##
 247     ## IEEE 80-bit float I/O taken from
 248     ##        ftp://ftp.mathworks.com/pub/contrib/signal/osprey.tar
 249     ##        David K. Mellinger
 250     ##        dave@mbari.org
 251     ##        +1-831-775-1805
 252     ##        fax       -1620
 253     ##        Monterey Bay Aquarium Research Institute
 254     ##        7700 Sandholdt Road
 255
 256     ## check the file magic header bytes
 257     arch = 'ieee-be';
 258     str = char(fread(file, 4, 'char')');
 259     if !strcmp(str, 'FORM')
 260       error(msg);
 261     end
 262     len = fread(file, 1, 'int32', 0, arch);
 263     str = char(fread(file, 4, 'char')');
 264     if !strcmp(str, 'AIFF')
 265       error(msg);
 266     end
 267
 268     ## skip to the "COMM" section, ignoring everything else
 269     while (1)
 270       if feof(file)
 271         error(msg);
 272       end
 273       str = char(fread(file, 4, 'char')');
 274       len = fread(file, 1, 'int32', 0, arch);
 275       if strcmp(str, 'COMM')
 276         break;
 277       end
 278       fseek(file, len, SEEK_CUR);
 279     end
 280
 281     ## read the "COMM" section
 282     channels = fread(file, 1, 'int16', 0, arch);
 283     frames = fread(file, 1, 'int32', 0, arch);
 284     bits = fread(file, 1, 'int16', 0, arch);
 285     exp = fread(file, 1, 'uint16', 0, arch);    % read a 10-byte float
 286     mant = fread(file, 2, 'uint32', 0, arch);
 287     mant = mant(1) / 2^31 + mant(2) / 2^63;
 288     if (exp >= 32768), mant = -mant; exp = exp - 32768; end
 289     exp = exp - 16383;
 290     rate = mant * 2^exp;
 291     fseek(file, len-18, SEEK_CUR);
 292
 293     ## skip to the "SSND" section, ignoring everything else
 294     while (1)
 295       if feof(file)
 296         error(msg);
 297       end
 298       str = char(fread(file, 4, 'char')');
 299       len = fread(file, 1, 'int32', 0, arch);
 300       if strcmp(str, 'SSND')
 301         break;
 302       end
 303       fseek(file, len, SEEK_CUR);
 304     end
 305     offset = fread(file, 1, 'int32', 0, arch);
 306     fread(file, 1, 'int32', 0, arch);
 307     fseek(file, offset, SEEK_CUR);
 308
 309     if bits == 8
 310       precision = 'uchar';
 311       sampleformat = 'uchar';
 312       samples = len - 8;
 313     elseif bits == 16
 314       precision = 'int16';
 315       sampleformat = 'int16';
 316       samples = (len - 8)/2;
 317     elseif bits == 32
 318       precision = 'int32';
 319       sampleformat = 'int32';
 320       samples = (len - 8)/4;
 321     else
 322       error(msg);
 323     endif
 324
 325   ## file extension unknown
 326   else
 327     error('auload(filename.ext) understands .wav .au and .aiff only');
 328   end
 329
 330   ## suck in all the samples
 331   if (samples <= 0) samples = Inf; end
 332   if (precision == 'int24')
 333     data = fread(file, 3*samples, 'uint8', 0, arch);
 334     if (arch == 'ieee-le')
 335       data = data(1:3:end) + data(2:3:end) * 2^8 + cast(typecast(cast(data(3:3:end), 'uint8'), 'int8'), 'double') * 2^16;
 336     else
 337       data = data(3:3:end) + data(2:3:end) * 2^8 + cast(typecast(cast(data(1:3:end), 'uint8'), 'int8'), 'double') * 2^16;
 338     endif
 339   else
 340     data = fread(file, samples, precision, 0, arch);
 341   endif
 342   fclose(file);
 343
 344   ## convert samples into range [-1, 1)
 345   if strcmp(sampleformat, 'alaw')
 346     alaw = [ \
 347              -5504,  -5248,  -6016,  -5760,  -4480,  -4224,  -4992,  -4736, \
 348              -7552,  -7296,  -8064,  -7808,  -6528,  -6272,  -7040,  -6784, \
 349              -2752,  -2624,  -3008,  -2880,  -2240,  -2112,  -2496,  -2368, \
 350              -3776,  -3648,  -4032,  -3904,  -3264,  -3136,  -3520,  -3392, \
 351             -22016, -20992, -24064, -23040, -17920, -16896, -19968, -18944, \
 352             -30208, -29184, -32256, -31232, -26112, -25088, -28160, -27136, \
 353             -11008, -10496, -12032, -11520,  -8960,  -8448,  -9984,  -9472, \
 354             -15104, -14592, -16128, -15616, -13056, -12544, -14080, -13568, \
 355               -344,   -328,   -376,   -360,   -280,   -264,   -312,   -296, \
 356               -472,   -456,   -504,   -488,   -408,   -392,   -440,   -424, \
 357                -88,    -72,   -120,   -104,    -24,     -8,    -56,    -40, \
 358               -216,   -200,   -248,   -232,   -152,   -136,   -184,   -168, \
 359              -1376,  -1312,  -1504,  -1440,  -1120,  -1056,  -1248,  -1184, \
 360              -1888,  -1824,  -2016,  -1952,  -1632,  -1568,  -1760,  -1696, \
 361               -688,   -656,   -752,   -720,   -560,   -528,   -624,   -592, \
 362               -944,   -912,  -1008,   -976,   -816,   -784,   -880,   -848 ];
 363     alaw = ([ alaw,-alaw]+0.5)/32767.5;
 364     data = alaw(data+1);
 365   elseif strcmp(sampleformat, 'ulaw')
 366     data = mu2lin(data, 0);
 367   elseif strcmp(sampleformat, 'uchar')
 368     ## [ 0, 255 ] -> [ -1, 1 ]
 369     data = data/127.5 - 1;
 370   elseif strcmp(sampleformat, 'int16')
 371     ## [ -32768, 32767 ] -> [ -1, 1 ]
 372     data = (data+0.5)/32767.5;
 373   elseif strcmp(sampleformat, 'int32')
 374     ## [ -2^31, 2^31-1 ] -> [ -1, 1 ]
 375     data = (data+0.5)/(2^31-0.5);
 376   end
 377   data = reshape(data, channels, length(data)/channels)';
 378
 379 endfunction
 380
 381 %!demo
 382 %! [x, fs] = auload(file_in_loadpath("sample.wav"));
 383 %! auplot(x,fs);