Forward and Backward Filtering (MATLAB's `filtfilt()`)

Question

I have an input sequence $x(n)$ and impulse sequence

\begin{equation} h(n) = \begin{cases} 0 & \quad \text{when } n = 0 \\ .239052 & \quad \text{when } n = 1 \\1.3357\exp(-0.650286n) & \quad \text{when } n > 1\\ \end{cases} \end{equation}

I did a convolution and obtained $y(n)=x(n)*h(n)$. I trimmed some of the data in the end from $y(n)$ and then I flipped $h(n)$ and I did reverse filtering $w(n)=y(n)*h(-n)$. When I researched about this I found filtfilt in matlab is an easier option but I am trying to implement this in C. So my question is without trimming the first output sequence $y(n)$ I am getting a drifted version of the result (there is a delay).

Is it trimming really necessary or am I doing something wrong here.
One more question is I can see some transients both at the start and end of $w(n)$ sequence when I use filtfilt I don't see that transients.The reason is matlab is using some initial states based on the below paper to get rid of the transients.

Gustafsson, F. "Determining the initial states in forward-backward filtering." IEEE® Transactions on Signal Processing. Vol. 44, April 1996, pp. 988–992.

To summarize :

1. Is it trimming really necessary ? PS the figure attached without trimming the the resultant data has a lag..Trimming the $y(n)$ sequence before backward filtering aligns the result.
2. how to determine the initial states for getting rid of transients.

Answer 1

My opinion is that the trimming isn't necessary.

The method in the paper is based on a simple idea, if we filter the signal $ x [n] $ forward and then backward, the result should be equal to the result of doing it backward and then forward.

The paper shows how to chose the initial conditions which minimize the $ {\ell}_{2} $ norm between the the different approaches.

Answer 2

For the first question, I would like to propose a very rudimentary method that can help you find out what is happening. One reason is that you do not provide too many details about your implementation. This method can even be used in practice when one cannot easily estimate filter delays.

Since your two-pass filter will be symmetric and monomodal (similar in shape to a two-sided exponential), you can estimate your "computational delay" (including all your windowing, padding, trimming, faulty allocations, subsampling...) by giving an input in the form of a Dirac-delta function, with its maximum 1 at a given location $k$, for instance far away from the edges. It could be like $$0\,0\,0\,0\,1\,0\,0\,0\,0$$ Then filter it your double-way, and find the index $l$ of the maximum in the output filtered signal. It could be longer that the original one: $$\epsilon_1\,\epsilon_2\,\epsilon_3\,\epsilon_4\,\epsilon_5\,\epsilon_6\,\epsilon_7\,\textrm{max}\,\epsilon_7\,\epsilon_6\,\epsilon_5\,\epsilon_4\,\epsilon_3\,\epsilon_2\,\epsilon_1\,$$ That $l-k$ (here $(8-5))$ gives you an integer quantity, a system delay, since your are using an LTI (linear, time-invariant) filter. Then, you can investigate where this delay comes from in your filtering process, and decide how to trim or (circularity) shift your "real signal".

For the second part, since your filter is partly exponential I would first suggest you to look at recursive implementations of exponential filters. You can split your filter $h$ as $h = h_1+h_2$ with $$h_1 = \{0\,0.239052\,0\,0\,0\, \dots\}$$ and $h_2$ the exponential part. The latter can be implemented with simple recursive formulas as exponential moving averages. Those can save you time, and treat borders with less pain.