function [output] = hsfilter(theta, Fs, input)
% hsfilter(theta, Fs, input)
% 
% filters the input signal according to head shadowing
% theta is the angle with the frontal plane 
% Fs is the sample rate

theta = theta + 90;
theta0 = 150 ;alfa_min = 0.05 ;		
c = 334;  % speed of sound
a = 0.08; % radius of head												
w0 = c/a;
alfa = 1+ alfa_min/2 + (1- alfa_min/2)* cos(theta/ theta0* pi) ;	

B = [(alfa+w0/Fs)/(1+w0/Fs), (-alfa+w0/Fs)/(1+w0/Fs)] ;	
    % numerator of Transfer Function
A = [1, -(1-w0/Fs)/(1+w0/Fs)] ;	                      
    % denominator of Transfer Function
if (abs(theta) < 90)
 gdelay = - Fs/w0*(cos(theta*pi/180) - 1)  
else
 gdelay = Fs/w0*((abs(theta) - 90)*pi/180 + 1) 
end;
a = (1 - gdelay)/(1 + gdelay);
    % allpass filter coefficient
out_magn = filter(B, A, input);
output = filter([a, 1],[1, a], out_magn);