function tonyc2(n,wo)
% tonyc(n,wo) forms coefficients of two path polyphase filter
% n must be odd (3,5,7.....) 
% wo is normalized band edge frequency .5>wo>.25 
% coefficient list is a0 a1 a2 a3... 
% even indices belong to path 0, odd indices belong to path 1  
% based on paper "Digital Signal Processing Schemes for Efficient Interpolation and Decimation"
% by Valenzuela and Constantinides, IEE Proceedings, Dec 1983

if (n-2*floor(n/2)) ==0
   disp(' filter order must be odd, please try again')
  return;
end;
%step 1
wt=4*pi*(wo-0.25);
k=tan(0.25*(pi-wt));
k=k*k;

kk=sqrt(1-k*k);

e=0.5*(1-sqrt(kk))/(1+sqrt(kk));

q=e+2*(e^5)+15*(e^9)+150*(e^13);
ww=zeros(1,(n-1)/2);
aa=ww;
cc=ww;
%step 2
for ii=1:(n-1)/2

ww(ii)=2*(q^0.25)*(sin(pi*ii/n)-(q^2)*sin(3*pi*ii/n));
ww(ii)=ww(ii)/(1-2*(q*cos(2*pi*ii/n)-(q^4)*cos(4*pi*ii/n)));

wwsq=ww(ii)*ww(ii);
aa(ii)=sqrt(((1-wwsq*k)*(1-wwsq/k)))/(1+wwsq);
cc(ii)=(1-aa(ii))/(1+aa(ii));
end
coef=cc
ordr1=floor((n-1)/4);
ordr2=ordr1;
if n-1-4*ordr1 ~= 0 
	ordr1=ordr1+1;
end
den1=zeros(ordr1,3);
den2=zeros(ordr2,3);

for ii=1:ordr1
den1(ii,:)=[1 0 cc(2*ii-1)];
end

for ii=1:ordr2
den2(ii,:)=[1 0 cc(2*ii)];
end

h0=[1];
for ii=1:ordr1
h0=conv(h0,den1(ii,:));
end

h1=[1];
for ii=1:ordr2
h1=conv(h1,den2(ii,:));
end
h1=[h1 0];
g0=fliplr(h0);
g1=fliplr(h1);

[tp,ww]=freqz(g0,h0,512);
[bt,ww]=freqz(g1,h1,512);

hh=0.5*(tp+bt);
clg;
plot(0:1/1024:.5-1/1024,unwrap(angle(tp))/(pi));axis([0 .5 -(2*ordr1+1) 0]);hold;
plot(0:1/1024:.5-1/1024,unwrap(angle(bt))/(pi),'c');axis([0 .5 -(2*ordr1+1) 0]);grid;
title('phase of each path in two path filter');
xlabel('normalized frequency (f/fs)');
ylabel('normalized phase (theta/pi)');
hold;pause

plot(0:1/1024:.5-1/1024,20*log10(abs(hh)));axis([0 .5 -100 10]);grid;
title('Magnitude response of two path filter');
xlabel('normalized frequency (f/fs)');
ylabel('log magnitude (dB)');
pause

plot(0:1/1024:.5-1/1024,unwrap(angle(hh))/(2*pi));grid;%axis([0 .5 -3 0]);
title('Phase response of two path filter');
xlabel('normalized frequency (f/fs)');
ylabel('normalized phase (THETA/2*pi)');
pause

gd=conv([-1 1],unwrap(angle(hh)));
id=find(gd<0);
gd(id)=gd(id-1);
plot(0:1/1024:.5,gd);grid;

title('group delay of two path filter');
xlabel('normalized frequency (f/fs)');
ylabel('normalized delay (TAU/Ts)');
pause

[hh0,t]=impz(g0,h0,10*n);
[hh1,t]=impz(g1,h1,10*n);
resp=0.5*(hh0+hh1);
%resp=impz([1 2.5 2.5 1],[1 0 .4 0],20);
stem(resp);

title('impulse response of filter');
xlabel('normalized time')
ylabel('amplitude')
pause

rr=xcorr(resp,resp);
stem(rr)
title('impulse response of composite filter');
xlabel('normalized time')
ylabel('amplitude')
pause
fyy=20*log10(abs(fft(rr,512)));
plot(0:1/512:.5,fyy(1:257)-max(fyy));
axis([0 .5 -80 10]);grid;