% MTX = steer2HarmMtx(HARMONICS, ANGLES, REL_PHASES)
%
% Compute a steering matrix (maps a directional basis set onto the
% angular Fourier harmonics). HARMONICS is a vector specifying the
% angular harmonics contained in the steerable basis/filters. ANGLES
% (optional) is a vector specifying the angular position of each filter.
% REL_PHASES (optional, default = 'even') specifies whether the harmonics
% are cosine or sine phase aligned about those positions.
% The result matrix is suitable for passing to the function STEER.

% Eero Simoncelli, 7/96.

function mtx = steer2HarmMtx(harmonics, angles, evenorodd)

%%=================================================================
%%% Optional Parameters:

if (exist('evenorodd') ~= 1)
evenorodd = 'even';
end

% Make HARMONICS a row vector
harmonics = harmonics(:)';

numh = 2*size(harmonics,2) - any(harmonics == 0);

if (exist('angles') ~= 1)
angles = pi * [0:numh-1]'/numh;
else
angles = angles(:);
end

%%=================================================================

if isstr(evenorodd)
if strcmp(evenorodd,'even')
evenorodd = 0;
elseif strcmp(evenorodd,'odd')
evenorodd = 1;
else
error('EVEN_OR_ODD should be the string EVEN or ODD');
end
end

%% Compute inverse matrix, which maps Fourier components onto
%% steerable basis.
imtx = zeros(size(angles,1),numh);
col = 1;
for h=harmonics
args = h*angles;
if (h == 0)
imtx(:,col) = ones(size(angles));
col = col+1;
elseif evenorodd
imtx(:,col) = sin(args);
imtx(:,col+1) = -cos(args);
col = col+2;
else
imtx(:,col) = cos(args);
imtx(:,col+1) = sin(args);
col = col+2;
end
end

r = rank(imtx);
if (( r ~= numh ) & ( r ~= size(angles,1) ))
fprintf(2,'WARNING: matrix is not full rank');
end

mtx = pinv(imtx);