function out = vsp(in, hrirl,hrirr, pos)

% vsp           - virtual source positioning, using HRIRs
%
% out = vsp(in, hrirl, hrirr, pos)
%
% input parameters:
% in:   input signal (mono)
% hrirl: matrix with hrir's in azimuth plane (left ear)
% hrirr: matrix with hrir's in azimuth plane (right ear)
% pos:  vector with positions of the virtual source [deg]
%
% Try:  
%  load kemar_l; load kemar_r; in=wavread('???);
%  pos=interp1(0:18:72, [0 90 0 -90 0], 0:72);
%  out=vsp(in,ir_l,ir_r,pos);

% (c) 2002 Piotr Majdak

N=length(in);
win=size(pos,2);    % get the number of windows
winS=ceil(N/win);   % get the window size
M=size(hrirl,1)-1;  % get the length of impulse response
in=[ zeros(M,1); in; zeros(winS*win-N,1)]; % zero padding of input
pos=round(mod(pos,360)*(size(hrirl,2)-1)/360+1); % deg to index
yl=zeros(win*winS+M,1); 
yr=zeros(win*winS+M,1);

h=waitbar(0, 'virtual sound positioning in progress...');
for ii=1:win
    y=conv(hrirl(:,pos(ii)),in((ii-1)*winS+1:ii*winS)); % filter
    yl((ii-1)*winS+1:ii*winS+M)=yl((ii-1)*winS+1:ii*winS+M)+y; % overlap and add
    y=conv(hrirr(:,pos(ii)),in((ii-1)*winS+1:ii*winS));
    yr((ii-1)*winS+1:ii*winS+M)=yr((ii-1)*winS+1:ii*winS+M)+y;
    waitbar(ii/win,h);
end
close(h);

out(:,1)=yl;
out(:,2)=yr;
out=out/max(max(abs(out))); % normalize