How can I implement a triangular filter in MatLab, given it's Z-transform representation?

Question

I have to implement an anti-aliasing filter for a certain processing step in MatLab. While searching in literature for inspiration, I came accross a paper in which the authors wrote the following:

We implement the lowpass filters as triangular smoothing. The Z-transform representation of an arbitrary N-point filter of this type is:

$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}$$

where the triangle length $N=2k+1$, and the filter scaling coefficient $\alpha=(k+1)^2$.

I was thinking about how to implement this filter, so following the approach taken in this DSP SE post, I rewrote the expression for the filter as follows:

$$G(z)=\frac{-z^{-(k+1)}+2-z^{k+1}}{\alpha(1-z)(1-z^{-1})}=\frac{Y(z)}{X(z)}$$ $$X(z)\Bigl(-z^{-(k+1)}+2-z^{k+1}\Bigr) = Y(z)\Bigl(\alpha(1-z)(1-z^{-1})\Bigr)$$ let's assume for now we use a filter of length $N=3$, which gives us $k=1$, we get: \begin{align} X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(1-z)(1-z^{-1})\Bigr)\\ X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(4(2-z-z^{-1})\Bigr)\\ X(z)\Bigl(-z^{-2}+2-z^{2}\Bigr) & = Y(z)\Bigl(8-4z-4z^{-1})\Bigr)\\ -x[n-2]+2x[n]-x[n+2] &= 8y[n]-4y[n+1]-4y[n-1]\\ \end{align} which I finally rewrite as: $$y[n+1]=\Bigl(8y[n]-4y[n-1]-2x[n]+x[n-2]+x[n+2]\Bigr)/4$$

I implemented this result in matlab as follows:

n = 1:1000;
x = sin(0.10*n) + (-0.5+rand(1,1000));
out=zeros(size(x));

for n_idx = 3:998
    out(n_idx+1) = (8*out(n_idx) - 4*out(n_idx-1) - 2*x(n_idx) + x(n_idx-2) + x(n_idx+2) )/4;
end

subplot(2,1,1); plot(n,x);
subplot(2,1,2); plot(n,out);

but unfortunately, the output looks absolutely nothing like a lowpass filtered version of the input. The authors do not go in detail about how they implement this filter, the only clue they give is this line:

we note that the denominator of $g(z)$ represents a causal integration of the data with $1/(1-z)$, followed by an anticausal integration with $1/(1-z^{-1})$.

which I can't make sense of.

My questions are therefore: -1 Is my implementation of the filter $G(z)$ correct? -2 Would anyone know what the authors mean in their comment on (anti)causal integration?

Answer 1

This MATLAB/OCTAVE code makes a causal implemtation the system $G(z)$ you mention. Note that $G(z)$ is non-causal, so you have to assume delays to make it causal. You can convert it to non-causal if you have the full data available before processing.

k = 5;
N = 2*k+1;

K = (k+1)^(-2);

b = K*[1, zeros(1,k) , -2, zeros(1,k), 1];
a = [1,-2,1];


[Hk,w]=freqz(b,a,1024*8,'whole');
figure,plot(w/pi,20*log10(abs(Hk)));title('frequency response');
ylabel('Magnitude (dB)');
xlabel('Frequency ( x \pi) ');

hn = impz(b,a,N);
figure,stem(0:N-1,hn);title('(causal) impulse response h[n]');