## Fred Harris polyphase recursive all-pass filter function tony_des2

Started by 3 weeks ago4 replieslatest reply 3 weeks ago126 views

Hi All,

Thanks,

Chris

[ - ]

Here it is

function tony_des_2(action)
% tony_des_2 forms coefficients of two path recursive polyphase filter
% filter can be halfband quadrature mirror, polynomials in Z^2
% filter can be bandwidth tuned,        Z^-1 => (1-b*z)/(z-b),       b=(1-tan(bw))/(1+tan(bw))
% filter can be center frequency tuned, z^-1 => (-1/z)*(1-cz)/(z-c), c=cos(theta_c)
% filter can be zero-packed             Z^-1 ->  z^-k,               k even
% order of filter is odd (3,5,7.....)
% wo is normalized band edge frequency .5>wo>.25
%
% %%%%%%%%% this program calls tonycxx.m %%%%%%%%%%%%%
%
% based on paper "Digital Signal Processing Schemes for Efficient Interpolation and Decimation"
% by Valenzuela and Constantinides, IEE Proceedings, Dec 1983
% Script file written by fred harris of SDSU. Copyright 2002

global N_
global ws_
global wb_
global wc_
global kk_
if nargin<1

figure(1)
clf reset
set(1,'name','                Two-Path Allpass Filter Design by fred harris, SDSU','numbertitle','off');

responseplot=axes('units','normalized','position',[.05 .4 .9 .55]);

N_=uicontrol('style','edit','units','normalized','position',[0.18 0.1 0.07 0.07],'string','5');
Ntext=uicontrol('style','text','units','normalized','position',[0.04 0.1 0.13 0.065],'string','Prototype Order (Odd)');

ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string','0.300');
wstext=uicontrol('style','text','units','normalized','position',[0.04 0.17 0.13 0.065],'string','Stopband Edge (0.25-0.5)');

wb_=uicontrol('style','edit','units','normalized','position',[0.42 0.17 0.07 0.07],'string','0.250');
wbtext=uicontrol('style','text','units','normalized','position',[0.28 0.17 0.13 0.065],'string','Passband Edge');

wc_=uicontrol('style','edit','units','normalized','position',[0.66 0.17 0.07 0.07],'string','0.000');
wctext=uicontrol('style','text','units','normalized','position',[0.52 0.17 0.13 0.065],'string','Center Frequency');

kk_=uicontrol('style','edit','units','normalized','position',[0.90 0.17 0.07 0.07],'string','2');
wctext=uicontrol('style','text','units','normalized','position',[0.76 0.17 0.13 0.065],'string','Order (even)');

half_=uicontrol('style','push','units','normalized','position',[0.04 0.25 0.21 0.07],'string','Half-Band Low-Pass','callback','tony_des_2(''start1'')');

tlow_=uicontrol('style','push','units','normalized','position',[0.28 0.25 0.21 0.07],'string','Tune Bandwidth','callback','tony_des_2(''start1'')');

tband_=uicontrol('style','push','units','normalized','position',[0.52 0.25 0.21 0.07],'string','Tune Center Freq','callback','tony_des_2(''start1'')');

fordr_=uicontrol('style','push','units','normalized','position',[0.76 0.25 0.21 0.07],'string','Polynomial Degree','callback','tony_des_2(''start1'')');

tprint=uicontrol('style','push','units','normalized','position',[0.30 0.1 0.65 0.05],'string','Print Coefficients','callback','tony_des_2(''start2'')');

readtext=uicontrol('style','text','units','normalized','position',[0.280 0.03 0.65 0.05],'string','(Highlight Parameter Value, Enter Desired Value, Depress Button)');
action='start1';
end

if strcmp(action,'start1')
coef=0;
end
if strcmp(action,'start2')
coef=1;
end

