Forward/reverse median filter. Reference needed

Question

I am using the double median filtering of my data. For this, I perform median filtering of the signal in forward direction, then in reverse direction (from end to start). Results are averaged. Example of the MATLAB code:

Signal=rnd(1,100);
Mfw=medfilt1(Signal,10);
Mrev=medfilt1(Signal(end:-1:1),10);
M=mean(Mfw,Mrev(end:-1:1));

I am happy with this approach, but now it is time to write a report... I understand that I am not the first who use this kind of filtering.

Could you help me to find some references to this filter?

The problem is I even don't know the official name for this approach. "bidirectional median filter" shows me something else.

Answer 1

Globally, I do not know of a global name, because it is either a median or a (locally) linear filter with fluctuating tap locations: a median (odd length), a 3-tap Gaussian approximation (2-length) or a weighted average of at most 4 samples.

• when the length is odd ($2l+1$): a simple median filter
• when the length is two: a simple Pascal averaging linear filter with taps $[1\,,2\, ,1]/4$

• when the length is even ($2l$), greater than or equal to 4: a linear filter with up to 4 non-zero coefficients, but with different locations (depending on the order of the samples) across an $l+1$ window, with one of the following shapes:

  • $[\cdots,1\,\cdots,1\,\cdots,1\,\cdots, 1,\cdots]/4$
  • $[\cdots,1\,\cdots,2\,\cdots, 1,\cdots]/4$
  • $[\cdots,1\,\cdots,1\,\cdots]/2$

I could not have your initial code work on my computer. Let me offer a more compact version.

data = Signal'; 
lFilter = 10; 
dataDoubleMedianFilter = mean([medfilt1(data,length),flipud(medfilt1(flipud(data),length))],2);

Due to the symmetry in the median filter with odd length, your filter is just a simple median filter in that case. In other words, you average the two same median filtered signals.

If you set an even length, as you did, it is a little more involved. Morally, I do not understand the rationale. Even-length median "invents" novel values, stemming from the mean on the two centered samples.

Let us first try with a 2-median. On a centered set of values $[x_{-1}\,x_{0}\,x_{1}]$, one media will give $(x_{-1}+x_{0})/2$, the other $(x_{0}+x_{1})/2$, and their mean is the linear Pascal filter of length $3$, with taps $[1\,,2\,, 1]/4$:

$$(x_{-1}+2x_{0}+x_{1})/4$$

More generally, for an even number $2l$ of samples, the direct and the reverse frames share a center of $2l-1$ samples $x_1,\ldots,x_{2l-1}$, and differ from the 2 outer samples $x_0$ and $x_{2l}$. The resulting filter will yields the mean of the two most-centered samples from the left and the right frame. So it will be linear in aspect. However, the filter taps may have different locations. For instance, if the three-most centered samples (call them $x_{\sigma{(1)}},x_{\sigma{(2)}},x_{\sigma{(3)}}$ after reordering) of the $2l+1$ total frame belong to the $2l-1$ center frame, then they will be averaged with weights $1/4$, $1/2$, $1/4$:

$$1/4x_{\sigma{(1)}} + 1/2x_{\sigma{(2)}}+1/4x_{\sigma{(3)}}$$

If the outer samples reshuffle the order of the center frame, other combinations are obtained, with coefficients $1/2$, $1/2$ or $1/4$, $1/4$, $1/4$, $1/4$.

Answer 2

I've held off on answering this because I don't have MATLAB to test any of this.

I've employed the same trick with exponential smoothing. The purpose of this is to have the lags in the opposite direction cancel each other out, so the resultant smoothed signal is centered instead of shifted from the source signal. For sinusoids, the lag cancellation is exact.

So, I thought this might be the same motivation in this. It does not seem to be the case. (Side note: In my searching I found the MATLAB function "movmedian". I don't immediately see how it is that different from "medfilt1".)

I agree with Laurent Duval that if an odd number is used that the forward and backward filtered signals will be the same. However, OP is using an even number. It seems to me that the backward filtered signal will be just one index off of the forward one since the window size is the same and the median for each window should give the same result. Thus, excepting the end cases, you should be able to just do the forward one and then define:

$$ M[x] = ( Mfw[x] + Mfw[x+1] ) / 2 $$

Which is a very simple FIR.

Two cents minus inflation worth.

Ced