Transfer function of recursive moving average filter

[Slev1n]

Hey guys,

I have the attached moving average filter. Here the picture of the blockdiagramm:

What I understand: output = (sum of all values) - (sum of all values M-points before); M = #points used for averaging which is equal to the control "decimate".

1, So is H(z)=(1+z^-1+z^-2+...+z^-n-1)^stages?

2, Is y[n]=y[n-1]+x[n]-x[n-M]? Or is it some kind of other recursive ma algorithm?

Kind regards

Slev1n

[WCR]

Slev1n,

I could not load the VI, because I do not have the .xnode files or some of the sub VIs. Thus my comments are based on the text of your message. You use "n" in your equations but you use M in your description.  I assume M=n.

The fact that M is called "decimate" makes me wonder.  If this VI is doing a "moving average decimation" (which I cannot tell), then (a) that's a bad thing to do from a sigal processing point of view* and (b) it is a pretty opaque and complicated VI for doing that simple task.

First: single pass transfer function.  Then: multi-pass.

1. The transfer function H1(z)=Y1(z)/X(z) for a SINGLE-pass moving average of n points:

y1(k) = (x(k) + x(k-1) + x(k-2) + ... + x(k-n+1)) / n  (time domain)      {1}

therefore Y1(z) = (X(z) + X(z)/z + X(z)/z^2 + ... + X(z)/z^n-1) / n (z transform of above)

therefore Y1(z) = X(z) * (1 + 1/z + 1/z^2 + ... + 1/z^n-1) / n      (algebra)

therefore H1(z) = Y1(z)/X(z) = (1 + z^-1 + z^-2 + ... + z^(-n+1)) / n    {2}

2. Multi-pass, with P passes:

If you really cascade the above filter m times, then the overall transfer function will be the mth power of the single-pass transfer function:

Y(z) = [H1(z)]^m = (1/n)^m * (1 + z^-1 + z^-2 + ... + z^(-n+1))^P     {3}

which is similar but not identical to your equation 1.  It is a little tricky to implement the above scheme correctly.  The trickiness has to do with the fact that the output y updates "instantaneously" with the input. Thus the output value y(k) reflects all P passes of filtering of [input up through and including x(k)].

Now on to your second equation: When the moving average filter width, n, is large, the simple-minded version of the single-pass MA filter (described in eq.{1} above) is not as efficient as a recursive implementation.  The recursive version is

y1(k) = y1(k-1) + x(k)/n - x(k-n)/n   {4}

It is not obvious that filters {1} and {4} would have the same transfer function. To find transfer funciton for {4}, note that

y1(k) - y1(k-1) = (x(k) -x(k-n)) / n

therefore

Y1(z) * (1 - z^-1) = X(z) * (1 - z^-n) / n

therefore

H1(z) = Y1(z)/X(z) = [(1 - z^-n)  / (1 - z^-1)] / n   {5}

The equivalence of {2} and {5} can be demonstrated by using the formula for a geometric series of finite length.

*If you are decimating, then, to avoid aliasing, you must low-pass-filter the original signal to remove all frequencies above the decimated signal's Nyquist frequency, before decimating.  A moving average that takes n points at a time does not do it, even if the MA filter is applied sequentially to al the points before decimation.

[WCR]

Oops. I introduced P, instead of m, for the number of passes, when I realized you used M for the decimation factor.  I did not change all my m's to Ps. I should have written:

If you really cascade the above filter P times, then the overall transfer function will be the Pth power of the single-pass transfer function:

Y(z) = [H1(z)]^P = (1/n)^P * (1 + z^-1 + z^-2 + ... + z^(-n+1))^P     {3}

[Slev1n]

Hey WCR,

thank you for your detailed answer.

This VI is not from myself but I found it here: https://decibel.ni.com/content/docs/DOC-1762 LIA_IPnet_LabVIEW2009.zip

Here some points I think I have to mention:

- As I am decimating with a CIC filter before this MA I think it is not really a decimation, the author just named it like this.

- I only make one addition and one substraction, so for me it has to be a recursive MA but I can't get to the equation y(n)=y(n-1)+x(n)-x(n-M); (M=decimation) from the code itself

- If I understand you correctly, you say that the frequency response of simple and recursive MA are the same. That's what I found out so far, too.

- I also want to mention, that it works fine for my application (Lock in amplifier) so far, only if I increase the number of stages/cascading above 4, I get issues as the values get unstable

- The division by M, the amount of points used for averaging, is done on the host!

Here you have the program part of the other case. If you want you can rebuild it from the pictures (BRAM has 4096 cells and I64 format). Actually there are no SubVis inside this code. In the zip folder you can find a MA.vi which should be the same!

So my questions are:

1, How to get from this code to the recursive equation (I think it is one)  y(n)=y(n-1)+x(n)-x(n-M)