N=str2num(get(N_,'string'));
if rem(N,2)~=1
N=N+1;
N_=uicontrol('style','edit','units','normalized','position',[0.18 0.1 0.07 0.07],'string',N);
end
ws=str2num(get(ws_,'string'));
if ws<0.25
ws=0.300;
ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string',ws);
end
if ws>0.5
ws=0.300;
ws_=uicontrol('style','edit','units','normalized','position',[0.18 0.17 0.07 0.07],'string',ws);
end
wb=str2num(get(wb_,'string'));
if wb>0.5
wb=0.25;
wb_=uicontrol('style','edit','units','normalized','position',[0.42 0.17 0.07 0.07],'string','0.250');
end
wc=str2num(get(wc_,'string'));
if wc>0.5
wc=0.0;
wc_=uicontrol('style','edit','units','normalized','position',[0.66 0.17 0.07 0.07],'string','0.000');
end

kk=str2num(get(kk_,'string'));
if rem(kk,2)~=0;
kk=2;
kk_=uicontrol('style','edit','units','normalized','position',[0.90 0.17 0.07 0.07],'string','2');
end

if wb==0.25 & wc==0
mode=1;
else
mode=2;
end

aa=1;
if wc==0;
aa=-aa;
end

[rts1,rts2,coef0,coef1,ff1,ff2]=tonycxx(N,ws,wb,wc,kk,mode);

figure(2)
plot(roots(rts1),'x','linewidth',2,'markersize',8);
hold on;
plot(roots(rts2),'x','linewidth',2,'markersize',8);
plot(roots(conv(rts1,fliplr(rts2))-aa*conv(fliplr(rts1),rts2)),'o','linewidth',2)
plot(exp(j*2*pi*[0:.01:1]),'r');
plot([-1.1 1.1],[0 0],'k')
plot([0 0],[-1.1 1.1],'k')

hold off;
grid
axis([-1.2 1.2 -1.2 1.2]);
axis('square')
title('Roots of Two-Path Filter (Sum of Paths)')

figure(1)
plot(0:1/1024:.5-1/1024,20*log10(abs(ff1+0.0000001)));
hold
plot(0:1/1024:.5-1/1024,20*log10(abs(ff2+0.0000001)),'r--');
axis([0 .5 -100 10]);

hold

grid;
title('Magnitude Response of Two Path Filter');
xlabel('Normalized Frequency (f/fs)');
ylabel('Log Magnitude (dB)');

% write out coefficients
format long
if coef==1
disp(' ')
disp('Denominator Polynomials, path-0')
if wc==0
disp('   (a0 Z^2  +  a1 Z^1  +  a2)')
else
disp('   (a0 Z^4  +  a1 Z^3  +  a2 Z^2  +  a3 Z^1  +  a4)')
end
disp(' ')
disp(coef0)
disp(' ')

disp('Denominator Polynomials, path-1')
if wc==0
disp('    (b0 Z^2  +  b1 Z^1  +  b2)')
else
disp('    (b0 Z^4  +  b1 Z^3  +  b2 Z^2  +  b3 Z^1  +  b4)')
end
disp(' ')
disp(coef1)
end

[ - ]

Thanks!  Looks like it also needs the tonycxx() function.

[ - ]

function [rts1,rts2,coef0,coef1,ff1,ff2]=tonycxx(n_ord,wo,bw,fc,k_ordr, mode)
% mode = 1 for lowpass prototype
% mode = 2 to bandwidth and center frequency of prototype
% k_ordr for filter order
% function called by tony_des_2,
% computes paramters of two-path recursive all-pass filter
% also computes parameters when frequency transformed to two-path
% low-pass or band-pass filter
% Script file written by fred harris of SDSU. Copyright 2002

mode;
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_ord-1)/2);
aa=ww;
cc=[ww 0];
%step 2

for ii=1:(n_ord-1)/2

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

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

ordr1=floor((n_ord-1)/4);
ordr2=ordr1;
if n_ord-1-4*ordr1 ~= 0
ordr1=ordr1+1;
end

