After a lot of tedious code stepping, I'm fairly certain nothing is
overflowing at reasonable frequencies.  The 0-12 dB oscillation with 3dB
mean is something I should have expected, and is insignificant since it
pertains to a signal of manitude 0-30 while the majority are around
1x10^15 or much larger.

The initial issue had to do with the Goertzel output being noisy with a
mean of about 90dB and range of 20dB up and down.  As the signal
increased, the mean and variance first decreased, and then the mean began
increasing as expected.  Since I am performing detection based on
thresholding against ambient noise, this initial decrease limits the
sensitivity of the device.

By simulating a 60Hz input in software plus additive noise (approx.
white), I have been able to approximate the 90dB out with the 20dB range. 
However, as I raise the magnitude of the 60Hz signal, the mean rises
monotonically without the initial drop (once it overcomes the noise).  It
is looking less and less like a software issue and I may have to go back
to hunting on the hardware side of things.

>
>Note that 50 Hz and 60 Hz sinusoids will be orthogonal
>every 0.1 seconds of window length (every 400 samples
>in your example).  They will not be orthogonal over other
>durations of filtering, thus the Goertzel filter output
>will not be zero for those durations.  If you want only
>orthogonal correlations, then you can try running two
>Goertzel filters starting 400 samples apart and subtract.
>

I'm not certain how the orthogonality of the two is applied since the 50Hz
and 60Hz levels are being analyzed with two seperate implementations of the
goertzel and returning indepedant single values.  In this case are you
talking about using a goertzel alorithm's states to provide outputs of a
filter and then subtracting a two time (sample) delayed results?  I
haven't had any luck finding an explanation of this technique and would
appreciate if you could point me towards an example or literature.

Thank you,
Louis

superlou wrote:

   ...

> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
> 
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).  When I run the exact same
> test with an input 1/2 the size (16382 ampl.), the 60Hz energy retreats to
> 150dB (expected 6.1 dB decrease for 1/2 peak to peak range) but the 50Hz
> energy is constant at -108dB (no oscillation).  Decreasing the amplitude
> by another factor of 2 sees the 60Hz energy to 144dB and 50Hz energy to
> -114dB.  It looks almost like something is saturating in the algorithm
> with the full scale wave, but all the data types are double and I have not
> had any luck finding it yet.  Being fairly inexperience with DSP in general
> and the Goertzel algorithm in particular, I was wondering if anyone had run
> across a similar situation.

Have you verified that no intermediate calculations overflow or saturate 
at any frequency?

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

On Nov 15, 12:43 pm, "superlou" <lsim...@safeflight.com> wrote:
> >The gain for a Goertzel is the same as for a DFT (ignoring numerical
> >issues). So if you are feeding an input whose frequency matches that
> >of the Goertzel, then your output will by scaled up by N/2 (amplitude)
> >or (N/2)^2 for energy, where N is the number of samples. I don't think
> >you've told us how many samples you are feeding into the thing.
>
> >So expect your output to be 32767*N/2 for the expected amplitude.
>
> >Since you didn't give the number of samples, we can back calculate to
> >see how many your have
>
> >(if 156dB is amplitude)
> >156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
> >amplitude
>
> >(if 156dB is energy)
> >156 = 20 log(32767*N/2) -> N = 3851 samples
>
> >Are you using about 3850 samples?  If so then your numbers seem about
> >right.
>
> >Why don't you let us know: sample rate, number of samples, amplitude
> >of signal
>
> >Then we can check to see if the algo's output makes sense.
>
> >Clay
>
> The output of the Goertzel algorithm is in terms of amplitude squared, and
> so 10*log(output^2) = 20*log(output).  If I had received 157dB output
> it would have indicated 4321 data points and the output I gave was simply
> truncated down.  Altering the program to give me another two significant
> figure, I get 156.32dB which leads to 3996 points.  Since I'm using a 4000
> sample block of data for the Goertzel, this seems ok to me.
>
> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
>
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).

