Least squares FIR for type II

Question

Following this paper (see the Matlab example at the end), I am trying to make a least-squares algorithm in Octave, but for type II (I know about firls). For type I FIRs, it works flawlessly, but not type II. This is my attempt (in Octave):

M = floor(N/2);
oddN = 1 - mod(N, 2);
n = [1:2*M + (~oddN)] - 0.5*(~oddN);
m = [0:M];
if(oddN)
    q = [wp+K*(1-ws), wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
else
    q = [wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
end
Q = (toeplitz(q(m+1)) + hankel(q(m+1), q(m+M+1)));
b = wp*sinc(wp.*(m+0.5*(~oddN)))';
a = Q\b;
if(oddN)
    h = [a(M+1:-1:2); 2*a(1); a(2:M+1)];
else
    h = [a(M+1:-1:2); a(1); a(1); a(2:M+1)];
end

Could someone please point out the mistakes?

---

I found one of the mistakes: q derives from the integral whose second part is simplified as -K*ws*sinc(ws), but the 0.5 offset in sampling does not discard the K*sinc(k+0.5) term. That was just a bad omission from my part. So I added it:

n = [~oddN : 2*M] + (~oddN)/2;
q = wp*sinc(wp.*n) + K*(sinc(n) - ws*sinc(ws.*n));

The results changed, but they are still wrong. Even so, the a(1) terms in the middle seem to need multiplying by a constant. I used the hammer a bit and ~2.25 seems to be a reasonable value. It doesn't change the fact that the results are wrong (or that the hammer is not a solution).

Answer 1

I may have fond the answer. The $Q$ matrix, made from the sum of the Hankel and Toeplitz forms, stems from:

$$cos(k\omega)cos(n\omega) = \frac{1}{2}\left[cos((k-n)\omega)+cos((k+n)\omega)\right]$$

but with the additional $0.5$ for even orders, it comes out as:

$$cos((k+0.5)\omega)cos((n+0.5)\omega)=\frac{1}{2}\left[cos((k-n)\omega)+cos((k+n+1)\omega)\right]$$

and since the $k+n$ term makes the Hankel matrix, it needs an increment in its indices, and the $q$ vector longer by 1. The $b$ vector only needs the extra $0.5$. Now the relevant code lines look like this:

q = [wp + K*(1 - ws), wp*sinc(wp.*n) - K*ws*sinc(ws.*n)];
Q = (toeplitz(q(m + 1)) + hankel(q(m + 1 + oddN), q(m + M + 1 + oddN)));
b = wp*sinc(wp.*(m + 0.5*(~oddN)))';

and this is a comparison between a 30th and a 31st orders, with f=[0 0.3 0.4 1], A=[1 1 0 0], K=[1 10]:

What I don't understand now is why the $q$ vector is not made with the $0.5$ sampling offset.