function [y,Fs,bits,opts]=auread(file,ext)
%AUREAD Read NeXT/SUN (".au") sound file.
%   Y=AUREAD(AUFILE) loads a sound file specified by the string AUFILE,
%   returning the sampled data in y. The ".au" extension is appended
%   if no extension is given.  Amplitude values are in the range [-1,+1].
%   Supports multi-channel data in the following formats:
%   8-bit mu-law, 8-, 16-, and 32-bit linear, and floating point.
%
%   [Y,Fs,BITS]=AUREAD(AUFILE) returns the sample rate (Fs) in Hertz
%   and the number of bits per sample (BITS) used to encode the
%   data in the file.
%
%   [...]=AUREAD(AUFILE,N) returns only the first N samples from each
%       channel in the file.
%   [...]=AUREAD(AUFILE,[N1 N2]) returns only samples N1 through N2 from
%       each channel in the file.
%   SIZ=AUREAD(AUFILE,'size') returns the size of the audio data contained
%       in the file in place of the actual audio data, returning the
%       vector SIZ=[samples channels].
%
%   See also AUWRITE, WAVREAD.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.7 $  $Date: 1998/02/04 20:19:21 $

%   D. Orofino, 10/95

if nargin>2,
  error('Too many input arguments.');
end
% Append extension if it's missing:
if isempty(findstr(file,'.')),
  file=[file '.au'];  % or .snd
end
fid=fopen(file,'rb','b');   % Big-endian
if fid == -1,
  error('Can''t open NeXT/SUN sound file for input.');
end

% No optional information is read from this file format.
opts = [];

% Read header:
snd = read_sndhdr(fid);
if ~isstruct(snd), error(snd); end
% Parse structure info for return to user:
Fs = snd.rate; bits = snd.bits;

% Determine if caller wants data:
if nargin<2, ext=[]; end    % Default - read all samples
exts=prod(size(ext));
if strncmp(lower(ext),'size',exts),
  % Caller doesn't want data - just data size:
  fclose(fid);
  y = [snd.samples snd.chans];
  return;
elseif exts>2,
  error('Index range must be specified as a scalar or 2-element vector.');
elseif (exts==1),
  ext=[1 ext];  % Prepend start sample index
end

% Read data:
snd = read_sndata(fid,snd,ext);
if ~isstruct(snd), error(snd); end
y = snd.data; fclose(fid);

% end of auread()


% ------------------------------------------------------------------------
% Private functions:
% ------------------------------------------------------------------------

% READ_SNDHDR: Read sound file header structure
%   Assumes fid points to the start of an open file.
%   Returns a structure if successful.
function snd = read_sndhdr(fid)
% Read file header:
snd.magic = setstr(fread(fid,4,'char')');
if ~strcmp(snd.magic,'.snd'),
  snd='Not a NeXT/Sun sound file.'; return;
end
snd.offset = fread(fid,1,'uint32');
snd.databytes = fread(fid,1,'uint32');
snd.format = fread(fid,1,'uint32');
snd.rate = fread(fid,1,'uint32');
snd.chans = fread(fid,1,'uint32');

% Directly determine how long info string is:
info_len = snd.offset-24;
[info,cnt] = fread(fid,info_len,'char');
snd.info = deblank(setstr(info'));
if cnt~=info_len,
  snd='Error while reading sound file.'; return;
end

% For some file types, the .databytes field seems to have
% an invalid file length.  Directly determine file length
% for all cases:
fseek(fid,0,1);             % Go to end of file
file_len = ftell(fid);      % Get position in bytes
fseek(fid,snd.offset,-1);   % Reposition file pointer
snd.databytes = file_len-snd.offset;

% Interpret format:
if snd.format==1,
  snd.bits = 8;  % 8-bit mu-law
elseif snd.format==2,
  snd.bits = 8;  % 8-bit linear
elseif snd.format==3,
  snd.bits = 16; % 16-bit linear
elseif snd.format==5,
  snd.bits = 32; % 32-bit linear
elseif snd.format==6,
  snd.bits = 32; % Single precision
elseif snd.format==7,
  snd.bits = 64; % Double-precision
else
  snd='Unrecognized data format.'; return;
end
% Determine # of samples per channel:
snd.samples = snd.databytes*8/snd.bits/snd.chans;
if snd.samples~=fix(snd.samples),
  snd='Truncated data file.'; return;
end

return;

% READ_SNDATA: Read sound file header structure
%   Assumes fid points to the start of an open file.
%   Returns a structure if successful.
%
function new_snd = read_sndata(fid,snd,ext)
  SamplesPerChannel = snd.samples;
  BytesPerSample = snd.bits/8;

  % Interpret format:
  if snd.format==1,
    dtype='uchar';  % 8-bit mu-law
  elseif snd.format==2,
    dtype='int8';  % 8-bit linear
  elseif snd.format==3,
    dtype='int16'; % 16-bit linear
  elseif snd.format==5,
    dtype='int32'; % 32-bit linear
  elseif snd.format==6,
    dtype='float'; % Single precision
  elseif snd.format==7,
    dtype='double'; % Double-precision
  else
    new_snd='Unrecognized data format.'; return;
  end

  % Determine sample range to read:
  if isempty(ext),
    ext = [1 SamplesPerChannel];    % Return all samples
  else
    if prod(size(ext))~=2,
      new_snd='Sample limit vector must have 2 elements.'; return;
    end
    if ext(1)<1 | ext(2)>SamplesPerChannel,
      new_snd='Sample limits out of range.'; return;
    end
    if ext(1)>ext(2),
      new_snd='Sample limits must be given in ascending order.'; return;
    end
  end
  % Skip over leading samples:
  if ext(1)>1,
    % Skip over leading samples, if specified:
    status = fseek(fid,BytesPerSample*(ext(1)-1)*snd.chans,0);
    if status==-1,
      new_snd='Error in file format.'; return;
    end
  end
  % Read desired data:
  nSPCext = ext(2)-ext(1)+1; % # samples per channel in extraction range
  extSamples = snd.chans*nSPCext;
  [data,cnt] = fread(fid, [snd.chans nSPCext], dtype);
  if cnt~=extSamples,
    snd = 'Data file is truncated.'; return;
  end
  % Rearrange data into a matrix with one channel per column:
  data = data';

  % Convert and normalize data range:
  if snd.format==1, % 8-bit mu-law
    data = mu2lin(data);
  elseif snd.format==2, % 8-bit linear
    data = data*2^(-7);
  elseif snd.format==3, % 16-bit linear
    data = data*2^(-15);
  elseif snd.format==5, % 32-bit linear
    data = data*2^(-31);
  elseif (snd.format==6) | (snd.format==7), % Float/Double
    a=min(min(data)); b=max(max(data));
    data=(data-a)/(b-a)*2-1;
  end

  new_snd = snd;
  new_snd.data = data;
return;

% end of auread.m
