1
$\begingroup$

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?

$\endgroup$
1
$\begingroup$

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

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Floris SA Oct 31 '18 at 11:51
  • $\begingroup$ Never mind that last point, I see from the impulse response that the length is indeed 2k+1. Thanks again! $\endgroup$ – Floris SA Oct 31 '18 at 12:03
  • $\begingroup$ @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. $\endgroup$ – Fat32 Oct 31 '18 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.