The complex Goertzel is identical to computing one bin of
a DFT, or, more efficiently, an FFT.

Because of the inherent finite rectangular window, a DFT
will sample a Sinc function for a single frequency sinusoidal
time domain input. For a bin-centered frequency (a sinusoid
perfectly periodic in the filer aperture) that's one bin
(plus the negative bin if you count those).  For a non-bin
centered frequency, the Sinc splatters energy all over
the spectrum.

The DFT result is not quite a Sinc, since the sampling will
create an infinite train of Sinc functions, or in summation,
a Dirichlet kernel, which does look pretty close to a Sinc
for large N.

And, for a real signal, there is also a negative frequency
image, which differs in effect depending on phase of the
sinusoid relative to the center of the aperture, and has
greater influence if your sinusoid is near DC or Fs/2 in
frequency.

I did my best at a closed form description of the frequency
response of a complex Goertzel filter on this web page:
  http://www.nicholson.com/rhn/dsp.html#4


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

On 23 Feb, 09:23, Matthias Zenker <matthiaszen...@nurfuerspam.de>
wrote:

> Instead of doing more or less random numeric experiments, I would
> rather like to have an analytic formula showing those dependencies on
> number of samples and center frequency. Does anyone know of such a
> formula or a way to find it, or where to look for it?

You need to look up the sampling theorem, and maybe also
the signal reconstruction theorem. These two theorems
together describe how a set of discrete samples can
represent an analog sampled signal; and how an analog
signal can be reconstructed from a set of samples.

You should find these things in most texts on DSP.

Rune

Thank you for this explanation, Clay.

On Feb 20, 4:35&#4294967295;pm, c...@claysturner.com wrote:

> BUT just from an analysis viewpoint the frequency
> variable in the DFT doesn't have to be an integer. But when it is not
> an integer you end up with an integral number of cycles and a partial
> cycle.

This is my case. The sample rate is fixed, the frequency is changing.

> The integral number of cycles will integrate (sum up) to a nice
> value. But the partial cycle will contribute a phase dependent
> fraction to the total result.

Yes.

> For example if you feed a single frequency into your DFT that is tuned
> to the exact same frequency and this frequency is not an integral
> frequency and you look at a series of outputs (formed from sucessive
> DFTs) you will notice the magnitude fluctuates due to the phase
> dependent partial cycle.

That is what I see.

> If your data span includes many cycles, then
> you can dilute this phase dependent portion.

It doesn't.

> But if you stick with
> integral frequencies and accept the leakage, at least you will see
> sucessive magnitudes be the same.

I can't, because I am monitoring a signal with moving ground frequency
with a fixed sample rate.

> The way you go with this depends on
> how accurately you need your estimate to be.

Well, I have seen that when I center my Goertzel at a significantly
deviating frequency (which is NOT in the true spectrum), I still see a
significant amplitude.

Example:

f_true = 350 Hz
Amplitude = 1

f_sample = 25 kHz

N_samples for Goertzel: 72 (corresponding to 347 Hz)
response of Goertzel @ 350 Hz: 0.99
response of Goertzel @ 200 Hz: 0.91
response of Goertzel @ 450 Hz: 0.76

N_samples for Goertzel: 56 (corresponding to 446 Hz)
response of Goertzel @ 350 Hz: 1.04
response of Goertzel @ 450 Hz: 0.94

This is definitively too unspecific.

Instead of doing more or less random numeric experiments, I would
rather like to have an analytic formula showing those dependencies on
number of samples and center frequency. Does anyone know of such a
formula or a way to find it, or where to look for it?

Thank you,

Matthias

