% MATLAB function computes the filter coefficients of 
% a type 1 FIR filter using least-squares inversion
% AUTHOR: Balu Santhanam
% DATE:   02/09/01
% USAGE: 
% [hopt,edB] = lsfilt(freq,magn,ord);
% hopt : optimal filter coefficents
% edB  : Least squares design error measuure
% freq : vector containing frequencies in [0,pi];
% magn : vector containing normalized magnitudes 
% ord  : order of filter desired
%*******************************************************
function [hopt,edB] = lsfilt(freq,magn,ord);
% Error checking
if nargin < 3
   error('Insufficient input arguments')
elseif length(find(freq < 0 | freq > pi)) ~= 0
   error('Frequencies have to lie in [0,pi]')
elseif length(find(magn > 1 | magn < 0)) ~= 0
   error('Magnitudes have to lie in +ve integer')
elseif length(freq) == 0 | length(magn) == 0
   error('Null Input')
elseif isnumeric(freq) == 0 | isnumeric(magn) == 0 | isnumeric(ord) == 0
   error('Non numeric input')
elseif ((mod(ord,1) ~= 0) | (mod(ord,2) == 0))
   error('Filter order has to be a odd +ve integer')
elseif any(freq) == 0 | any(magn) == 0 
   error('Zero Input')
end
% Create A matrix
A = [];
for i = 1:1:length(freq)
    vec = cos(freq(i)*[0:1:round((ord-1)/2)]); 
    A = [A ; vec];
end 
% Least squares inversion 
b = magn(:); a = A\b;
%  Generate optimal filter coefficients
hopt(round((ord+1)/2)) = a(1);
hopt([1:1:round((ord-1)/2)]) = 0.5*a([round((ord+1)/2):-1:2]);
hopt([round((ord+3)/2):1:ord]) = hopt([round((ord-1)/2):-1:1]);
edB = 20*log10(norm(A*a - b));
return
