# Filter design to realize Cauchy product

I come from Computer Science so please pardon for my possibly wrong terminology.

I need to design a filter which has coefficients $$h_0, h_1, \ldots, h_n, \ldots \quad\text{such that}\quad h_0 > h_1> \ldots > h_n > \ldots$$

The input is $x_n$. The output is $y_n$. The relationship between $h_n$, $x_n$ and $y_n$ is a convolution: $$y_n = h_0 x_n + h_1 x_{n-1} + \ldots + h_n x_0$$

The problem is that I have to carry out a convolution to compute the output $y_n$ upon the arrival of every new input $x_n$. So I have to keep all $x_1, x_2, \ldots, x_n, \ldots$ in memory and do $n+1$ multiplications and $n$ additions every time.

I would like to ask whether there is any way to design the coefficient $h_n$ such that I can compute $y_n$ more efficiently. I know of one solution, which is $h_n=\lambda^n$ with $0<\lambda<1$ but this function decays too fast.

Any idea to approximate this procedure is also welcomed.

• I am curious. What is your application? With $\lambda$ close to $1^-$, well-chosen wrt $n$, $\lambda^k$ may decay slowly – Laurent Duval Feb 28 '18 at 15:01

## 2 Answers

You could implement some convolutions as recursive filters. Consider moving average as an example, to calculate the output signal, instead of summing over the averaging window, you could add the new input and subtract the part of input which exited the window from the output signal to obtain the output of next time step. Also some filters could be implemented as lattice.

But there is not a general method to implement all filters recursively or as lattice.

• Thank for the answer. I notice that most of the time domain impulse responses of IIR filters decay exponentially. This is too fast for my application. Anyway, if I use $h_n=\frac{1}{n}$. Do you think that this design can be efficiently implemented? – Ngoc Anh Huynh Jun 3 '17 at 1:48
• A general method for reducing calculation of convolution is to transfer signal to Fourier space and multiply them, then return the answer to the time space again. Calculating convolutions has order of MN where M and N are length of signals but using FFT to implement convolutions the order of calculation will be Nlog(N). The double-sided (containing both negative and positive side) convolution with 1/n function is called Hilbert transform and its Fourier transform would be the sign function. Also there are lots of efficient implementation of Hilbert transform. – Mohammad M Jun 3 '17 at 16:59
• Take a look at the below link. mathematica.stackexchange.com/questions/341/… – Mohammad M Jun 3 '17 at 17:05
• Thanks for the pointers. The Hilbert transform is really enlightening. I understand that it is more efficient to do convolution in Frequency Space. But my other difficulty is that I need to produce $y_n$ incrementally from streaming input $x_n$. That said, if I am to use Fourier Transform or some sorts, I still need to keep all $x_1, x_2, ..., x_n,...$ in memory and have to do frequency transformation every time. Do you have any idea to address this issue? – Ngoc Anh Huynh Jun 4 '17 at 2:47
• Hi: you can stay in time domain if you make some assumptions about the impulse response. I won't go into it here but you'll see what I mean very clearly if you google for "koyck distributed lag". It will allow you to calculate the output in real time quite easily because it ends up reducing to exponentially smoothing with one extra coefficient. I'll see if I can find a link to a useful paper and send in another comment. – mark leeds Jan 29 '18 at 7:18

Apart from filters with exponential decay, and assuming you allowing approximations, the question resides in what does cost

• if your $$h_k$$ decay like exponentials, you can fit one or two exponentials to them, and there is a huge image processing literature on fast implementations of Gaussian/exponential filters.
• if float multiplication costs, you can use coefficients coded as dyadic rationals ($$i/2^n$$) with $$i$$ and $$n$$ integers, or sums thereof (SOPOT: sums of powers of two), with some multiplierless designs.
• if you can parallelize, one factorization idea is based on the polynomial Horner factorization: $$y_n = h_0 (x_n + h_1/h_0 (x_{n-1} + \ldots + (h_n /h_{n-1}\ldots h_0 x_0)))$$ where the $$h_k/\prod h_{k-1}$$ could be small (with your monotonic rule) or coded as SOPOTs. I don't know if this has been used elsewhere.