On Feb 20, 6:17&#4294967295;am, Matthias Zenker <matthiaszen...@nurfuerspam.de>
wrote:
> First of all, thanks to all of you for your answers.
>
> Maybe I should attempt to explain more clearly what I am trying to
> find out.
>
> I have a signal consisting of a ground frequency and some of its
> harmonics. I want to determine the amplitude of one of the harmonics.
> We can assume in the first place that really only those frequencies
> are present in the signal. The signal may be continuous, but may also
> contain only one period of the ground frequency. The sampling
> frequency is fixed, but the ground frequency is not, so one period of
> the ground frequency may or may not contain an integer number of
> samples. The number of samples in one period is at least 50.
>
> At the moment, I determine the number of samples coming closest to the
> ground frequency, and then use a rectangular window and the goertzel
> algorithm to find the amplitude of the harmonic. This seems to work OK
> so far. Now I want to examine more closely inhowfar this method can be
> trusted, in particular:
>
> * What is the error induced by the non-integer period? It is clear
> that if the window length does not correspond to the period length, I
> have additional frequencies introduced by the small snippet of a wave
> at the end (or by the missing snippet of wave at the end). So it must
> be more than a normalization problem.
> * How sensitive is my calculation to a disturbing signal with a
> frequency near the disired harmonic? My impression is that this might
> be a major problem of this method.
> * I have observed that the calculated amplitude at the harmonic
> frequency is influenced by the amplitude of the ground frequency. (I
> have tested that with a signal like y1(t)=10*sin(2*pi*f0*t)+sin
> (2*pi*3*f0*t), compared to y2(t)=sin(2*pi*3*f0*t))Does this mean that
> there are nonlinearities here?
>
> So my question is if there is way to address this kind of problems
> analytically. Perhaps Rune is right, and all this is about spectral
> leakage, and not so much the Goertzel algorithm. So I maybe would need
> to learn more about spectral leakage, and how to overcome it.
> Any ideas, hints, pointers would be welcome.
>
> Thank you,
>
> Matthias

Well in a nutshell this is all about spectral leakage. The Goertzel
algo is simply an efficient way to calculate a DFT for a single
frequency. Normally the bin frequency for the DFT is chosen to have a
whole number of cycles within the span of data covered by the DFT. And
if you wish to involve Parseval's theorem then you will pretty much
need to stay with the bin frequencies that fit the whole number of
cycles thing. BUT just from an analysis viewpoint the frequency
variable in the DFT doesn't have to be an integer. But when it is not
an integer you end up with an integral number of cycles and a partial
cycle. The integral number of cycles will integrate (sum up) to a nice
value. But the partial cycle will contribute a phase dependent
fraction to the total result.

For example if you feed a single frequency into your DFT that is tuned
to the exact same frequency and this frequency is not an integral
frequency and you look at a series of outputs (formed from sucessive
DFTs) you will notice the magnitude fluctuates due to the phase
dependent partial cycle. If your data span includes many cycles, then
you can dilute this phase dependent portion. But if you stick with
integral frequencies and accept the leakage, at least you will see
sucessive magnitudes be the same. The way you go with this depends on
how accurately you need your estimate to be. If you know something
about your expected frequencies, then you can chose your sampling rate
and sample size to place your expected frequencies on or near the
integral bin frequencies.

IHTH,
Clay

First of all, thanks to all of you for your answers.

Maybe I should attempt to explain more clearly what I am trying to
find out.

I have a signal consisting of a ground frequency and some of its
harmonics. I want to determine the amplitude of one of the harmonics.
We can assume in the first place that really only those frequencies
are present in the signal. The signal may be continuous, but may also
contain only one period of the ground frequency. The sampling
frequency is fixed, but the ground frequency is not, so one period of
the ground frequency may or may not contain an integer number of
samples. The number of samples in one period is at least 50.

At the moment, I determine the number of samples coming closest to the
ground frequency, and then use a rectangular window and the goertzel
algorithm to find the amplitude of the harmonic. This seems to work OK
so far. Now I want to examine more closely inhowfar this method can be
trusted, in particular:

* What is the error induced by the non-integer period? It is clear
that if the window length does not correspond to the period length, I
have additional frequencies introduced by the small snippet of a wave
at the end (or by the missing snippet of wave at the end). So it must
be more than a normalization problem.
* How sensitive is my calculation to a disturbing signal with a
frequency near the disired harmonic? My impression is that this might
be a major problem of this method.
* I have observed that the calculated amplitude at the harmonic
frequency is influenced by the amplitude of the ground frequency. (I
have tested that with a signal like y1(t)=10*sin(2*pi*f0*t)+sin
(2*pi*3*f0*t), compared to y2(t)=sin(2*pi*3*f0*t))Does this mean that
there are nonlinearities here?

So my question is if there is way to address this kind of problems
analytically. Perhaps Rune is right, and all this is about spectral
leakage, and not so much the Goertzel algorithm. So I maybe would need
to learn more about spectral leakage, and how to overcome it.
Any ideas, hints, pointers would be welcome.

Thank you,

Matthias

On Feb 18, 7:58&#4294967295;am, matthiaszen...@nurfuerspam.de wrote:
> Hi,
>
> as a partial newbie in digital signal processing (I have read some
> basics in textbooks, namely the one by Proakis and Manolakis, but
> never done heavy applications in FFTs, filtering and such), I am
> trying to learn more about the Goetzel algorithm, namely it's spectral
> response, and found this thread (old, but not outdated ;) ):http://groups.google.com/group/comp.dsp/browse_thread/thread/68806908...
>
> My questions are:
>
> > "Clay S. Turner" <Phys...@Bellsouth.net> wrote in messagenews:4_Mte.149703$J25.60996@bignews6.bellsouth.net...
>
> > > The gain formula is
>
> > > X(f) = (1/N)*sin(pi*f*N/F)/sin(pi*f/F)
>
> 1.
> I have succeeded to reproduce Clay's formula above except for the
> leading 1/N. I did a (manual) DFT of
>
> W_N^-kn = exp(i*2*pi*k*n/N)
>
> which gives me
>
> X(f) = exp(i*pi*(N-1)*f/F) * sin(pi*f*N/F) / sin(pi*f/F)
> |X(f)| = sin(pi*f*N/F)/sin(pi*f/F)
>
> with f = f0-fx (f0: center frequency, fx: current frequency)
>
> but I may have made a math error since the 1/N is not there. Or does
> it come from an extra source I have missed in my starting point?
>
> 2.
> There was another solution by Tim Wescott:
>
> >Barring stupid math errors on my part such a filter will have a
> >frequency-domain shape of
>
> >H(theta) = sin((q0 - theta) * N/2) / (q0 - theta) +
>
> &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; &#4294967295; sin((q0 + theta) * N/2) / (q0 + theta)
>
>
>
> Why do Clay's and Tim's formulae differ?
>
> 2. The formula above (with leading 1/N) says that for f0=fx, the
> response is always unity. I would expect that if the window length is
> not an integer multiple of 1/fx, the gain would not be unity. I have
> run "manual" tests (i.e. apply a Goertzel to a test signal, varying
> the window length) which support this expectation. Why doesn't the
> gain formula reflect this result?
>
> Thank you for any hint,
>
> Matthias

Hello Matthias,

I added the (1/N) factor upfront to the Goertzel, since the Goertzel
is often used to measure the amount of a frequency in a signal.

So from the point of view of analysis, I used a 1/N factor to
normalize the output to be unity if you feed the algo a pure sine
exactly centered on a bin frequency.

Rune is correct that in an analysis - synthesis application the round
trip gain is N, so you are free to split the correction factors
between the analysis and synthesis stages in any way you choose.

Most often I see the forward DFT use a correction factor of 1, and
then a correction of 1/N is used on the back end. This is what
Goertzel used in his original paper.  The user needs to be aware of
the "gain" of the common (standard way it is implemented) gain of
Goertzel's method especially in a fixed point app.

A difference between mine and Tim's formulae is I gave you a periodic
sync function and Tim just gave you a pair of sinc functions. I
believe mine is the correct one although for small delta frequencies,
the difference is negligible.

Clay