Note that 50 Hz and 60 Hz sinusoids will be orthogonal
every 0.1 seconds of window length (every 400 samples
in your example).  They will not be orthogonal over other
durations of filtering, thus the Goertzel filter output
will not be zero for those durations.  If you want only
orthogonal correlations, then you can try running two
Goertzel filters starting 400 samples apart and subtract.




IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M

On Nov 15, 3:43 pm, "superlou" <lsim...@safeflight.com> wrote:
> >The gain for a Goertzel is the same as for a DFT (ignoring numerical
> >issues). So if you are feeding an input whose frequency matches that
> >of the Goertzel, then your output will by scaled up by N/2 (amplitude)
> >or (N/2)^2 for energy, where N is the number of samples. I don't think
> >you've told us how many samples you are feeding into the thing.
>
> >So expect your output to be 32767*N/2 for the expected amplitude.
>
> >Since you didn't give the number of samples, we can back calculate to
> >see how many your have
>
> >(if 156dB is amplitude)
> >156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
> >amplitude
>
> >(if 156dB is energy)
> >156 = 20 log(32767*N/2) -> N = 3851 samples
>
> >Are you using about 3850 samples?  If so then your numbers seem about
> >right.
>
> >Why don't you let us know: sample rate, number of samples, amplitude
> >of signal
>
> >Then we can check to see if the algo's output makes sense.
>
> >Clay
>
> The output of the Goertzel algorithm is in terms of amplitude squared, and
> so 10*log(output^2) = 20*log(output).  If I had received 157dB output
> it would have indicated 4321 data points and the output I gave was simply
> truncated down.  Altering the program to give me another two significant
> figure, I get 156.32dB which leads to 3996 points.  Since I'm using a 4000
> sample block of data for the Goertzel, this seems ok to me.
>
> The signal amplitude was 32767 (range [-32767 32767]) and the sampling
> rate is 4000Hz.  In order to provide more than one energy reading per
> second, I am shifting in data from a 40 point input buffer and
> recalculating the energy after each shift in.  The input buffer size in my
> testing has ranged from about 40 to 400 in order to free up processing time
> for filtering and other processing, but the goertzel output remains the
> same (just gets calculated less often).
>
> The interesting thing I am noticing is that when I run the algorithm with
> the 32767 ampl. wave, the 60Hz energy level is a constant output while
> another Goertzel algorithm with a target frequency at 50Hz is oscillating
> between 0 and 12 dB (with an occasional -112).  When I run the exact same
> test with an input 1/2 the size (16382 ampl.), the 60Hz energy retreats to
> 150dB (expected 6.1 dB decrease for 1/2 peak to peak range) but the 50Hz
> energy is constant at -108dB (no oscillation).  Decreasing the amplitude
> by another factor of 2 sees the 60Hz energy to 144dB and 50Hz energy to
> -114dB.  It looks almost like something is saturating in the algorithm
> with the full scale wave, but all the data types are double and I have not
> had any luck finding it yet.  Being fairly inexperience with DSP in general
> and the Goertzel algorithm in particular, I was wondering if anyone had run
> across a similar situation.
>
> Thanks,
> Louis- Hide quoted text -
>
> - Show quoted text -

Hello Louis,

I think your code is running fine now. When you input 60Hz and you see
the expected value coming from the 60 Hz detector and a nearly zero
value coming from the 50Hz detector - this is what you should get. The
confusion you have stems from the dB representation of the numbers. If
your 50 Hz level is over 100dB below that of the 60 Hz level, you are
essentially getting zero of the 50Hz.  You know doubles have
essentially 16 decimal places in their mantissas. So when you subtract
two of them with similar scaling, a lot of the digits go away.

So I'm not suprised that you are seeing some non-linear stuff because
of numerical noise. But your detector is working with well over a
100dB range. How much do you need? Try normalizing your numbers to
32767*(4000/2) and displaying the numbers in a linear scale. I think
this will put stuff into perspective. If you recall the limit of
log(x) as x->0 is -infinity, so your near zero outputs on a log scale
will be all over the place, but in contrast to a full scale output,
they indicate essentially zero out.

IHTH,
Clay

