Re: About the derivation in firls.m

[je suis]

I've made a version of firls.m that handles even-length FIRs, too. I
am not versed in Octave, so the code will look awful, at best. The way
I did it uses ifft for the b vector, but it needs fairly large sizes
in order to become close to accurate. It works, though. Then, since I
still don't understand how the integral got to have cos(w)/(n^2*w)
from cos(w)/n^2, I just adapted the code as I saw fit. This also means
I can't go further and make differentiators and Hilbert transformers,
but I'll try the old hammer. It works, but it will hurt your eyes. It
would be nice for someone to verify it against Matlab, I don't have
it. The name of the function is arbitrary. Here it is:

function h = lsfir2(N,f,A,K);
% prepare a few helpers
M = floor(N/2);
m = [1:M+1];
oddN = mod(N, 2);
bands = length(f);

% calculate q...
q = zeros(1, 2*M + 1 + oddN);
for i=bands:-1:1
        q += (-1)^i*K(floor((i - 1)/2)+1).*f(i)*sinc(f(i)*[0:length(q)-1]);
end

% ...and the Q matrix
Q = toeplitz (q(m)) + hankel (q(m+oddN), q(m+M+oddN));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% I can't figure out cos_ints2 in firls.m so I'll just brute-force this.
% Resize the frequency/amplitude/weights vectors for proper ifft,
% since Selesnick derives the b[] vector from fp, rather than the usual
% (fp+fs)/2. So double the frequency values except the margins, and
% make K[] the same size as f[] with 0 inbetween (as the graph shows
% in Selesnick's paper). Then make A*K to be used with f in ifft.

% weights
K2 = zeros(1, 2*length(K)-1);
for i=1:2:bands-1
        K2(i) = K(floor(i/2)+1);
end
K2 = [K2; K2](:)';

% amplitude and frequency
A2 = zeros(1, 2*bands-2);
f2 = A2;
for i=1:2*bands-2
        f2(i) = f(floor(i/2)+1);
        A2(i) = A(floor(i/2)+1)*K2(i);
end

% prepare the magnitude, then ifft, and make the b vector
W = interp1(f2, A2, linspace(0,1,1024+1));
if oddN
        w = real(ifft([W zeros(1, 2*1024) W(end:-1:2)]));
        b = 2*w(2:2:2*M+2)';
else
        w = real(ifft([W W(end:-1:2)]));
        b = w(1:M+1)';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%
% This is adapted from firls.m, to include the even length filters.
% I don't understand how the integral formula from the comment became
% (cos(w2)-cos(w1))/(n*(w2-w1))+sinc, so I can't go further to make
% differentiators and Hilbert transformers, and general type III and IV FIRs.

w = f * pi;
i1 = 1:2:length (f);
i2 = 2:2:length (f);
if oddN
        cos_ints = sin((0.5:1:M+1)'*w);
else
        cos_ints = [w; sin((1:M)' * w)];
end
if ~oddN
        cos_ints2 = (w(i1).^2 - w(i2).^2)./(2*(w(i2) - w(i1)))
else
        cos_ints2 = [];
end
cos_ints2 = [cos_ints2; cos(((1:M+oddN)-.5*oddN)' * w(i2)) - ...
        cos(((1:M+oddN)-.5*oddN)' * w(i1)) ./ ...
        (([1:M+oddN]-.5*oddN)' * (w(i2) - w(i1)))];
d = [-K .* A(i1); K .* A(i2)] (:);
b = [];
if ~oddN
        b = [1, 1./(1:M)]';
else
        b = [1./((1:M+1)-.5*oddN)]';
end
b = b.* ((kron (cos_ints2, [1, 1]) + cos_ints(1:M+1,:)) * d)/pi;

%%%%%%%%%%%%%%

% the rest of the algorithm
a = Q \ b;
if(oddN)
        h = [flipud(a); a];
else
        h = [ a(end:-1:2); 2*a(1); a(2:end) ];
end