function [out,f] = headmix( in, d)
% Return a dummy-head recording simulation
%
%  Usage:
%    [out,f] = headmix( in, d);
%
%  Inputs:
%    in: NxM matrix, containing one source in every row
%    d: N vector, nth element containing the azimuth of the nth source
%
%  Outputs:
%    out: If there is only one output it returns the mixing
%         filters in an NxN FIR matrix structure,
%         otherwise an NxM matrix containing the mixtures
%    f: Contains the mixing filters in an NxN FIR matrix structure

% Default directions
if nargin == 1
	d = []
else
	if length( d) ~= 2
		error( 'Direction vector size is different than number of sources');
	end
end

% We cannot have more than two mixes
[r,c] = size( in);
if r > 2
	error( 'Cannot have more than two ears now, can we?');
end

% Load head filters
f1 = hrtf( d(1));
f2 = hrtf( d(2));

% If we wan't both mix and mixing filters
if nargout == 2
	% Convolve and mix
	out(1,:) = fftfilt( f1(1,:), in(1,:)) + fftfilt( f2(1,:), in(2,:));
	out(2,:) = fftfilt( f1(2,:), in(1,:)) + fftfilt( f2(2,:), in(2,:));

	% Make mixing set
	f = {f1,f2};
end

% Just return mixing set
if nargout < 2
	f = {f1,f2};
end

%--------------------------------
% Load HRTF for a given azimuth
function f = hrtf( az)

% Bring in 0-360 degrees region
while az > 360
	az = az - 360;
end
while az < 0
	az = az + 360;
end

% 180+x = flipped x
flip = 0;
if az > 180 
	az = az - 180;
	flip = 1;
end

% Load filters
fname = sprintf( '../hrtf/H0e%.3da.dat', az);
fid = fopen( fname, 'r', 'ieee-be');
if fid == -1
	error( sprintf( 'Cannot open file %s', fname));
end
tmp = fread( fid, inf, 'short');
fclose( fid);

% Assign to output and scale
% filters are at 44100, decimate to the default 22050
if ~flip
	f = [decimate( tmp(1:2:end), 2)' ; decimate( tmp(2:2:end), 2)']/32767;
else
	f = [decimate( tmp(2:2:end), 2)' ; decimate( tmp(1:2:end), 2)']/32767;
end