>The gain for a Goertzel is the same as for a DFT (ignoring numerical
>issues). So if you are feeding an input whose frequency matches that
>of the Goertzel, then your output will by scaled up by N/2 (amplitude)
>or (N/2)^2 for energy, where N is the number of samples. I don't think
>you've told us how many samples you are feeding into the thing.
>
>So expect your output to be 32767*N/2 for the expected amplitude.
>
>Since you didn't give the number of samples, we can back calculate to
>see how many your have
>
>(if 156dB is amplitude)
>156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
>amplitude
>
>
>(if 156dB is energy)
>156 = 20 log(32767*N/2) -> N = 3851 samples
>
>Are you using about 3850 samples?  If so then your numbers seem about
>right.
>
>Why don't you let us know: sample rate, number of samples, amplitude
>of signal
>
>Then we can check to see if the algo's output makes sense.
>
>Clay

The output of the Goertzel algorithm is in terms of amplitude squared, and
so 10*log(output^2) = 20*log(output).  If I had received 157dB output
it would have indicated 4321 data points and the output I gave was simply
truncated down.  Altering the program to give me another two significant
figure, I get 156.32dB which leads to 3996 points.  Since I'm using a 4000
sample block of data for the Goertzel, this seems ok to me.

The signal amplitude was 32767 (range [-32767 32767]) and the sampling
rate is 4000Hz.  In order to provide more than one energy reading per
second, I am shifting in data from a 40 point input buffer and
recalculating the energy after each shift in.  The input buffer size in my
testing has ranged from about 40 to 400 in order to free up processing time
for filtering and other processing, but the goertzel output remains the
same (just gets calculated less often).

The interesting thing I am noticing is that when I run the algorithm with
the 32767 ampl. wave, the 60Hz energy level is a constant output while
another Goertzel algorithm with a target frequency at 50Hz is oscillating
between 0 and 12 dB (with an occasional -112).  When I run the exact same
test with an input 1/2 the size (16382 ampl.), the 60Hz energy retreats to
150dB (expected 6.1 dB decrease for 1/2 peak to peak range) but the 50Hz
energy is constant at -108dB (no oscillation).  Decreasing the amplitude
by another factor of 2 sees the 60Hz energy to 144dB and 50Hz energy to
-114dB.  It looks almost like something is saturating in the algorithm
with the full scale wave, but all the data types are double and I have not
had any luck finding it yet.  Being fairly inexperience with DSP in general
and the Goertzel algorithm in particular, I was wondering if anyone had run
across a similar situation.

Thanks,
Louis

On Nov 14, 3:59 pm, "superlou" <lsim...@safeflight.com> wrote:
> >Check to see what you get from your filter with _no_ input, or with some
> >very small input, or with some small initial value on the states.
> >Ideally you should be able to set the two states to some value (set them
> >both to the same one), run the filter for one cycle, and and see those
> >values back.
>
> After implementing one of the low-frequency algorithms to supplement the
> text-book goertzel I have been using, I have the following results and
> interpretations.  Please let me know if they make sense:
>
> 1) Given a constant input of 0 and states (q1,q2) reset to 0, both output
> 0 "dB".  Both q's always equal 0.  This makes intuitive sense following
> the math of the algorithms.
> 2) Inputting a full scale 60Hz sine wave (16bit, values -32768 to +32767,
> synthesized in software) leads to an output of 156dB and 6dB at 60 and 50
> Hz respectively.  If the states are instead set at each run to start at 1
> or 32767, the 50Hz value increases slightly to a max at 9dB while the 60Hz
> value stays exactly the same.  Both goertzel implementations behave
> virtually indentically for output, although the states for the textbook
> goertzel are similar positive or negative values while the alternative
> goertzel has one positive and one negative.
>
> In the above two cases, the algorithms were only run for one cycle (well,
> multiple cycles, but all information was reset each time).  Because my
> implementation is a sliding form and the real world will incur a random
> phase offset, I then tested using a synthasized 60Hz full scale wave and
> ran the algorithms over many cycles.  I found:
>
> 1) Given a fullscale wave, the output at 60Hz was a constant at 156dB
> while the 50Hz output oscillated 0 and 12 dB.
>
> 2) Given a 1/2 fullscale wave, the output at 60Hz is a constant 150.  This
> makes sense as I have reduced the wave amplitude by have and am observing a
> 6dB drop.  The 50Hz power oscillates around -173 to -161.
>
> I have examined both wave forms to verify there is no clipping.  In all
> situations the q values change largely (more so in case 1 than 2) but
> multiple combinations will result in the same power.  If my understanding
> is correct, this change is due to the phase information but I only care
> about the magnitude.  If it wasn't for the vast difference in the 50Hz and
> between (1) and (2), I would say it is behaving as I expect it.  Since the
> only change is input amplitude, I think the cause may be overflowing a
> variable somewhere and I am stepping through now trying to track it down.
>
> >Intuition tells me that with 4000 points, you need data paths that are
> >at least 12 bits wider than your ADC data path.  Mathematical analysis
> >would tell me more, but I'm not getting paid for this...
>
> I am gathering 16bit data and am using doubles (52-bit mantessa) as when I
> started the project, I experienced large quantization errors trying to
> perform a similar task using high order narrow bandpass IIR filters.
> Because the goertzel is an FFT in an IIR form to evaluate only the
> frequency of interest, I stuck with the double datatype.
>
> Since both forms of the goertzel algorithm perform almost identically, and
> the alternative algorithm is keyed to avoid quantization errors at low
> target frequencies, it also points against quantization errors.
>
> My analysis of the necessary data widths is very weak and I'm very
> interested in your rule of thumb saying I would need 12 bits wider than
> the ADC data path.  Based on your response, did you choose 12 bits as it
> is the lowest power of 2 to be greater than 4000?  If you know of any
> texts or have a moment to explain some of your intuition I would greatly
> appreciate it.
>
> Thank you,
> Louis

