Moving the notches of a Moving Average Filter
Started by 4 months ago●9 replies●latest reply 4 months ago●209 viewsHi!
Is there a way of slightly moving all the notches of a recursive MAV filter without changing the sampling period?
The order of my filter is 16, and changing it dinamically to either 15 or 14 would make achive the goal though I'd like more finesse.
Thanks fellas!
In the moving average filter, the notches are at the frequency location in Hz that is the reciprocal of the length of the filter in seconds, independent of the sampling rate. This works out to be notches at normalized frequencies given at (n 2 pi)/N radians/sample where N is the total number of coefficients in the filter with n=1,2,...N-1. So with that you can modify the notches by changing the number of coefficients as you suggest. (Normalized frequency in radians per sample results in 2 pi being the sampling rate; this results in the unit circle on a complex plane representing the frequency axis, and as we'll see below will be convenient for going from poles and zeros on the complex "z-plane", as the roots of the transfer function to filter implementation).
You can design an Mth order filter (with N=M+1 coefficients) to place M notches at any frequency by setting the zeros to be at those frequencies specifically, where the frequency axis is the unit circle in the z-plane. For example consider a 4 sample moving average given by:
Writing the transfer function with positive exponents so that we can more easily get to its poles and zeros:
If we factor the numerator we get:
We see from this that it has 3 poles at the origin (which ensures causality) and three zeros at z = j, -1, -j.
Note in general the moving average filter will have zeros equally spaced around the unit circle divided evenly into angles given by n 2 pi/N except for z=1 (DC), which is exactly the description for the location of the nulls as I gave in the first paragraph! We can see this more directly if I rewrite H(z) so that the zeros are given in magnitude angle form (where here I didn't simplify to emphasize the n 2 pi/N pattern for the roots in the numerator:
So if we wanted to modify those nulls, we simply change those zeros to whatever frequencies we want and reverse the process to get H(z) back in negative powers of z for the direct implementation of the filter.
Hi Dan!
My issue here is that the filter goes in a FPGA and it is implemented recursively, like an interpolator (comb stage + integrator stage). This means I cannot change its coefficients directly. I can only access the value of the order, which must be an integer.
I seek to introduce minimal changes.
Maybe there is some DSP wizardy that can be used to represent a fractional delay using as inputs the signal and the signal delayed by one TS. The naive approach if I need a total delay of 15.5 samples would be averaging my 15th ans 16th samples, and that kind of works, but also kills the depth of the notches.
Thanks anyway for you prompt response!
Yes, you can replace z^(-1) blocks with fractional delay filters, but that doesn't fall into the category of "minimal changes", as you found to achieve the accuracy in the delay (and maintain the deep nulls) the approach to creating the fractional delay element would itself be a higher order filter (what you described was the first order approach). Perhaps someone else has a clever idea to offer.
One thing to consider is implementing the fractional delay filter in the comb itself, thus you could still do the CIC approach but that simple subtraction would replaced with a differencing of a higher order "all pass" or FIR approximation of an all pass as the fractional delay. This would scale all of your nulls proportionally and not allow individual control.
The simplest case is if you have a 16th order filter with zeros on the unit circle.
Thus Hz = 1 + z^(-16) or z^(16) + 1, if you prefer positiive powers.
This gives 16 zeros equally spaced on the unit circle.
If you want to shift all of them by the same amount exp(1i*pi/f), then just
multiply the last coefficient by exp(1i*N*pi/f), where N is the filter order.
This will shift all the zeros by exp(1i*pi/f). Note that the filter now has a
complex coefficient.
If the filter has non-zero coefficients for all the other powers of z, then the i^th
coefficient would be multiplied by exp(1i* (i-1)*pi/f) for all coefficients of
2 to N.
For example, let H(z) = z^4 -z^3 + .25z^2 + .25z -.125 has roots at
-0.5000 + 0i
0.5000 + 0.5000i
0.5000 - 0.5000i
0.5000 + 0i
Rotating the roots by exp(1i*pi/7), the roots would be at
-0.4505 - 0.2169i
0.6674 - 0.2335i
0.2335 + 0.6674i
0.4505 + 0.2169i
which corresponds to a new filter
Hs(z) = z^4 + (-0.9010 - 0.4339i)z^3 + (0.1559 + 0.1955i)z^2 + (0.0556 + 0.2437i)z + (0.0278 - 0.1219i)
These coefficeints can be obtained by multiplying the original coefficients of H(z) by exp(1i*k*pi/7) for k=0:4
times the original coefficient of H(z).
H(z) has coefficients [1 -1 .25 .25 -.125]
Thus (-0.9010 - 0.4339i) = exp(1i*1*pi/7) * -1
( 0.1559 + 0.1955i) = exp(1i*2*pi/7) * .25
( 0.2335 + 0.6674i) = exp(1i*3*pi/7) * .25
and ( 0.4505 + 0.2169i0 = exp(1i*4*pi/7) * -.125
The first coefficient is always 1 since exp(1i*0*pi/7) equals 1.
Here's a MATLAB script illustrating this.
h = [ 1.0000 -1.0000 0.2500 0.2500 -0.1250];
rh = roots(h);
rhs = rh * exp(1i*pi/7)
hs = conv([1,-rhs(1)],[1,-rhs(2)]);
hs = conv(hs,[1,-rhs(3)]);
hs = conv(hs,[1,-rhs(4)])
roots(hs)
% coefficients of new filter
c = h(2:5) .* exp(1i*(1:4)*pi/7)
Do you mean this structure:
If so just bypass one stage and it becomes 15, 2 stages and it becomes 14 and so on.
Hi G Medrano,
For order 16 MA filter very few FPGA resources would be used ... could you implement orders 14 and 15 and run concurrently ? Then you can use the different outputs as needed - or even implement some type of "consensus" algorithm.
-Jeff
Hi Engineer68-
That's not what I mean. What I'm saying is implement three (3) actual MA filters in parallel: order 16, 15, and 14. Together this is still very few FPGA resources unless were talking about a really dinky device.
Or as Kaz points out, implement order 16 but then tap off at stages 13 and 14 as needed -- assuming the OP's FPGA tools let him do that.
-Jeff