On Feb 19, 10:41&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 19 Feb, 17:07, maury <maury...@core.com> wrote:
>
> > On Feb 19, 6:50&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:
> > > This question has nothing to do with the Goertzel algorithm,
> > > but with spectral leakage.
> > Rune,
> > I'm going to agree and disagree at the same time :). &#4294967295;I think there
> > are two things happening.
>
> As I understand it, Goertzel is just an algorithm
> to compute DFT coefficients, just like the FFT [*].
> So discussing the output of the Goertzel with an N
> point data series is no different than discussing
> the DFT with an N point data series.
>
> Except for run-time or numerical issues, but that's
> not an issue in this context.
>
> >&#4294967295;Matthias is correct that the output of the
> > Goertzel will vary as the number of samples analyzed is varied.
>
> Agreed.
>
> Again, my point was that this amplitude scaling effect
> will happen even if one uses some other algorithm than
> Goertzel, and that the basic effect has nothing to do
> with the Goertzel as such. The OP stated quite clearly
> that he used specific number of samples, in which case
> he actually computes N-pt DFT coefficients.
>
> It's far easier to understand the causes and effects
> when the obfuscating factors have been removed.
>
> Rune
>
> [*] Yes, I know: The FFT computes *all* the DFT coefficients
> &#4294967295; &#4294967295; while Goertzel computes only one. Both of them compute
> &#4294967295; &#4294967295; DFT coefficients, though.



I agree.  I got the impression the OP was looking for a way to
mitigate the output level problem as a function of the number of
samples.

Maurice

On 19 Feb, 17:07, maury <maury...@core.com> wrote:
> On Feb 19, 6:50&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:

> > This question has nothing to do with the Goertzel algorithm,
> > but with spectral leakage.

> Rune,
> I'm going to agree and disagree at the same time :). &#4294967295;I think there
> are two things happening.

As I understand it, Goertzel is just an algorithm
to compute DFT coefficients, just like the FFT [*].
So discussing the output of the Goertzel with an N
point data series is no different than discussing
the DFT with an N point data series.

Except for run-time or numerical issues, but that's
not an issue in this context.

>&#4294967295;Matthias is correct that the output of the
> Goertzel will vary as the number of samples analyzed is varied.

Agreed.

Again, my point was that this amplitude scaling effect
will happen even if one uses some other algorithm than
Goertzel, and that the basic effect has nothing to do
with the Goertzel as such. The OP stated quite clearly
that he used specific number of samples, in which case
he actually computes N-pt DFT coefficients.

It's far easier to understand the causes and effects
when the obfuscating factors have been removed.

Rune

[*] Yes, I know: The FFT computes *all* the DFT coefficients
    while Goertzel computes only one. Both of them compute
    DFT coefficients, though.

