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.