function [h,g] = HBLPF(No,fp,fs,m);

% [h, g] = HBLPF(No,fp,fs,m);
%
% Purpose:  	Generate Coefficients for a 1/2 band filter (even type I)
%           	of order No (length No+1) using an order No/2 (odd type II)
%           	LPF.
%
% Input Args:
%              No => Filter Order (N must be even & N/2 must be odd)
%              fp => passband width (fp/fs < 1/4)
%              fs => input sample frequency    
%              m  => coeffcient resolution in bits
%
% Output Args:
%              h => actual coefficients
%              g => transform filter coefficients
%
%
%          References:
%          [1]  P. P. Vaidyanathan and T. Q. Nguyen, "A 'Trick' for
%               the design of FIR half-band filters", IEEE Trans. on
%               Circuits & Systems, vol. CAS-34, no. 3, pp. 297-300
%               March, 1987.
%          [2]  R. E. Crochiere and L. R. Rabiner, "Interpolation and
%               decimation of digital signals-A tutorial review", IEEE
%               Proceedings, vol. 69, no. 3, pp. 300-331, March, 1981.
%          [3]  A. Oppenheim, R. Schafer, "Discrete-time signal processing"
%               pp. 469-488, Prentice Hall, 1990. 
%

n=2^14;         	% fft length
xf=(0:n-1)*fs/n;
m=2^m;
h=1;
g=1;

% Check Input Arguments
if fs/fp < 4
        disp('Error #1.1::  fs/fp is not > 4!!!')
        disp('HBLPFprototype failed')
        return;
elseif fs <=0 | fp <= 0 | No <= 0
        disp('Error #2.1:: Negative or zero calling arguements not allowed!!!')
        disp('HBLPFprototype failed')
        return;
elseif rem(No,2) ~= 0  | (rem(No,4) == 0 & rem(No,2) == 0)
        disp('Error #3.1:: No must be EVEN & No/2 must be ODD !!!')
        disp('HBLPFprototype failed')
        return;
end

% Set Up Half-Band Filter Parameters
wp = 2*fp/fs;
Fo = [0 2*wp 1 1];
Mo = [1 1 0 0];
Ng = No/2;

% Compute g(z) Coefficients
g = remez(Ng,Fo,Mo);

% Form h(z) = T[g(z)]
h = zeros(size(1:2*Ng+1));

for k = 1:Ng+1
   h(2*k-1) = g(k)/2;
end

%figure(1)
%subplot(2,1,1)
%zplane(h);
%title('zeros before addition of center coefficient ')
%axis([-2.5 2.5 -2.5 2.5]);
%grid;

% Round coefficients to m bit
h(Ng+1) = .5;
h=round(m*h)/m;
%subplot(2,1,2)
%zplane(h);
%title('zeros after  addition of center coefficient ')
%axis([-2.5 2.5 -2.5 2.5]);
%grid;

% Coefficient conversion to 2's complement
le=ceil(length(h)/2);
hd=h(1:le);
h2=m*hd;
for i=1:le
  if h2(i)<0
  h2(i)=2*m+h2(i);
  end
end
hb=dec2bin(h2);
hh=dec2hex(h2);

% Print Filter Coefficients
format long
Decimal_Filter_Coefficients=hd;
format short
TwosComp_Bin_Coefficients=hb;
TwosComp_Hex_Coefficients=hh;

% Plot Filter Response
%figure(2)
subplot(3,1,1)
m1=min(h);
m2=max(h);
if No < 78
  impz(h);
end
if No >= 78
  x=0:No;
  plot(x,h);
  end
title('Impulse Response of Half-Band Filter')
axis([0 No 1.2*m1 1.05*m2]);
grid;
 
subplot(3,1,2);
f1=20*log10(abs(fft(h,n)));
plot(xf,f1);
axis([0 fs/2 -150 10]);
title('Frequency Response of Half-Band Filter')
xlabel('Frequency (kHz)');
ylabel('Magnitude (dB)');
grid;
 
subplot(3,1,3);
plot(xf,f1);
axis([0 fs/4 -.8 .2]);
title('Passband Response of Half-Band Filter')
xlabel('Frequency (kHz)');
ylabel('Magnitude (dB)');
grid;

