Software & Algorithms

[moodz]

Noise removal from the samples is a core requirement for the DSP processing ....

Typical Averaging algorithms.

Box Car.
Smoothing algorithms involving unweighted averages can be
performed in various ways. One method, known as boxcar averaging, involves collecting a
specified number of points and averaging them to produce a single point.

Moving Average. Another method
involves a "moving" average, where a specified number of successive points (n) are collected
and averaged, then the next measurement is averaged with the previous n-1 measurements, and
this process continues through the data set. Uses alot or RAM if more than 1 sample point is to be averaged.

Moving Average Box. ( my name for it ... cant find a citation )
The Box Car and Moving average both either throw away information and/or utilise excessive amounts of ram and cpu cycles. In this case the problem is to extract the averaged waveform from noise and without distortion of the waveform features at high speed.
Funnily enough I think this algorithm had its origins in financial software for extracting market trends from huge amounts of raw stock data in real time but I cannot lay my hand on the reference ....
The algorithm goes thus ....

CONDITIONS / ASSUMPTIONS
A periodic waveform is sampled at a specific time for S samples and the requirement is to average the wave form for N complete sets of S samples ( ie waveforms) .
The sampling MUST commence at exactly the same time ( ie synchronous) with the waveform
N should be a binary number 4,8,16,32,64,128,256,512,1024 etc etc for optimal code performance. The higher the number the more noise will be removed from the waveform but it will take longer to for each cycle to complete.
This will require N locations of ram for storage noting that these locations will also be used for summing ( dependant on bit size of ADC multiplied by averaging N bit size ).
The initial value of the ram locations (eg after CPU reset / initialisation ) is arbitrary ... though for this case we will set it to 0.

EXECUTION
1. Sample each point along the waveform from 1 to S and ADD the ADC value to corresponding RAM location ( 0 to S-1). ie keep a running sum in each RAM location.

2. Increment the Average count and repeat step 1 till N cycles have been done.

3. The MSB ( most significant bits + oversampling gain bits ) of each RAM location now contain the average value ... grab or process the results of the averaging NOW. ( ie you will do this every N cycles of waveform conversions)

4. Now REINITIALISE each ram location by dividing the current value by 2 ...( most efficient to shift the value 1 bit right.)

5.Set the average count N to 0 and GOTO step 1.

END.....

For example if I use a long integer ( 4 bytes ) in DSPIC to store 512 samples ( ie I am digitising the RX waveform in a PI circuit for 512 us ) then I need 2048 bytes of RAM + a couple of counter variables for index pointer and cycle counter.
Utilising this algorithm I dont have to remember the discrete values as for example with the traditional moving average I have to pop out the oldest value and insert the new for averaging and I would have to do this for each point in the waveform ( I would need N x S x no_of_bytes_to_store or 1024 x 512 x 2 bytes !!!)

The dspic code for the diff front end averages 64 x 64bit discrete samples 5300 times a second with an averaging cycle of 1024 for each of the 64 bit sample points on the waveform. I have run it out to 512 samples but this is stretching the friendship if the pulse repetition rate is to remain high ( > 4 Khz ).

On an FPGA with 24 bit ADCs the results using this algorithm are spectacular !

moodz.

[scuba840]

4. Now REINITIALISE each ram location by dividing the current value by 2 ...( most efficient to shift the value 1 bit right.)

Great post. I am going to try it. However I am a little confused. On step #4 (above), why divide by 2 to reinitialize?

Can you help me out here?

[moodz]

Great post. I am going to try it. However I am a little confused. On step #4 (above), why divide by 2 to reinitialize?

Can you help me out here?

Ok .... you could calculate a new average for each sample point by initialising to zero however this would not remove low frequency noise less than twice the sample run frequency. ( nyquist sampling ).
Also it kind of makes sense that unless you reduce by some amount each sample buffer would overflow. By removing half we are shifting the sum right one bit and thus each sample will be the result of the current sample + 1/2 the previous sample + 1/4 the one before that + 1/8 the one before that etc etc.. so you may be able to see that the current average at any point depends on the n previous samples to some degree. Probably a bad explanation but it works like a charm. Try it and see......

