Long time lurker and first time poster - but unfortunately I haven't had any joy untangling this on my own.
I've been studying Mathias Lang's thesis, Algorithms for the Constrained Design of Digital Filters with Arbitrary Magnitude and Phase Response, with particular interest in the method for least squares design of stable IIR filters.
However, I'm struggling to understand the Jacobian calculation stage within the
Looking through the code, I'm comfortable with the calculations up to the section immediately prior, where the computation of
H and the complex response error is completed. After this, however, I fail to follow the Jacobian derivation.
I've included an excerpt as far as the line where my understanding runs out.
function [b,a,l2error] = mpiir_l2(M,N,om,D,W,r,a0) % MPIIR_L2: [b,a,l2error] = mpiir_l2(M,N,om,D,W,r,a0) % Least Squares Digital IIR Filter Design with Arbitrary Magnitude % and Phase Responses and Specified Maximum Pole Radius % % OUTPUT: % b numerator polynomial coefficients % a denominator polynomial coefficients % l2error approximation error % % INPUT: % M order of numerator polynomial % N order of denominator polynomial % om frequency grid, 0<=om<=pi % D complex desired frequency response on the grid om % W positive weighting function on the grid om % r maximum pole radius % a0 initial denominator guess; optional % % EXAMPLE: % Lowpass filter with approximately linear passband phase % (passband group delay = 19 samples, maximum pole radius = 0.97) % om=pi*[linspace(0,.2,20),linspace(.25,1,75)]; % D=[exp(-j*om(1:20)*19),zeros(1,75)]; % W=[ones(1,20),100*ones(1,75)]; % [b,a,e]=mpiir_l2(16,6,om,D,W,.97); % % Author: Mathias C. Lang, Vienna University of Technology, Oct. 98 % uses: lslevin, levin, update, locmax, % qp (Optimization Toolbox) om=om(:); D=D(:); W=W(:); srW = sqrt(W); EM = exp(-j*om*(0:max([M,N]))); tol = 1e-4; alpha = 0.5; MAXIT = 100; loopcnt = 0; fprintf('\n ERROR MAX. RADIUS STEP SLOPE\n'); fprintf('------------------------------------------\n'); % initial solution ini=0; if nargin == 7, a=a0(:); if length(a)~=N+1, ini=1; elseif a(1)==0, ini=1; elseif max(abs(roots(a))) > r, ini=1; else if a(1)~=1, a=a/a(1); end A = freqz(a,1,om); b = lslevin(M+1,om,A.*D,W./(abs(A).^2)); end else, ini=1; end if ini, % compute FIR solution a = [1;zeros(N,1)]; b = lslevin(M+1,om,D,W); end x = [a(2:N+1);b]; delta = ; % iterate (outer loop) while 1, % compute complex error, Jacobian, and objective function value A = EM(:,1:N+1)*a; B = EM(:,1:M+1)*b; H = B./A; E = srW.*(D - H); l2error = E'*E; if N<1, break; end vec1 = srW./A; vec2 = -vec1.*H; J = [vec2(:,ones(1,N)).*EM(:,2:N+1),vec1(:,ones(1,M+1)).*EM(:,1:M+1)];
Could some kind soul perhaps provide an explanation or an alternative expression of the final two lines? I can't quite wrap my head around this being the computation of a matrix of partial derivatives as it stands.
Having had a little play around with implementing a numerical approach to the Jacobian calculation, I think I can get that approach to work reasonably well. However, I still struggle to grasp the thesis implementation, which is presumably more direct/efficient.
My numerical approach, for reference:
% Compute Jacobian numerically. PERTURB_STEP = 1e-4; numericalJacobian = zeros(numel(H), nDenominator + 1 + mNumerator); for iXTerm = 1 : nDenominator + 1 + mNumerator % Perturb target x term. x1 = x; x2 = x; x1(iXTerm) = x(iXTerm) - PERTURB_STEP; x2(iXTerm) = x(iXTerm) + PERTURB_STEP; % Separate a/b coefficient sets. a1 = [1; x1(1:nDenominator)]; b1 = x1(nDenominator + 1 : nDenominator + 1 + mNumerator); a2 = [1; x2(1:nDenominator)]; b2 = x2(nDenominator + 1 : nDenominator + 1 + mNumerator); % Get H terms. A1 = zPowMinusN * a1; B1 = zPowMinusM * b1; H1 = B1 ./ A1; A2 = zPowMinusN * a2; B2 = zPowMinusM * b2; H2 = B2 ./ A2; % Fill jacobian. numericalJacobian(:, iXTerm) = (H2 - H1) / (2 * PERTURB_STEP); end % Weight Jacobian. numericalJacobian = numericalJacobian .* rootErrorWeight;