2, Is the frequency response the following: |H(w)| = [sin(w*M/2)/M*sin(w/2)]^P

3, The unstabelness beyond 4 stages, although stages*decimation is less than 4096 ( I think this point can be discussed later)

[4 stages and 1024=decimation works, but 8 stages and 512=decimation works not. My guess so far is, that I64 is not big enough and therefor it cannot be stored correctly in the BRAM]

Thank you for help so far WCR!!

kind regards

SLev1n

[WCR]

SLev1n,

 Q1: I can't prove or disprove that the code yields the recursive difference equation which you include.  When I try loading the moving average VI from the zip file you included, I had 16 instances of subVIs, xnodes, and controls which Labview could not find.  The documentaiton you included says that the code requires the FPGA package and several hardware items, none of which I have.

You mention that you decimate with a CIC filter prior to this MA filter.  That sounds surprising and I wonder if this IS the CIC filter, and maybe decimation too. Typically one uses a CIC filter (which is a recursive MA filter) to remove high frequencies from the input signal, then one deicmates.

Q2. Yes the equatoin you gave is the theoretical magnitude portion of the frequency response for a moving average filter with P passes, IF you add parentheses as shown, to clarify what goes together in the denominator:

|H(w)| = [ sin(w*M/2) / (M*sin(w/2)) ] ^P

I added () around M*sin(w/2).

Q3. I don't know why it would be unstable with more than 4 passes.  I have certainly seen digital filters become unstable as the filter order increases. But those were floating point filters with delicately balanced coefficients.  This is a fixed point filter with no multiplications. I think your guess is plausile that it is due to hardware limits.

[Slev1n]

Hey WCR,

I just noticed, that I didnt mention that this is a FPGA code...sorry for that.

I've been talking to one of the guys, who is responsible for the "decibel" page I got the code from, but he is not the engineer who build the code, so his support is limited. He told me it is a recursive moving average.

I am not sure if you are into lock in technique, but the main CIC purpose is downsampling, although the lowpass characteristic is a positive side effect and as I only want the DC value I dont care about bad passband characteristics.

As the signal has a lot of noise which is significantly bigger than the DC signal itsself (SNR: -50dB) I use a MA for filtering the noise, too. Later I want to compare simple moving average to exponentially weighted average.

I still dont think that the MA is decimating the signal, because decimation would cause downsampling and as after every loop iteration a value is at at the output (single pass mode) the signal is not sampled down, only averaged.

That's how I understand the code so far:

• If the #stage=1 than we have single pass and only one loop iteration per input. If stage=2 the value is passed  two times through the filter and so on...

• The value "decimation" is responsible for the adresses used in the BRAM and the amount of input values used for averaging. So if it is 10 and stage=1, only 10 cells are used and averaging is performed over 10 values.

• The first adder is accumulating the input, overflow is no problem if "wrap" is used, the added values are stored in the BRAM.

• If decimation=10, the input value which is accumulated at the first adder and stored in the BRAM, will stay in the BRAM for 10 iterations and then will be substracted from the most recent accumulated value. y[n]=x[n]+x[n*]-x[n*-10];  

x[n*] is the value of the accumulator/first adder;

x[n] is the most recent input value, so x[n]+x[n*] is the "most recent accumulated value"

x[n*-10] is the "accumulated value 10 inputs before"

As the code is working for my application (locking a -54dB SNR signal) it has to be a moving average due to the noise suppression. As there is only one adder and one subtractor I also think it is recursive. I just wanted to understand it and derive the recursive equation 🙂

Maybe someone who has FPGA package can help, too?

Kind regards

Slev1n

[Slev1n]

Maybe someone can also tell me the bit growth of a recursive moving average filter. As a CIC is such a filter, can I really just use this equation?

Bmax=Bin-1+N*log2(decimate*D); (N=number of integrator and comb stages, which should be 1 here and D=differential delay which is also 1 here)

[Slev1n]

I think I found the crux!

It is the position of the first feedback node!

The first picture is my MA the second the CIC. In Both cases the feedback nodes are, in my opinion, at the wrong position. If I change the position as implied, the filter do not work anymore but than I could derivethe eqution: y[n]=y[n-1]+x[n]-x[n-M] and therefor proof that it is a recursive moving average.

Can anyone help me?

I really dont get it...

kind regards

Slev1n

[Slev1n]

Here the non fpga version. so everybody should be able to test it.

[WCR]

Thanks Slev1n.  I can run it now that the FPGA stuff is gone, but I don't get it. The input and output of the posted VI are scalars.  It takes one integer and outputs one integer. Is that a mistake?  This is supposed to be a filter for a 1 dimensional signal.  Therefore I expect the input and output to be 1D arrays.