function leastSquaresFIR()
    close all;
    
    % **************************************************
    % construct filter requirements
    % only relative weights between different frequency
    % regions matter for the LS solver. 
    % Here, they are for convenience written in terms of 
    % _absolute_ requirements (i.e. x dB ripple, y dB stopband
    % rejection).
    % The solution from the IRLS algorithm will either pass 
    % or fail them all simultaneously by a similar margin.
    % **************************************************    
    om = linspace(0, 1, 20000);

    if true
        % passband: 0 dB gain, 0.2 dB ripple
        omega1 = om((om >= 0) & (om < 0.3));
        [HTarget1, weight1] = passbandSpecHw(0, 0.2, numel(omega1));
        
        % first stopband: 100 dB attenuation
        omega2 = om(om >= 0.5 & om < 0.6);
        [HTarget2, weight2] = stopbandSpecHw(100, numel(omega2));

        % second stopband: 80 dB attenuation
        omega3 = om(om >= 0.7 & om < 0.8);
        [HTarget3, weight3] = stopbandSpecHw(80, numel(omega3));
    else
        % first passband: 0 dB gain, 0.2 dB ripple
        omega1 = om((om >= 0) & (om < 0.3));
        [HTarget1, weight1] = passbandSpecHw(0, 0.2, numel(omega1));
        
        % second passband: -60 dB gain, 0.2 dB ripple
        omega2 = om(om >= 0.5 & om < 0.6);
        [HTarget2, weight2] = passbandSpecHw(-60, 0.2, numel(omega2));
        
        % stopband: 80 dB attenuation
        omega3 = om(om >= 0.7 & om < 0.8);
        [HTarget3, weight3] = stopbandSpecHw(80, numel(omega3));
    end

    % nominal group delay (usually half the FIR length)
    delay_samples = 15;
    % -pi..pi covers negative and positive frequencies
    omega = [omega1, omega2, omega3] * pi;
    HTarget = [HTarget1, HTarget2, HTarget3];
    
    % frequency-dependent phase shift from nominal group delay
    HTarget = HTarget .* exp(-1i * omega * delay_samples);    
    weight = [weight1, weight2, weight3];
    
    % **************************************************
    % example coefficient routing
    % **************************************************    
    var = 1;
    switch var
      case 1
        % ordinary FIR: One coefficient per tap
        map =  eye(31);
      case 2
        % symmetric FIR: one coefficient handles two delay taps
        map = zeros(16, 31);
        map(sub2ind(size(map), 1:16, 1:16)) = 1;
        map(sub2ind(size(map), 15:-1:1, 17:31)) = 1;
      case 3
        % FIR with selected taps shared (second and second-last row)
        X = 1;
        map = [
            X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X
            0 X X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X X 0
            0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0
            0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0
            0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0
            0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0
            0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0 0 0 0 X X 0 X X 0 0 0 0 0 0 0 0 0 0 0 0 0
            0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 X 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
              ];
    end
    % dlmwrite('map.txt', map, 'delimiter', ' ');
    
    obj = struct(...
        'iter', 150, ...
        'tol', 1e-6, ...
        'complex', false, ...
        'omega', omega, ...
        'HTarget', HTarget, ...        
        'weight', weight, ...
        'map', map); 
    
    % run IRLS solver
    h = run_IRLS(obj);

    doPlot(h, obj.omega, obj.HTarget, obj.weight);
end

% **************************************************
% IRLS algorithm
% **************************************************
function h = run_IRLS(obj)
    