coef0=zeros(ordr1,3);
coef1=zeros(ordr2,3);

zz=zeros(1,k_ordr-1);

for ii=1:ordr1
den0(ii,:)=[1 zz cc(2*ii-1)];
end
coef0=den0;

for ii=1:ordr2
den1(ii,:)=[1 zz cc(2*ii)];
end
coef1=den1;

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

h1=[1];
for ii=1:ordr2
h1=conv(h1,den1(ii,:));
end

zz2=zeros(1,k_ordr/2);
h1=[h1 zz2];

g0=fliplr(h0);
g1=fliplr(h1);

rts1=h0;
rts2=h1;

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

ff1=0.5*(tp+bt);
ff2=0.5*(tp-bt);

if mode==1
return

else

tt=tan(pi*bw);
b=(1-tt)/(1+tt);

c=cos(2*pi*fc);

if fc==0

den0=zeros(ordr1,k_ordr+1);
den1=zeros(ordr2,k_ordr+1);
zz=zeros(1,(k_ordr/2)-1);
for nn=1:ordr1
c00=1+cc(2*nn-1)*b*b;
c01=-2*b*(1+cc(2*nn-1));
c01=c01/c00;
c02=cc(2*nn-1)+b*b;
c02=c02/c00;
den0(nn,:)=[1 zz c01 zz c02];
end
coef0=den0;

for nn=1:ordr2
c10=1+cc(2*nn)*b*b;
c11=-2*b*(1+cc(2*nn));
c11=c11/c10;
c12=cc(2*nn)+b*b;
c12=c12/c10;
den1(nn,:)=[1 zz c11 zz c12];
end
zz2=zeros(1,k_ordr/2-1);
den1(ordr2+1,1:1+k_ordr/2)=[1 zz2 -b];
coef1=den1;

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

h1=[1 zz2 -b];
for ii=1:ordr2
h1=conv(h1,den1(ii,:));
end

g0=fliplr(h0);
g1=fliplr(h1);
rts1=h0;
rts2=h1;

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

ff1=0.5*(tp+bt);
ff2=0.5*(tp-bt);

else

den0=zeros(ordr1,2*k_ordr+1);
den1=zeros(ordr2,2*k_ordr+1);

zz=zeros(1,(k_ordr/2)-1);

for nn=1:ordr1
c00=1+cc(2*nn-1)*b*b;
c01=-2*c*(1+b)*(1+cc(2*nn-1)*b);
c01=c01/c00;
c02=(1+cc(2*nn-1))*(c*c*(1+b*b)+2*b*(1+c*c));
c02=c02/c00;
c03=-2*c*(1+b)*(cc(2*nn-1)+b);
c03=c03/c00;
c04=cc(2*nn-1)+b*b;
c04=c04/c00;
den0(nn,:)=[1 zz c01 zz c02 zz c03 zz c04];
end
coef0=den0;

for nn=1:ordr2
c10=1+cc(2*nn)*b*b;
c11=-2*c*(1+b)*(1+cc(2*nn)*b);
c11=c11/c10;
c12=(1+cc(2*nn))*(c*c*(1+b*b)+2*b*(1+c*c));
c12=c12/c10;
c13=-2*c*(1+b)*(cc(2*nn)+b);
c13=c13/c10;
c14=cc(2*nn)+b*b;
c14=c14/c10;
den1(nn,:)=[1 zz c11 zz c12 zz c13 zz c14];
end
zz2=zeros(1,k_ordr/2-1);
den1(ordr2+1,1:1+k_ordr)=[1 zz2 -c*(1+b) zz2 b];
coef1=den1;

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

h1=[1 zz -c*(1+b) zz b];
for ii=1:ordr2
h1=conv(h1,den1(ii,:));
end

g0=fliplr(h0);
g1=fliplr(h1);
rts1=h0;
rts2=h1;

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

ff1=0.5*(tp-bt);
ff2=0.5*(tp+bt);
end

end

[ - ]