Regards,

Moodz.

[scuba840]

4. Now REINITIALISE each ram location by dividing the current value by 2 ...( most efficient to shift the value 1 bit right.)

5.Set the average count N to 0 and GOTO step 1.

Now I understand a little more.
Another question...If the purpose of step #4 (above) is to retain some past history (the 1.2, 1/4, 1/8, etc) data, then why do we set AVERAGE=0 in step 5 (above)? It is almost like each RAM location will contain the remnants of at lease 1 sample, so why not set AVERAGE=1 in step 5???

[moodz]

Now I understand a little more.
Another question...If the purpose of step #4 (above) is to retain some past history (the 1.2, 1/4, 1/8, etc) data, then why do we set AVERAGE=0 in step 5 (above)? It is almost like each RAM location will contain the remnants of at lease 1 sample, so why not set AVERAGE=1 in step 5???

Hi ... in step 5 we are resetting a counter not the average value ... will post a spreadsheet example of the code. Moodz

[scuba840]

Now I understand a little more.
Another question...If the purpose of step #4 (above) is to retain some past history (the 1.2, 1/4, 1/8, etc) data, then why do we set AVERAGE=0 in step 5 (above)? It is almost like each RAM location will contain the remnants of at lease 1 sample, so why not set AVERAGE=1 in step 5???

I look forward to the spreadsheet! That would be fantastic.
I need to correct my question.....

Another question...If the purpose of step #4 (above) is to retain some past history (the 1.2, 1/4, 1/8, etc) data, then why do we set COUNT=0 in step 5 (above)? It is almost like each RAM location will contain the remnants of at lease 1 sample, so why not set COUNT=1 in step 5???

[moodz]

I look forward to the spreadsheet! That would be fantastic.
I need to correct my question.....

Another question...If the purpose of step #4 (above) is to retain some past history (the 1.2, 1/4, 1/8, etc) data, then why do we set COUNT=0 in step 5 (above)? It is almost like each RAM location will contain the remnants of at lease 1 sample, so why not set COUNT=1 in step 5???

Ahh...now I think I see why you are asking ... the code explanation was based on a binary implementation. It is important that the sample counts be a power of two if you are storing results as binary numbers.... you could use sample counts as power of 10 if you use floating point ... like in a pentium processor for example ...however for PIC type CPUs the binary implementation is most efficient. If you are working in binary ... counts start at 0 not 1. Thus if you have an 8 bit count for example to include all possible 256 counts the count must run from 0 to 255. Otherwise you would need a 9 bit counter to count from 1 to 256 which would work but not strictly compatible with binary logic.
Each RAM location does contain 1/2 the previous "history" of the last average cycle + 1/4 + 1/8 + 1/16 etc of the ones before that at the start of each summing cycle.
If you implement the code it is quite easy to compare the difference between initialising each ram location with half the previous sum vs 0 and you should see that the noise reduction is much improved with the half residual averaging method. ( also if you are say storing 16 bits worth of sum on each cycle and you add half the previous sum then the sum division number is 17 bits not 16 bits as it would be if the average was initialised to 0 .. so you get 1 bit extra of averaging resolution for the half sum residual method )

moodz.
Regards,

moodz.

[maca404]

I have heard that refered to as Ensemble Average, I have just startedwith fpga by interfacing an existing hammer head circuit. What I am interested in is how you intend to extract the information from the wave form. I am particularly interesting in gold and rejecting the bad ground we have here. I was thinking of selecting certain periods on the decay curve and calculating the rate of change, then comparing the rate of change at eg. 1us to 10us to 50us If they are all the same its something with a long decay or ground and you could keep a running average of the ground signal for rejection, Then if there is a sudden spike in ROC at around 10us you could have a gold signal mixed in there. I realize it is oversimplified but thats all I can think to do.

[simonbaker]

Noise removal from the samples is a core requirement for the DSP processing ....

Typical Averaging algorithms.

Box Car.
Smoothing algorithms involving unweighted averages can be
performed in various ways. One method, known as boxcar averaging, involves collecting a
specified number of points and averaging them to produce a single point.

