Polyphase decimation filter plot in Octave

Question

I'm using de Soras hiir library for decimation and I'm trying to verify the filter design in Octave. The filter is a polyphase 2x downsampling filter with 96dB attenuation and 0.01 transition band, computed like:

PolyphaseIir2Designer::compute_coefs(c2x, 96.0f, 0.01f);

The result is:

c2x[12] = { 0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 };

Now, according to the code, the polyphase filters transfer function is: $$ H(z) ={ \frac{1}{2}(\prod_{k=0}^{K/2-1} \frac{a_{2k+1}+z^{-2}}{1+a_{2k+1}z^{-2}} +z^{-1}\prod_{k=0}^{(K-1)/2} \frac{a_{2k}+z^{-2}}{1+a_{2k}z^{-2}} } ) $$

Which I would like to verify by using the coefficients to draw the frequency response for each added allpass section, until reaching the final response. This is the program I wrote in Octave:

pkg load signal;
pkg load ltfat;
close all, clear all
%c = [0.10717745346023573,0.53091435354504557];
%c = [0.041893991997656171, 0.16890348243995201, 0.39056077292116603, 0.74389574826847926];
c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
H=0; 
i=0; 
%subplot(211)
while (i<K)
  beta = c(i+1);
  disp(beta);
  b0=[beta 0 1];
  a0=[1 0 beta];

  beta = c(i+2);
  b1=[0 beta 0 1];  % note, extra zero!
  a1=[1 0 beta];

  [A0,w]=freqz(b0,a0);
  [A1,w]=freqz(b1,a1);

  A = 0.5*(A0 + A1); % Computing the IIR halfband filter frequency response

  if (H!=0)
    H = H.*A;
  else
    H=A;
  endif
  plot(w/pi, 20*log10(abs(H))); % amplitude plot in decibel  
  hold on;
  ylabel("Gain (dB)"), axis([0 1 -100 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end

And this is the output:

Which looks nothing like what I would expect, as it doesn't meet the requirements Any idea what could be wrong? I had to swap the coefficients between the two allpass sections compared to the math to make it look even half way correct), but perhaps the allpass filters are the problem or the cascading? Thanks.

---

So I've looked more into this, but still havn't solved it.

After looking at the actual code, it seems that these coefficients are not for second order filters after all. The two parallel filters are both of the form:

b*(x(n) - y(n-1)) + x(n-1) => b*x(n) + 1*x(n-1) - b*y(n-1)

which means that these are 1st order allpass with: b0=b, b1=1, a1=b and transfer function:

$$ H(z)= \frac{b+z^{-1}}{1+bz^{-1}} $$

The input to the second allpass chain is one-sample delayed. (Although, I'm a little confused as the input seems to be at the source-rate, rather than the target-rate.) Still, I've put this into Octave/Matlab code, using a different approach where I construct the parallel filters and run the passes on the impulse response, and analyze the result:

c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
H=0; 
i=0; 

y = zeros(1,1024);
y(1)=1;

while (i<K)
  beta0 = c(i+1);
  b0=[beta0 1];
  a0=[1 beta0];

  beta1= c(i+2);
  b1=[0 beta1 1]; 
  a1=[1 beta1];

  y0 = filter(b0, a0, y);
  y1 = filter(b1, a1, y);

  y = 0.5*(y0 + y1); 

  [A,w]=freqz(y);

  H=A;

  plot(w/pi, 20*log10(abs(H))); % amplitude plot in decibel  
  hold on;
  ylabel("Gain (dB)")%, axis([0 1 -100 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end

And this is the result:

Which still doesn't look anything like what I would expect of a 2x halfband decimation filter?? If anyone has any idea of what could be wrong, I would appreciate it.

Answer 1

I think I might have finally solved it myself.

Afaik, while it is true that the actual allpass filters are 1st order, the trick is to see that the original formula actually just says that the two chains of allpasses run interleaved. That is the reason for the z^-2 and the reason for the extra multiply by z^-1. Hence it is perfectly reasonable to simply use z^-1 filters in the implementation, and feed the allpasses the data interleaved.

However, for plotting the response, I had to first compute the cascaded filter coefficients, and then finally apply the z^-1 delay to the A1 chain.

This is the result:

Which is exactly what I wanted. This is the final code:

c = [0.036681502163648017, 0.13654762463195794, 0.27463175937945444, 0.42313861743656711, 0.56109869787919531, 0.67754004997416184, 0.76974183386322703, 0.83988962484963892, 0.89226081800387902, 0.9315419599631839, 0.96209454837808417, 0.98781637073289585 ];
K = length(c);
i=0; 

b0=[1];
a0=[1];
b1=[1];
a1=[1];

while (i<K)

  beta0 = c(i+1);
  b0=conv(b0,[beta0,0,1]);
  a0=conv(a0,[1,0,beta0]);

  beta1= c(i+2);
  b1=conv(b1,[beta1,0,1]); 
  a1=conv(a1,[1,0,beta1]);

  b1z1 = [0,b1];

  [A0,w] = freqz(b0,a0,512,2); 
  [A1,w] = freqz(b1z1,a1,512,2);

  A=(A0+A1)/2;

  plot(w, 20*log10(abs(A))); % amplitude plot in decibel  
  hold on;
  ylabel("Gain (dB)"), axis([0 1 -120 20])
  xlabel("Normalized frequency");
  grid on

  i+=2;
end

I think this should be the accepted answer, but if anyone has anything to add, please do.