% force column vectors
    obj.omega = obj.omega(:);
    obj.HTarget = obj.HTarget(:);
    obj.weight = obj.weight(:);
    
    % phase shift for one sample delay in the frequency domain
    zInv = exp(-1i * obj.omega);
    
    weight = obj.weight .^ 2;
    h = [];
    for ixIter = 1:obj.iter
        hPrev = h;
        h = run_LS(obj.omega, obj.HTarget, weight, obj.complex, obj.map);
        
        % evaluate frequency response
        % the first vector element in polyval() scales the highest
        % power of z^-1 => use fliplr() on h
        HDesign = polyval(fliplr(h), zInv);
        
        % status report
        ixIter
        % the weighted least-squares error that is minimized by the iteration
        evalErr = max(abs((HDesign - obj.HTarget) .* obj.weight))

        % determine error to nominal freq. response
        err = abs(HDesign - obj.HTarget);
        
        % check convergence
        if ~isempty(hPrev)
            deltaH = norm(h - hPrev);
            if deltaH < obj.tol
                disp('change in impulse response below tolerance - exiting');
                break;
            end
        end
        
        % iteratively adjust weight
        weightedError = obj.weight .* err;
        f = weightedError / (weight' * weightedError);
        weight = weight .* f;
    end
    if ixIter >= obj.iter
        disp('iteration limit reached');
    end
    
    % plot weight
    figure(34); grid on; hold on; 
    logw = 20*log10(weight) .';
    logw = logw - max(logw); % Align max. weight with 0 dB line
    plot(obj.omega/2/pi, logw, 'r+');
end

% **************************************************
% least-squares algorithm
% **************************************************
function h = run_LS(omega, HTarget, weight, isComplexFlag, map)
    
    [nCoeff, nDelay] = size(map);
    nEquations = nCoeff;
    if false && nCoeff == nDelay && norm(map - eye(size(map))) < eps        
        % map is square identity matrix
        disp('using Toeplitz matrix structure');
        map = []; % flag use of faster Toeplitz method (special case)
    end
    
    outSample = [];
    inSample = [];

    % force column vectors
    omega = omega(:); 
    HTarget = HTarget(:); 
    weight = weight(:); 

    % waveCoeffIn: 
    %     Magnitude and phase coefficients of an input signal 
    %     to the unknown filter. 
    %     The solver will find the filter that, applied to 
    %     samples of this signal, causes the smallest possible
    %     error (least-squares).
    %     The time domain waveform at a time t could be
    %     evaluated (sampled) as 
    %     waveIn(t) = waveCoeffIn * exp(1i * omega * t/tSample);
    %     At t = 0, this simplifies to
    %     waveIn(t) = sum(waveCoeffIn), as exp(0)==1.
    %
    %     For evaluation at t = k tSample with k = 0, 1, 2, ... :
    %     waveIn(k tSample) = sum(waveCoeffIn .* exp(1i * omega) .^ k)
    waveCoeffIn = weight;
    
    %     multiplication with z^-1 (zInv) causes a phase shift that 
    %     is equivalent to a unity delay. 
    %     multiplication with z does the opposite, it results in a later
    %     sample
    z = exp(1i*omega);
    zInv = exp(-1i*omega);
     
    % waveCoeffOut: 
    %     An ideal (zero-error) filter operating on the signal described by
    %     waveCoeffIn would result in the signal described by waveCoeffOut.
    waveCoeffOut = weight .* HTarget;
    
    % build one input signal sample for each element of the delay line
    for ixDelay = 1:nDelay+(nEquations-nCoeff) 
        % evaluate time domain waveform at t=0 
        % => all exp() expressions are 1
        inSample(ixDelay, 1) = sum(waveCoeffIn);

        % inSample represents the FIR delay line, and higher-indexed
        % entries are older samples => use zInv
        waveCoeffIn = waveCoeffIn .* zInv;
    end

    % build one output signal sample for each row of the equation system
    for ix = 1:nEquations
        % every row in the equation system corresponds to the input signal, 
        % at a later time, with the FIR chain shifted. 
        % an ideal filter would transform inSample into outSample.
        outSample(ix, 1) = sum(waveCoeffOut);
        % multiply with z, as we move towards more recent samples
        waveCoeffOut = waveCoeffOut .* z;
    end

    % normalize average power (not mandatory)
    scale = 1 / sqrt(weight .' * weight);    
    inSample = inSample * scale;
    outSample = outSample * scale; 
    
    if ~isComplexFlag
        % a real-valued design could mirror omega and H around 0 Hz, resulting
        % in real-valued coefficients.
        % Alternatively, remove the imaginary parts of the rotating phasors.        
        % note: real(exp(ix)) = cos(x) = 1/2 (exp(ix) + exp(-ix))
        inSample = real(inSample);
        outSample = real(outSample);
    end
    
    if isempty(map)
        M = toeplitz(inSample);
    else
        % each row in M is the signal in the FIR delay chain, as it appears at the
        % input of the coefficient multiplier, creating one value in outSample.
        % The first row should match outSample(1) etc.
        % Coefficients are then optimized to match all values as closely as possible (LMS).
        ixMap = find(map);
        nEntries = numel(ixMap);
        
        % one index pair for each delay <=> coefficient assignment
        [ixCoeff, ixDelay] = ind2sub(size(map), ixMap);
        
        % repeat, as the map will be applied to each row in M individually
        ixMap = repmat(ixMap, nEquations, 1);
        ixCoeff = repmat(ixCoeff, nEquations, 1);
        ixDelay = repmat(ixDelay, nEquations, 1);
        ixRow = repmat(1:nEquations, nEntries, 1);
        ixRow = ixRow(:);
        
        % The lower rows are "snapshots" of the FIR delay chain taken at
        % a later time. Therefore, increasingly newer input samples 
        % contribute.
        ixInSample = ixDelay - (ixRow-1);
        
        % inSample has been evaluated only for one FIR length into the past, relative 
        % to the first row. Evaluating the other side would be tedious; an 
        % alternative is to exploit the same symmetry as explained in the "Toeplitz" section.
        ixPast = find(ixInSample >= 1);
        ixFuture = find(ixInSample < 1);
        inSampleLookup = zeros(numel(ixInSample), 1);
        inSampleLookup(ixPast) = inSample(ixInSample(ixPast));
        inSampleLookup(ixFuture) = conj(inSample(2-ixInSample(ixFuture)));
        
        % at this point:
        % input waveform value "inSampleLookup" appears at coefficient "ixCoeff"
        % after shifting the delay line (ixRow-1) times.
        % When multipliers are shared (i.e. symmetric FIR), there will be multiple
        % contributing values to each coefficient with different delays.
        % The value in map gives the relative contribution. 
        % For example, map(3, 4) = 0.5; map(3, 6) = 1 routes output from delay tap 
        % 4 and 6 through coefficient 3 to the output, the former scaled down by 
        % 0.5 (maybe using an efficient bit shift).
        
        M = zeros(nEquations, nCoeff);
        for ix = 1:numel(ixMap)
            M(ixRow(ix), ixCoeff(ix)) = M(ixRow(ix), ixCoeff(ix)) + map(ixMap(ix)) * inSampleLookup(ix);
        end
    end
 
    if isempty(map)
        disp('replace with algorithm that exploits Hermitian Toeplitz structure');
        % see the last line in program lslevin on page 20,
        % http://www.nt.tuwien.ac.at/fileadmin/users/gerhard/diss_Lang.pdf
    end

    % Using pinv as it allows modifications that include more
    % equations than coefficients
    hPre = pinv(M) * outSample;
    
    % map solver coefficients to FIR coefficients
    if isempty(map)
        h = hPre;
    else
        % convert to unity-delay-based impulse response
        h = zeros(1, nDelay);
        for ix = find(map) .'
            [ixCoeff, ixDelay] = ind2sub(size(map), ix);
            h(ixDelay) = h(ixDelay) + hPre(ixCoeff) * map(ix);
        end
    end
end

function doPlot(h, omega, HTarget, weight)
    % **************************************************
    % impulse response
    % **************************************************    
    figure(1); hold on;
    if norm(imag(h)) > eps
        stem(real(h), 'k');
        stem(imag(h), 'b');
        legend('real(h)', 'imag(h)');
    else
        handle = stem(h);
        set (handle, 'lineWidth', 3);
        %legend('h');
        handle = plot(h, 'o');
        set (handle, 'lineWidth', 3);
        set (handle, 'markerSize', 2);
    end
    
    % **************************************************
    % frequency response
    % **************************************************    
    figure(34); hold on;
    % zero padding for higher frequency resolution
    h = [h(:) .', zeros(1, 10000)];
    
    % FFT frequency basis
    n = numel(h);
    fb = 0:(n - 1);
    fb = fb + floor(n / 2);
    fb = mod(fb, n);
    fb = fb - floor(n / 2);
    fb = fb / n; % now [0..0.5[, [-0.5..0[

    handle = plot(fftshift(fb), fftshift(20*log10(abs(fft(h)) + 1e-15)));
    set (handle, 'lineWidth', 3);
    xlabel('normalized frequency');
    ylabel('|H(f)|');

    if false
        % zoom
        xlim([0.25, 0.41]);
        ylim([-110, -55])
    elseif true        
        % full range
        xlim([0, 0.5]);
        ylim([-180, 1]);
    else
        % include negative frequencies for complex-valued
        % filter
        xlim([-0.5, 0.5]);
    end
    legend('weights', 'frequency response');
end

function [H, w] = passbandSpecHw(gain_dB, ripple_dB, n)
    H = ones(1, n) * 10 .^ (gain_dB / 20);
    err = abs(H) * (1 - 10 ^ (-ripple_dB / 20));
    w = 1 ./ err;
end

function [H, w] = stopbandSpecHw(attenuation_dB, n)
    H = zeros(1, n);
    err = ones(1, n) * 10 ^ (-attenuation_dB / 20);
    w = 1 ./ err;
end