Hello Louis,

The gain for a Goertzel is the same as for a DFT (ignoring numerical
issues). So if you are feeding an input whose frequency matches that
of the Goertzel, then your output will by scaled up by N/2 (amplitude)
or (N/2)^2 for energy, where N is the number of samples. I don't think
you've told us how many samples you are feeding into the thing.

So expect your output to be 32767*N/2 for the expected amplitude.

Since you didn't give the number of samples, we can back calculate to
see how many your have

(if 156dB is amplitude)
156 = 10 log(32767*N/2) ->  N = 2.4E11 samples if 156dB is your
amplitude


(if 156dB is energy)
156 = 20 log(32767*N/2) -> N = 3851 samples

Are you using about 3850 samples?  If so then your numbers seem about
right.

Why don't you let us know: sample rate, number of samples, amplitude
of signal

Then we can check to see if the algo's output makes sense.

Clay

>superlou wrote:
>
>>   ... a full scale 60Hz sine wave (16bit, values -32768 to +32767
>
>has a DC component because of the asymmetry. -32767 to +32767 avoids
that.
>

I have modified the program to generate the wave over -32767 to 32767. 
The results seem the same, although I will examine this more and make sure
the real-world signals I have don't have a DC bias.  As far as predicting
it's impact, I thought that in the frequency domain, the DC offset would
create an impulse at 0 Hz and the narrow filter shouldn't react to it. 
When working with DSPs, will there be an effect?  I have actually noticed
a small DC offset of the signals I am sending to the ADC.  I am going to
try removing this DC offset by subtracting the mean from each set of data
before passing it to the Goertzel algorithm.

I apologize for posting that with zero input to the goertzel I receive 0
dB out.  Given that the output before converting to dB is 0, I should have
output -Inf and it is likely an issue with the Log10 function in that
situation.

superlou wrote:

>   ... a full scale 60Hz sine wave (16bit, values -32768 to +32767

has a DC component because of the asymmetry. -32767 to +32767 avoids that.

   ...

Jerry
-- 
Engineering is the art of making what you want from things you can get.
�����������������������������������������������������������������������

>Check to see what you get from your filter with _no_ input, or with some 
>very small input, or with some small initial value on the states. 
>Ideally you should be able to set the two states to some value (set them

>both to the same one), run the filter for one cycle, and and see those 
>values back.

After implementing one of the low-frequency algorithms to supplement the
text-book goertzel I have been using, I have the following results and
interpretations.  Please let me know if they make sense:

1) Given a constant input of 0 and states (q1,q2) reset to 0, both output
0 "dB".  Both q's always equal 0.  This makes intuitive sense following
the math of the algorithms.
2) Inputting a full scale 60Hz sine wave (16bit, values -32768 to +32767,
synthesized in software) leads to an output of 156dB and 6dB at 60 and 50
Hz respectively.  If the states are instead set at each run to start at 1
or 32767, the 50Hz value increases slightly to a max at 9dB while the 60Hz
value stays exactly the same.  Both goertzel implementations behave
virtually indentically for output, although the states for the textbook
goertzel are similar positive or negative values while the alternative
goertzel has one positive and one negative.