On Feb 19, 6:50&#4294967295;am, Rune Allnor <all...@tele.ntnu.no> wrote:
> On 19 Feb, 09:20, Matthias Zenker <matthiaszen...@nurfuerspam.de>
> wrote:
>
>
>
>
>
> > On 18 Feb., 15:55, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > > The spectral response function is a difficult animal to
> > > discuss at the best of times, since it kind of plays up
> > > to human intuition when f = k/N, but not for f =/= k/N.
> > > I must admit I have lost track of what your question really
> > > was, once the 1/N scaling terms were taken care of.
>
> > I will repeat the example I gave in a previous post, and try to
> > elaborate to make it clearer:
>
> > A 500 Hz sinusoidal signal is sampled with 1 MHz. One period
> > corresponds to 2000 points. A Goertzel filter with 500 Hz center
> > frequency and a 2000 points rectangular window would give me a
> > normalized output
> > of 1, i.e. it will give me the amplitude of the sine (the value a in y
> > = a*sin(2*pi*f*t)).
> > With a 1900 points window, the response is 1.049*a (i.e. 5% more than
> > the original amplitude), and with 2100
> > points, it is 0.956*a (4.4% less than originally). That is what I try
> > to find a formula for. The
> > naively calculated spectral response formula does not give this, it
> > always gives unity at the center frequency.
>
> This question has nothing to do with the Goertzel algorithm,
> but with spectral leakage. Below is a matlab script where I
> reproduce the cases you describe as plain DFT coefficients.
>
> As you see, the signals are all on the same form,
>
> x[n] = cos(2*pi*F/Fs*n)
> &#4294967295; &#4294967295; &#4294967295;= 1/2(exp(j2*pi*F/Fs*n)+exp(-j2*pi*F/Fs*n))
>
> (Euler's formula). So you expect the inner product
>
> y = <x[n], exp(-j2*pi*F/Fs*n)>/N
>
> to give
>
> y = 1/2.
>
> This actually happens in case 0, when N0 = 2000, since
> in that case the normalized frequency f becomes
>
> f = F/Fs = 500/1e6 = &#4294967295;1/2000 = 1/N0
>
> which is a rational number.
>
> In case 0, with N=2000 samples, this rational number
> *coincidentially* happens to be expressed in terms of
> the sequence length N = 2000. In which case the physical
> frequency F happens to coincide with a spectrum coefficient.
>
> In the other two case, however, the normalized frequency
> f still equals 1/2000 but now N = N1 =1900.
>
> In this case there is no integer ratio bewteen f and N1,
>
> f = k/N1
> 1/2000 = k/1900
> k = 1900/2000 = 19/20
>
> k is not an integer, so in this case the spectrum leakage
> effect does its work. The same happens for N2 = 2100.
>
> Rune
>
> %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
> Fs = 1e6;
> N0 = 2000;
> N1 = 1900;
> N2 = 2100;
> F = 500;
>
> x0= cos((0:N0-1)'*2*pi*F/Fs);
> x1= cos((0:N1-1)'*2*pi*F/Fs);
> x2= cos((0:N2-1)'*2*pi*F/Fs);
>
> X0= exp(i*(0:N0-1)'*2*pi*F/Fs)'*x0/N0;
> X1= exp(i*(0:N1-1)'*2*pi*F/Fs)'*x1/N1;
> X2= exp(i*(0:N2-1)'*2*pi*F/Fs)'*x2/N2;
>
> r1 = X1/X0; &#4294967295; % Amplitude ratio for N = 1900 and N = 2000
> r2 = X2/X0; &#4294967295; % Amplitude ratio for N = 2100 and N=2000- Hide quoted text -
>
> - Show quoted text -

Rune,
I'm going to agree and disagree at the same time :).  I think there
are two things happening.  Matthias is correct that the output of the
Goertzel will vary as the number of samples analyzed is varied.  If he
has 2000 sample points, and the input is 500 Hz, the output level will
vary as a function of the number of sample points.  Below is what I
use to normalize that problem and have the algorithm give me the
actual input level IF THE INPUT SIGNAL IS AT THE DETECTION.  If the
frequency is off, say 510 Hz instead of 500 Hz, then the spectral
leakage is apparent.

Maurice Givens


function y = goert(x,f,N,Fs)

% function y = goert(x,f,N,Fs)
%
% goertzel frequency detection algorithm
%
% x is input signal
% f is frequency to detect
% N is number of samples
% Fs is the sampling frequency
%
for m = 1:N:length(x)
k = fix(0.5 + N*f/Fs);
w = (2*pi/N)*k;
cosine = cos(w);
sine = sin(w);
coeff = 2 * cosine;
Q1 = 0;
Q2 = 0;
for j = 1:N
    Q0 = coeff * Q1 - Q2 + x(j + m - 1);
    Q2 = Q1;
    Q1 = Q0;
end
real_y = Q1 - Q2 * cosine;
imag_y = Q2 * sine;
mag_y = sqrt(real_y^2 + imag_y^2);
y(m:m+N-1) = sqrt(Q1^2 + Q2^2 - (Q1 * Q2 * coeff))*2/N;
end