Moving Average. Another method
involves a "moving" average, where a specified number of successive points (n) are collected
and averaged, then the next measurement is averaged with the previous n-1 measurements, and
this process continues through the data set. Uses alot or RAM if more than 1 sample point is to be averaged.

Moving Average Box. ( my name for it ... cant find a citation )
The Box Car and Moving average both either throw away information and/or utilise excessive amounts of ram and cpu cycles. In this case the problem is to extract the averaged waveform from noise and without distortion of the waveform features at high speed.
Funnily enough I think this algorithm had its origins in financial software for extracting market trends from huge amounts of raw stock data in real time but I cannot lay my hand on the reference ....
The algorithm goes thus ....

CONDITIONS / ASSUMPTIONS
A periodic waveform is sampled at a specific time for S samples and the requirement is to average the wave form for N complete sets of S samples ( ie waveforms) .
The sampling MUST commence at exactly the same time ( ie synchronous) with the waveform
N should be a binary number 4,8,16,32,64,128,256,512,1024 etc etc for optimal code performance. The higher the number the more noise will be removed from the waveform but it will take longer to for each cycle to complete.
This will require N locations of ram for storage noting that these locations will also be used for summing ( dependant on bit size of ADC multiplied by averaging N bit size ).
The initial value of the ram locations (eg after CPU reset / initialisation ) is arbitrary ... though for this case we will set it to 0.

EXECUTION
1. Sample each point along the waveform from 1 to S and ADD the ADC value to corresponding RAM location ( 0 to S-1). ie keep a running sum in each RAM location.

2. Increment the Average count and repeat step 1 till N cycles have been done.

3. The MSB ( most significant bits + oversampling gain bits ) of each RAM location now contain the average value ... grab or process the results of the averaging NOW. ( ie you will do this every N cycles of waveform conversions)

4. Now REINITIALISE each ram location by dividing the current value by 2 ...( most efficient to shift the value 1 bit right.)

5.Set the average count N to 0 and GOTO step 1.

END.....

For example if I use a long integer ( 4 bytes ) in DSPIC to store 512 samples ( ie I am digitising the RX waveform in a PI circuit for 512 us ) then I need 2048 bytes of RAM + a couple of counter variables for index pointer and cycle counter.
Utilising this algorithm I dont have to remember the discrete values as for example with the traditional moving average I have to pop out the oldest value and insert the new for averaging and I would have to do this for each point in the waveform ( I would need N x S x no_of_bytes_to_store or 1024 x 512 x 2 bytes !!!)

The dspic code for the diff front end averages 64 x 64bit discrete samples 5300 times a second with an averaging cycle of 1024 for each of the 64 bit sample points on the waveform. I have run it out to 512 samples but this is stretching the friendship if the pulse repetition rate is to remain high ( > 4 Khz ).

On an FPGA with 24 bit ADCs the results using this algorithm are spectacular !

moodz.

Hi moodz:

Where do the good old Markov processes fit into this concept?

For example, the simple filter:

S(t) = A * S(t - 1) + B * x(t)

where A + B = 1 and are coefficients. Typically A is much larger than B.

S is the "running average"
x is a sampled signal at sample time (t).
I don't know how to do subscripts, so I used (t -1) to mean the sample just prior to sample (t).

I think this is very close to a first-order low-pass filter. There are other equivalents for higher order filters.

-SB

[moodz]

Hi moodz:

Where do the good old Markov processes fit into this concept?

For example, the simple filter:

S(t) = A * S(t - 1) + B * x(t)

where A + B = 1 and are coefficients. Typically A is much larger than B.

S is the "running average"
x is a sampled signal at sample time (t).
I don't know how to do subscripts, so I used (t -1) to mean the sample just prior to sample (t).

I think this is very close to a first-order low-pass filter. There are other equivalents for higher order filters.

-SB

