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

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?

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]');

• Thanks for your reply, this was very helpful. One small question: the filter length is defined as $N=2k+1$, but the length of b is actually $2(k+1)+1$, or $2k+3$. Do you think this is a mistake by the authors, or is filter length something different from what I think it is? – Floris Oct 31 '18 at 11:51
• Never mind that last point, I see from the impulse response that the length is indeed 2k+1. Thanks again! – Floris Oct 31 '18 at 12:03
• @FlorisSA sorry for late reply. Yes that's right, the filter length is the impulse response length, whereas the order refers to LCDDE representation; the maximum delay. For simple, causal, FIR filters length = order +1. For IIR filters length is infinite. – Fat32 Oct 31 '18 at 19:14