tDFiltFIR (LogFile)

function tDFiltFIR (LogFile)

% $Id: tDFiltFIR.m,v 1.8 2009/06/17 15:25:04 pkabal Exp $

addpath('..');

if (nargin == 0)
  LogFile = [];
end
StopDiaryFile = StartDiaryFile(LogFile, 'tDFiltFIR.ref');

fprintf('==================================\n');
fprintf('BPF1: N=24 lowpass filter\n');
clear B
[B(1:2).Freq]  = deal([0 0.08], [0.16 0.5]);
[B(1:2).Value] = deal(1, 0);
h = DFiltFIR(24, B, 'bpf');
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF2: N=32 bandpass filter\n');
clear B
[B(1:3).Freq]   = deal([0 0.1], [0.2 0.35], [0.425 0.5]);
[B(1:3).Value]  = deal(0, 1, 0);
[B(1:3).Weight] = deal(10, 1, 10);
h = DFiltFIR(32, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF3: N=50 bandpass filter, unequal weighting\n');
clear B
[B(1:3).Freq]   = deal([0 0.15], [0.2 0.3], [0.35 0.5]);
[B(1:3).Value]  = deal(0, 1, 0);
[B(1:3).Weight] = deal(1/0.1, 1/1, 1/0.01);   % 1 / deviation
h = DFiltFIR(50, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF4: N=31 bandstop filter\n');
clear B
[B(1:3).Freq]   = deal([0 0.1], [0.15 0.35], [0.42 0.5]);
[B(1:3).Value]  = deal(2, 0, 2);
[B(1:3).Weight] = deal(1, 50, 1);
h = DFiltFIR(31, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF5: N=55 multiband filter\n');
clear B
[B(1:5).Freq]   = deal([0 0.05], [0.1 0.15], [0.18 0.25], [0.3 0.36], ...
                       [0.41 0.5]);
[B(1:5).Value]  = deal(0, 1, 0, 1, 0);
[B(1:5).Weight] = deal([20 10], 1, 3, 1, [20 30]);
[B(1:5).WeightInt] = deal ('log', [], 'linear', [], 'log');
h = DFiltFIR(55, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF6: N=32 constrained bandpass filter (c.f. BPF2)\n');
clear B
[B(1:3).Freq]   = deal([0 0.1], [0.2 0.35], [0.425 0.5]);
[B(1:3).Value]  = deal(0, 1, 0);
[B(1:3).Weight] = deal(10, 1, 10);
B(1).LLimit = 0;
B(1).ULimit = 1;
h = DFiltFIR(32, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF7: N=37 constrained bandstop filter\n');
clear B
[B(1:5).Freq]   = deal([0 0.1], [0.15 0.1995], [0.2 0.3], ...
                       [0.3005 0.35], [0.4 0.5]);
[B(1:5).Value]  = deal(1, 0, 0, 0, 1);
[B(1:5).LLimit] = deal(0.9, 0,    0,    0,    0.98);
[B(1:5).ULimit] = deal(1,   0.03, 0.01, 0.03, 1.02);
h = DFiltFIR(37, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF8: N=76 constrained multiband filter\n');
clear B
[B(1:7).Freq]   = deal([0 0.05], [0.0505 0.0995], [0.1 0.15], ...
                       [0.18 0.20 0.22 0.23 0.24 0.25], [0.2505 0.2995], ...
                       [0.3 0.36], [0.41 0.5]);
[B(1:7).Value]  = deal(0, 0, 1, 0, 0, 1, 0);
[B(1:7).Weight] = deal(200, 0.05, 1, 60, 0.05, 1, 400);
[B(1:7).LLimit] = deal(0, 0, 0.92, [0  0 -Inf -Inf 0 0], 0, -Inf, 0);
[B(1:7).ULimit] = deal(0.0004, 1, 1, ...
                       [0.0013333 0.0013333 Inf Inf 0.0013333 0.0013333], ...
                       1, Inf, 0.0002);
h = DFiltFIR(76, B);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF9: N=24 lowpass filter (Density 10, c.f. BPF2)\n');
clear B
[B(1:2).Freq]  = deal([0 0.08], [0.16 0.5]);
[B(1:2).Value] = deal(1, 0);
h = DFiltFIR(24, B, 'bpf', 1, 10);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF10: N=24 lowpass filter (High grid density, c.f. BPF2)\n');
clear B
[B(1:2).Freq]  = deal([0 0.08], [0.16 0.5]);
[B(1:2).Value] = deal(1, 0);
h = DFiltFIR(24, B, 'bpf', 1, 100);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF11: N=33 highpass filter (zero dc response)\n');
clear B
[B(1:2).Freq]  = deal(0, [300 4000]);
[B(1:2).Value] = deal(0, 1);
B(1).LLimit = 0;
B(1).ULimit = 0;
h = DFiltFIR(33, B, 'bpf', 8000);
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF12: N=999 filter (highly oversampled)\n');
clear B
[B(1:2).Freq]   = deal([0 3200], [4800 256000]);
[B(1:2).Value]  = deal(1, 0);
[B(1:2).Weight] = deal(1, 60);
h = DFiltFIR(999, B, 'bpf', 512000);
PrMatrix(h');

fprintf('==================================\n');
fprintf('DIF1: N=32 differentiator\n');
clear B
B.Freq  = [0 0.5];
B.Value = 1;
h = DFiltFIR(32, B, 'dif');
PrMatrix(h');

fprintf('==================================\n');
fprintf('HIL1: N=20 Hilbert transformer\n');
clear B
B.Freq  = [0.05 0.5];
B.Value = 1;
h = DFiltFIR(20, B, 'Hil');
PrMatrix(h');

fprintf('==================================\n');
fprintf('BPF13: log deviations\n');
InvdBV = inline('10.^(x/20)');

Value = [1 0];
dBdev = [0.2 -40];
dev(1) = Value(1) .* (InvdBV(2 * dBdev(1)) - 1) ./ ...
                     (InvdBV(2 * dBdev(1)) + 1);
dev(2) = InvdBV(dBdev(2));
Weight = 1 ./ dev;
Weight = Weight / min(Weight);
clear B
[B(1:2).Freq]   = deal([0 0.2], [0.3 0.5]);
[B(1:2).Weight] = deal(Weight(1), Weight(2));
[B(1:2).Value]  = deal(Value(1), Value(2));
h = DFiltFIR(19, B);
PrMatrix(h');

% ----------
StopDiaryFile();

return