```On Nov 16, 5:31 am, Richard Owlett <rowl...@atlascomm.net> wrote:
> Could someone point me to source code for Goertzel in BASIC ... ?

Here's one which Clay Turner posted in another thread, converted
to generic BASIC:

k = 2 * sin(pi * f0 / sr)
y1 = 0
y2 = 0
for j = 1 to n
y2 = y2 - k * y1 + my_signal(j)
y1 = y1 + k * y2
next j
goertzel_result = 2 * sqrt(y1*y1 + y2*y2 - k*y1*y2) / n

Works in Chipmunk Basic: http://www.nicholson.com/rhn/basic

IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
```
```"Tim Wescott" schrieb
> > I'm not fluent in Scilab, what does the "%" in
> >       analysis_real = sin(2*%pi*time)
> > mean?
> > Shouldn't there be a frequency term somewhere
> > in this formula?
> > What does sin(12 seconds) mean?
> >
>
> I didn't see the "sin(12 seconds)", so I can't comment.
>
As Richard wrote, the frequency is implicitly assumed to be 1.
I was stumbling over the sin(time) term, so I was wondering
whether the frequency was missing.

> Scilab uses '%' to denote built-in variables;
> %pi is Pi: 3.14159 (etc). Similarly, %i is sqrt(-1),
>  and %e is the base of the natural logarithm.
>
Ahh, thanks.

> The frequency term in the above example is 2 Pi (2*%pi).
>
> The action of the Goertzel algorithm is to run the signal through
a
> bandpass filter with a resonant frequency at the frequency of
interest and
> damping ratio of zero for some period of time, then stop and look
at the
> amount of energy that's been collected by the filter.
>
> You'll also see the filter referred to as a "forced digital
oscillator" or
> as having "poles on the unit circle" -- all of these statements
are true,
> and mean the same thing. What they mean is that the filter has an
impulse
> response that's a perfect sine wave.
>
> The action of an IIR filter is to take it's input and convolve it
by it's
> impulse response.  Convolution is, well, convoluted, but in this
case it
> means multiplying the input by a sine wave, point by point, and
> the result.  That is what Richard's process is doing:  he's
multiplying
> the signal by a sine wave, and by the same thing in quadrature to
capture
> _all_ the energy, then at the end he's measuring the energy he's
captured
> in his filter.  The FIR filter part of his algorithm is a
one-frequency
> discrete Fourier transform, or DFT. The difference between his
filter and a
> 'true' Goertzel is that the Goertzel uses an IIR filter as it's
kernel,
> where Richard's DFT method uses an FIR filter.  Taken as a black
box, the
> two algorithms work the same -- they just require different
processor
> resources to achieve their goals with reasonable accuracy.
>
Thanks for the explanation.

Regards
Martin

```
```On Sat, 17 Nov 2007 17:02:54 +0100, Martin Blume wrote:

> "Tim Wescott"  schrieb
>>
>> > pseudo code in Scilab notation
>> >
>> > time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>> > signal = S[1:N];
>> > analysis_real=sin(2*%pi*time);
>> > analysis_imag=cos(2*%pi*time);
>> > output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>> >
>> > Is that Goertzel?
>> >
>> > or
>> > Could someone point me to source code for Goertzel
>> > in BASIC or FORTRAN?
>> > I've found C/C++ source but I don't read C.
>>
>> The Goertzel algorithm uses an IIR filter with special properties
> to
>> emulate an FIR filter.  You're implementing the FIR filter itself,
> which
>> is arguably a better approach to the task.
>>
>> Absent of odd numerical effects, your algorithm will do the
>> same as the Goertzel.
>>
> Sorry for being so dumb, could you elucidate a little bit more just
> how
> this is equivalent to a Goertzel algorithm? For instance, where does
> the frequency of interest for the filter come in?
>
> I'm not fluent in Scilab, what does the "%" in
>>  analysis_real = sin(2*%pi*time)
> mean? Shouldn't there be a frequency term somewhere in this formula?
> What does sin(12 seconds) mean?
>
> Martin

I didn't see the "sin(12 seconds)", so I can't comment.

Scilab uses '%' to denote built-in variables; %pi is Pi: 3.14159 (etc).
Similarly, %i is sqrt(-1), and %e is the base of the natural logarithm.

The frequency term in the above example is 2 Pi (2*%pi).

The action of the Goertzel algorithm is to run the signal through a
bandpass filter with a resonant frequency at the frequency of interest and
damping ratio of zero for some period of time, then stop and look at the
amount of energy that's been collected by the filter.

You'll also see the filter referred to as a "forced digital oscillator" or
as having "poles on the unit circle" -- all of these statements are true,
and mean the same thing. What they mean is that the filter has an impulse
response that's a perfect sine wave.

The action of an IIR filter is to take it's input and convolve it by it's
impulse response.  Convolution is, well, convoluted, but in this case it
means multiplying the input by a sine wave, point by point, and adding up
the result.  That is what Richard's process is doing:  he's multiplying
the signal by a sine wave, and by the same thing in quadrature to capture
_all_ the energy, then at the end he's measuring the energy he's captured
in his filter.  The FIR filter part of his algorithm is a one-frequency
discrete Fourier transform, or DFT. The difference between his filter and a
'true' Goertzel is that the Goertzel uses an IIR filter as it's kernel,
where Richard's DFT method uses an FIR filter.  Taken as a black box, the
two algorithms work the same -- they just require different processor
resources to achieve their goals with reasonable accuracy.

--
Tim Wescott
Control systems and communications consulting
http://www.wescottdesign.com

Need to learn how to apply control theory in your embedded system?
"Applied Control Theory for Embedded Systems" by Tim Wescott
Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html
```
```Martin Blume wrote:

> "Tim Wescott"  schrieb
>
>>>pseudo code in Scilab notation
>>>
>>>time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>>>signal = S[1:N];
>>>analysis_real=sin(2*%pi*time);
>>>analysis_imag=cos(2*%pi*time);
>>>output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>>>
>>>Is that Goertzel?
>>>
>>>or
>>>Could someone point me to source code for Goertzel
>>>in BASIC or FORTRAN?
>>>I've found C/C++ source but I don't read C.
>>
>>The Goertzel algorithm uses an IIR filter with special properties
>
> to
>
>>emulate an FIR filter.  You're implementing the FIR filter itself,
>
> which
>
>>is arguably a better approach to the task.
>>
>>Absent of odd numerical effects, your algorithm will do the
>>same as the Goertzel.
>>
>
> Sorry for being so dumb, could you elucidate a little bit more just
> how
> this is equivalent to a Goertzel algorithm? For instance, where does
> the frequency of interest for the filter come in?
>
> I'm not fluent in Scilab, what does the "%" in
>
>> analysis_real = sin(2*%pi*time)
>
> mean? Shouldn't there be a frequency term somewhere in this formula?
> What does sin(12 seconds) mean?
>
> Martin
>

first.

Scilab uses % as a prefix for predefined constants. In this case %pi is
3.14....
Frequency has an implied value of 1. It would have been clearer if I had
written
analysis_real=sin(2*%pi*analysis_frequency*time)
[ Tim automatically translates from Owlett to reality ;o>

```
```On Sat, 17 Nov 2007 07:02:10 -0600, Richard Owlett wrote:

> Tim Wescott wrote:
>> Richard Owlett wrote:
>>
>>> Tim Wescott wrote:
>>>
>>>> On Fri, 16 Nov 2007 07:31:39 -0600, Richard Owlett wrote:
>>>>
>>>>
>>>>> pseudo code in Scilab notation
>>>>>
>>>>> time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>>>>> signal = S[1:N];
>>>>> analysis_real=sin(2*%pi*time);
>>>>> analysis_imag=cos(2*%pi*time);
>>>>> output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>>>>>
>>>>> Is that Goertzel?
>>>>>
>>>>> or
>>>>> Could someone point me to source code for Goertzel in BASIC or FORTRAN?
>>>>> I've found C/C++ source but I don't read C.
>>>>
>>>>
>>>>
>>>> The Goertzel algorithm uses an IIR filter with special properties to
>>>> emulate an FIR filter.  You're implementing the FIR filter itself, which
>>>> is arguably a better approach to the task.
>>>>
>>>> Absent of odd numerical effects, your algorithm will do the same as the
>>>> Goertzel.
>>>>
>>>
>>> I take it then that its "bandwidth" would the same dependence on
>>> sample rate and number of samples used to do the calculation.
>>>
>> Yes.  Instead of saying "Absent of ..." I should have said "Your
>> algorithm does _exactly_ what the Goertzel does".  Conclusions for A
>> apply to B, and visa versa.
>>
>
> Hmmm. Is it valid to stretch things a little further.
>
> In my pseudo code, I wrote "signal = S[1:N]" with the explicit
> presumption that S[1:N] was a pure sinusoid. What if I replaced it with
> the DFT of a desired frequency response? If viable, how would I specify
> the frequency response to account for the folding around Nyquist
> frequency that is observed when going from time --> frequency domain?

Do you mean if you used the DFT instead of a Goertzel, or if you replaced

In the former case, the DFT and Goertzel will still be identical.  They'll
act on S[1:N], _as sampled_, which means the aliasing has already
happened.  If you want to show the response going from frequencies in the
time domain to the output of your DFT or Goertzel, then figure out the
effects of any anti-alias filtering, then the effects of sampling, then
the effect of the sampled-time filtering.

For more about sampling see this:
http://www.wescottdesign.com/articles/Sampling/sampling.html.  For even
more, see a good book on DSP (such as Rick Lyons' book).

--
Tim Wescott
Control systems and communications consulting
http://www.wescottdesign.com

Need to learn how to apply control theory in your embedded system?
"Applied Control Theory for Embedded Systems" by Tim Wescott
Elsevier/Newnes, http://www.wescottdesign.com/actfes/actfes.html
```
```"Tim Wescott"  schrieb
>
> > pseudo code in Scilab notation
> >
> > time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
> > signal = S[1:N];
> > analysis_real=sin(2*%pi*time);
> > analysis_imag=cos(2*%pi*time);
> > output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
> >
> > Is that Goertzel?
> >
> > or
> > Could someone point me to source code for Goertzel
> > in BASIC or FORTRAN?
> > I've found C/C++ source but I don't read C.
>
> The Goertzel algorithm uses an IIR filter with special properties
to
> emulate an FIR filter.  You're implementing the FIR filter itself,
which
> is arguably a better approach to the task.
>
> Absent of odd numerical effects, your algorithm will do the
> same as the Goertzel.
>
Sorry for being so dumb, could you elucidate a little bit more just
how
this is equivalent to a Goertzel algorithm? For instance, where does
the frequency of interest for the filter come in?

I'm not fluent in Scilab, what does the "%" in
>  analysis_real = sin(2*%pi*time)
mean? Shouldn't there be a frequency term somewhere in this formula?
What does sin(12 seconds) mean?

Martin

```
```Tim Wescott wrote:
> Richard Owlett wrote:
>
>> Tim Wescott wrote:
>>
>>> On Fri, 16 Nov 2007 07:31:39 -0600, Richard Owlett wrote:
>>>
>>>
>>>> pseudo code in Scilab notation
>>>>
>>>> time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>>>> signal = S[1:N];
>>>> analysis_real=sin(2*%pi*time);
>>>> analysis_imag=cos(2*%pi*time);
>>>> output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>>>>
>>>> Is that Goertzel?
>>>>
>>>> or
>>>> Could someone point me to source code for Goertzel in BASIC or FORTRAN?
>>>> I've found C/C++ source but I don't read C.
>>>
>>>
>>>
>>> The Goertzel algorithm uses an IIR filter with special properties to
>>> emulate an FIR filter.  You're implementing the FIR filter itself, which
>>> is arguably a better approach to the task.
>>>
>>> Absent of odd numerical effects, your algorithm will do the same as the
>>> Goertzel.
>>>
>>
>> I take it then that its "bandwidth" would the same dependence on
>> sample rate and number of samples used to do the calculation.
>>
> Yes.  Instead of saying "Absent of ..." I should have said "Your
> algorithm does _exactly_ what the Goertzel does".  Conclusions for A
> apply to B, and visa versa.
>

Hmmm. Is it valid to stretch things a little further.

In my pseudo code, I wrote "signal = S[1:N]" with the explicit
presumption that S[1:N] was a pure sinusoid. What if I replaced it with
the DFT of a desired frequency response? If viable, how would I specify
the frequency response to account for the folding around Nyquist
frequency that is observed when going from time --> frequency domain?

```
```Richard Owlett wrote:
> Tim Wescott wrote:
>
>> On Fri, 16 Nov 2007 07:31:39 -0600, Richard Owlett wrote:
>>
>>
>>> pseudo code in Scilab notation
>>>
>>> time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>>> signal = S[1:N];
>>> analysis_real=sin(2*%pi*time);
>>> analysis_imag=cos(2*%pi*time);
>>> output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>>>
>>> Is that Goertzel?
>>>
>>> or
>>> Could someone point me to source code for Goertzel in BASIC or FORTRAN?
>>> I've found C/C++ source but I don't read C.
>>
>>
>> The Goertzel algorithm uses an IIR filter with special properties to
>> emulate an FIR filter.  You're implementing the FIR filter itself, which
>> is arguably a better approach to the task.
>>
>> Absent of odd numerical effects, your algorithm will do the same as the
>> Goertzel.
>>
>
> I take it then that its "bandwidth" would the same dependence on sample
> rate and number of samples used to do the calculation.
>
Yes.  Instead of saying "Absent of ..." I should have said "Your
algorithm does _exactly_ what the Goertzel does".  Conclusions for A
apply to B, and visa versa.

--

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
```
```On Nov 16, 5:31 am, Richard Owlett <rowl...@atlascomm.net> wrote:
> pseudo code in Scilab notation
>
> time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
> signal = S[1:N];
> analysis_real=sin(2*%pi*time);
> analysis_imag=cos(2*%pi*time);
> output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>
> Is that Goertzel?

Goertzel seems to refer to the filter method above, but
with calculating the sine and cosine coefficients using
a recurrence formula in the form of an IIR resonator (but
I've seen using various other recurrences referred to here
as Goertzel-XYZ for various XYZ).

They work roughly the same as a 1-bin DFT, but have
slightly different accuracy and efficiency characteristics.

IMHO. YMMV.
--
rhn A.T nicholson d.0.t C-o-M
```
```Tim Wescott wrote:

> On Fri, 16 Nov 2007 07:31:39 -0600, Richard Owlett wrote:
>
>
>>pseudo code in Scilab notation
>>
>>time=[start_time:(1/sample_rate):(start_time+N/sample_rate)];
>>signal = S[1:N];
>>analysis_real=sin(2*%pi*time);
>>analysis_imag=cos(2*%pi*time);
>>output=((signal*analysis_real')^2 +(signal*analysis_imag')^2)^.5
>>
>>Is that Goertzel?
>>
>>or
>>Could someone point me to source code for Goertzel in BASIC or FORTRAN?
>>I've found C/C++ source but I don't read C.
>
>
> The Goertzel algorithm uses an IIR filter with special properties to
> emulate an FIR filter.  You're implementing the FIR filter itself, which
> is arguably a better approach to the task.
>
> Absent of odd numerical effects, your algorithm will do the same as the
> Goertzel.
>

I take it then that its "bandwidth" would the same dependence on sample
rate and number of samples used to do the calculation.

```