In the above two cases, the algorithms were only run for one cycle (well,
multiple cycles, but all information was reset each time).  Because my
implementation is a sliding form and the real world will incur a random
phase offset, I then tested using a synthasized 60Hz full scale wave and
ran the algorithms over many cycles.  I found:

1) Given a fullscale wave, the output at 60Hz was a constant at 156dB
while the 50Hz output oscillated 0 and 12 dB.

2) Given a 1/2 fullscale wave, the output at 60Hz is a constant 150.  This
makes sense as I have reduced the wave amplitude by have and am observing a
6dB drop.  The 50Hz power oscillates around -173 to -161.

I have examined both wave forms to verify there is no clipping.  In all
situations the q values change largely (more so in case 1 than 2) but
multiple combinations will result in the same power.  If my understanding
is correct, this change is due to the phase information but I only care
about the magnitude.  If it wasn't for the vast difference in the 50Hz and
between (1) and (2), I would say it is behaving as I expect it.  Since the
only change is input amplitude, I think the cause may be overflowing a
variable somewhere and I am stepping through now trying to track it down.

>Intuition tells me that with 4000 points, you need data paths that are 
>at least 12 bits wider than your ADC data path.  Mathematical analysis 
>would tell me more, but I'm not getting paid for this...

I am gathering 16bit data and am using doubles (52-bit mantessa) as when I
started the project, I experienced large quantization errors trying to
perform a similar task using high order narrow bandpass IIR filters. 
Because the goertzel is an FFT in an IIR form to evaluate only the
frequency of interest, I stuck with the double datatype.

Since both forms of the goertzel algorithm perform almost identically, and
the alternative algorithm is keyed to avoid quantization errors at low
target frequencies, it also points against quantization errors.

My analysis of the necessary data widths is very weak and I'm very
interested in your rule of thumb saying I would need 12 bits wider than
the ADC data path.  Based on your response, did you choose 12 bits as it
is the lowest power of 2 to be greater than 4000?  If you know of any
texts or have a moment to explain some of your intuition I would greatly
appreciate it.

Thank you,
Louis

