Example: Filters and banded matrices

This example illustrates the filters used in the algorithms for LPF/TVD and LPF/CSD. The filter is a non-causal zero-phase filter defined as H = inv(A)*B where A and B are banded matrices. The filter is implemented efficiently using fast solvers for banded systems.

Ivan Selesnick,
selesi@poly.edu,
NYU-Poly

Contents

Start

clear
close all

printme = @(filename) print('-dpdf', sprintf('figures/filter_Example_%s', filename));

Define high-pass filter H

N = 100;                    % N : signal length

fc = 0.05;                  % fc : cut-off frequency in cycles/sample (0 < fc < 0.5)

d = 2;                      % d : degree of filter is 2d (d must be small positive integer: 1, 2, 3, etc.)

[A, B, B1, D, a, b, b1] = ABfilt(d, fc, N);  % Construct banded matrices (sparse data type in MATLAB)

H = @(x) A\(B*x);           % H : function for filter H

Verify banded system solver

We verify that MATLAB calls a banded system solver to implement H.

spparms('spumoni', 3);      % Set sparse monitor flag to obtain diagnostic output

x = rand(N, 1);
H(x);

%  Diagnostic output:
%
%  sp\: bandwidth = 2+1+2.
%  sp\: is A diagonal? no.
%  sp\: is band density (1.00) > bandden (0.50) to try banded solver? yes.
%  sp\: is LAPACK's banded solver successful? yes.

spparms('spumoni', 0);      % Reset sparse monitor flag to no diagnostic output
sp\: bandwidth = 2+1+2.
sp\: is A diagonal? no.
sp\: is band density (1.00) > bandden (0.50) to try banded solver? yes.
sp\: is LAPACK's banded solver successful? yes.

Display high-pass filter H

Frequency response of high-pass filter H

[Hf, om] = freqz(b, a);     % Compute frequency response

figure(1)
clf
plot(om/pi/2, abs(Hf))
line(fc, 0.5, 'marker', 'o')
ylim([0 1.2])
title('HPF frequency response')
xlabel('Frequency (cycles/sample)')

printme('hpf_freqresp')

Impulse response of high-pass filter H

no = 40;
x = zeros(N,1); x(no) = 1;      % input signal
y = H(x);                       % output signal

figure(1)
clf
stem((1+d):N-d, y, '.')
title(sprintf('HPF impulse response (n_o = %d)', no))
xlabel('Time (samples)')

printme('hpf_impresp')

Pole-zero diagram of high-pass filter H

figure(2)
clf
zplane(b, a);
title('HPF pole-zero diagram')

printme('hpf_zplane')

Null space of H

Compute the null space of the high-pass filter H Result: null space is 4-dimensional as expected. The null space consists of polynomials of order 3.

u = null(full(A\B));
size(u)

figure(2)
clf
plot(u)
title('High-pass filter null space')
ans =

   100     4

Define low-pass filter L

Define low-pass fitler L as I - H. Note: because H reduces the length of the signal, I must also reduce the length of the signal.

L = @(x) x - [nan(d, 1); H(x); nan(d, 1)];

Impulse response of low-pass filter L

no = 40;
x = zeros(N,1); x(no) = 1;      % input signal
y = L(x);                       % output signal

figure(1)
clf
stem(1:N, y, '.')
title(sprintf('LPF impulse response (n_o = %d)', no))
xlabel('Time (samples)')

printme('lpf_impresp')

Frequency response of low-pass filter L

b_lpf = a - b;
a_lpf = a;

[Hf_lpf, om] = freqz(b_lpf, a_lpf);

figure(1)
clf
plot(om/pi/2, abs(Hf_lpf))
line(fc, 0.5, 'marker', 'o')
ylim([0 1.2])
title('LPF frequency response')
xlabel('Frequency (cycles/sample)')

printme('lpf_freqresp')

Pole-zero diagram of low-pass filter L

figure(2)
clf
[hz, hp, hl] = zplane(b_lpf, a_lpf);
title('LPF pole-zero diagram')

printme('lpf_zplane')

Polynomial preservation property of L

The low-pass filter L exactly preserves polynomial signals of order 3 (2d-1), i.e. x(n) = c0 + c1 n + c2 n^2 + c3 n^3.

n = (1:N)';
x = rand + rand * n + rand * n.^2 + rand * n.^3;
err = x - L(x);
max(abs(err))
% ~zero!
ans =

   2.0955e-09

Polynomial signals of higer order are not preserved:

x = rand * n.^4;
err = x - L(x);
max(abs(err))
% large value!
ans =

   2.3739e+03

LPF example

Apply LPF to some noisy data.

x = cos(0.05*pi*n) + 0.2*randn(N,1);        % x : a noisy data
y = L(x);                                   % y : low-pass filtered data

figure(1)
clf
plot(1:N, y, 'b',  1:N, x,'k.')

% Note: Sometimes, the filter does not work so well at the start and
% end of the signal. The output of the filter is too high or
% too low. This is due to transient effects.
% These start- and end-transients can be addressed by
% appropriate polynomial extrapolation of the data
% prior to filtering.