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