superlou wrote:
> Clay and Tim,
> 
>> You are measuring variation in dB?  Relative to what?  This may be why
> you
>> are seeing an apparent change in the noise 'variance' -- the term
>> 'variance' is an absolute one, so your variance should be calculated in
>> ADC counts, not dB relative to the signal or some such.  The filter is a
>> linear one, so with the same amount of noise you should always see the
>> same absolute variance.
>>
> 
> I apologize for the vagueness of how I explained my problem and will start
> from the beginning.  I am new to DSP (a little theoretical background from
> texts) and welcome any corrections to my terminology.
> 
> I did not mean variance (stat. definition) but the range of the power
> signals in the time domain.  The dB levels were obtained simply by taking
> the magnitude squared output of a Goertzel algorithm and finding 10*log10
> of that value.  My naive approach was to find the maximum value when no
> intentional 50/60 Hz signal was applied and refer to this as a noise
> floor.  I would then measure any increase in the 50 or 60 Hz goertzel
> outputs compared to this noise floor to identify the magnitude of each
> 50/60 Hz signal.
> 
>>> With a
>>>
>>> small signal, the variance decreases and the mean drops to 80 dB. 
>> Hold it.  You just said with no signal the mean is 80dB.  It changes
> from
>> 80dB to 80dB?  80dB relative to what?
>>
> 
> The "dB" values are only valid for comparing power changes over time as
> they are really only relative to the ADC values.  With a signal generator,
> I verified that increasing an input 50 or 60 Hz by n dB resulted in an
> increase in the Goertzel's output by n dB.
> 
>>> As the
>>>
>>> desired signal in the noisy environment increases, the measured levels
>>>
>>> converge and match the changes in the input signal.
>> Certainly if you're expressing your output in dB this is what would
> appear
>> to be the case; as the signal swamps out the noise the dB variation due
> to
>> noise will, indeed, go down.
>>
>>> Since I am trying to determine the existance and strength of the
> desired
>>> 50/60 Hz tone by simply its rise above the noise floor, the
> convergance
>>> to
>>>
>>> underneath the apparent (no signal) noise floor is frustrating and
>>>
>>> puzzling. 
>> This is the first time that you've mentioned the filter variance being
>> less than the noise floor.  Explain things thoroughly, please.
>>
> 
> Given a test case starting out with no applied 50 or 60 Hz signal, the
> Goertzel power outputs vary between 60 and 100 dB with a mean near 80 dB. 
> Again, these dB are relative only to themselves.  Using my earlier
> assumption that with no applied signal, the values at 50 and 60 Hz will be
> a result of ambient noise, I said that 100 dB would be the top of this
> noise floor.  
> 
> What suprised me is that when a 60 Hz signal was applied, the 60 hz power
> first converged to 80 dB (<100dB) before increasing with increases in
> applied 60 Hz power.  From my limited noise analysis, I had thought that
> additive white noise would be a constant component of the power output.
> 
> While the system is definitly identifying a 60 Hz signal at that 80 dB
> mark, by simple thresholding, I can't say it exists until it passes 100dB
> and this significantly detracts from the sensitivity of the device.  I
> could identify it earlier by monitoring the system variance, but I was
> wondering if there is a way to eliminate this effect and avoid the
> overhead of the additional calculations.
> 
>>> I was wondering if anyone had any experience with this effect
>>>
>>> or had any suggestions to try.
>> How are you measuring noise?  If you have broadband noise, then you
> would
>> expect that a narrow filter such as you implement using the Goertzel
>> algorithm would, indeed, have output "below the noise floor".  It's why
>> you're using the thing, for goodness sake!
>>
>> You should expect that in the presence of white noise, the noise output
>> power of a narrow filter is the filter's bandwidth times the spectral
> power
>> density of the noise.
>>
>>> While I have had minimal luck with a
>>>
>>> moving average or windowed minema filter, the delay incurred is much
>>>
>>> larger than allowable.  I have tried various forms of windowing on the
>>>
>>> input data, but the system still displays the same characteristics. 
> The
>>> Geortzel algorithm is performed on a 4000 point data buffer with new
>>> data
>>>
>>> shifted in in 40 point segments.
>>>
>> It sounds like the algorithm is doing just as it should.
>>
> 
> While the implementation I have can currently detect 50 or 60 Hz signals,
> I would like to improve the sensitivity on the device.  I believe that
> there is room for improvement since an analog bandpass filtering device is
> able to detect a 60 Hz signal of a lower power level in the same
> environment.  My sampling rate is as low as I could get it via the
> available hardware, 4kHz.  I have not been using decimation in order to
> avoid the delay required for acquiring/filtering samples.  The goal
> requires a bandwidth of 1 Hz.  Clay, thank you for the link and I'm
> working on digesting those posts.

Check to see what you get from your filter with _no_ input, or with some 
very small input, or with some small initial value on the states. 
Ideally you should be able to set the two states to some value (set them 
both to the same one), run the filter for one cycle, and and see those 
values back.

Intuition tells me that with 4000 points, you need data paths that are 
at least 12 bits wider than your ADC data path.  Mathematical analysis 
would tell me more, but I'm not getting paid for this...

-- 

Tim Wescott
Wescott Design Services
http://www.wescottdesign.com

Do you need to implement control loops in software?
"Applied Control Theory for Embedded Systems" gives you just what it says.
See details at http://www.wescottdesign.com/actfes/actfes.html