Hi Simon,
The aim of the filtering is to recover the flyback + target waveform - noise.
If a wideband PI signal from 10 kHz to 1Mhz is filtered through any sort of low / high / bandpass filter then some of the wanted signal is attenuated along with the noise. The algorithm I described recovers the waveform minus the noise. It does this because it is a synchronous filter. The difference between the wanted signal and the unwanted noise is that the noise is not synchronised to the sample rate. Thus the time term for the synchronous equation is removed. The performance of the filter is only dependant on the number of samples. An infinite number of samples would remove all the noise. ( not possible of course ). An oscillator is a type of synchronous filter. In most oscillators the wanted signal is much higher than the noise.

Regards,

moodz.

[simonbaker]

Hi Simon,
The aim of the filtering is to recover the flyback + target waveform - noise.
If a wideband PI signal from 10 kHz to 1Mhz is filtered through any sort of low / high / bandpass filter then some of the wanted signal is attenuated along with the noise. The algorithm I described recovers the waveform minus the noise. It does this because it is a synchronous filter. The difference between the wanted signal and the unwanted noise is that the noise is not synchronised to the sample rate. Thus the time term for the synchronous equation is removed. The performance of the filter is only dependant on the number of samples. An infinite number of samples would remove all the noise. ( not possible of course ). An oscillator is a type of synchronous filter. In most oscillators the wanted signal is much higher than the noise.

Regards,

moodz.

So if you look at each sample point (which are treated the same), you have:

Y(i) = .5 * Y(i - n) + Average(S(i-n) + S(i-n+1)... S(i))

Where Y(i) is your filtered value at sample time i, and S(i) is a sample at time i.

So you kind of take an average of n samples, then use a Markov process to filter the average with the previous average.

But something looks a little off with the coefficients of the Markov process -- they should add to 1 I think for an unbiased estimate. More like this:

Y(i) = .5 * Y(i - n) + .5 * Average(S(i-n) + S(i-n+1)... S(i))

Based on your algorithm, you would need to average each sample divided by two.

-SB

[moodz]

So if you look at each sample point (which are treated the same), you have:

Y(i) = .5 * Y(i - n) + Average(S(i-n) + S(i-n+1)... S(i))

Where Y(i) is your filtered value at sample time i, and S(i) is a sample at time i.

So you kind of take an average of n samples, then use a Markov process to filter the average with the previous average.

But something looks a little off with the coefficients of the Markov process -- they should add to 1 I think for an unbiased estimate. More like this:

Y(i) = .5 * Y(i - n) + .5 * Average(S(i-n) + S(i-n+1)... S(i))

Based on your algorithm, you would need to average each sample divided by two.

-SB

I think you may be throwing data away if you halve before summing ?

moodz.

[simonbaker]

I think you may be throwing data away if you halve before summing ?

moodz.

I think, looking again, your algorithm is equivalent to mine if you consider the "reset" (after dividing the sum by two) values as your signal estimate. I was considering the point after the sum as the signal estimate -- just a different timing practicality.

I agree, implementation-wise I wouldn't want to divide each sample by two before averaging, I was just trying to fit the algorithm into my picture of the signal estimate. But again, if the signal estimate is the final "reset" value, then it looks right to me.

It is kind of a funky algorithm though -- not really a moving average at all. Well, it seems to be a moving average of non-overlapping "chunk" averages. I'm sure fine for practical purposes, but I wonder what the advantage is over a simple digital filter. I still may not have a clear understanding of the details though. And if N = 4, hardly an issue.

Regards,

-SB

[simonbaker]

I think, looking again, your algorithm is equivalent to mine if you consider the "reset" (after dividing the sum by two) values as your signal estimate. I was considering the point after the sum as the signal estimate -- just a different timing practicality.

I agree, implementation-wise I wouldn't want to divide each sample by two before averaging, I was just trying to fit the algorithm into my picture of the signal estimate. But again, if the signal estimate is the final "reset" value, then it looks right to me.

It is kind of a funky algorithm though -- not really a moving average at all. Well, it seems to be a moving average of non-overlapping "chunk" averages. I'm sure fine for practical purposes, but I wonder what the advantage is over a simple digital filter. I still may not have a clear understanding of the details though. And if N = 4, hardly an issue.

Regards,

-SB

to be a moving average of non-overlapping "chunk" averages.

Correction: I guess I should say a simple first-order digital filter of non-overlapping "chunk" averages, if you see what I mean.

-SB