Goertzel spectral response

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/688069085a5312fe

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) +
            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

On 18 Feb, 13:58, 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?

The 1/N term is placed according to preference. You need one
1/N scaling term somewhere in order to reproduce x[n] as IDFT{DFT{x
[n]}}
but you can assign it to the DFT, the IDFT or as two 1/sqrt(N) terms
as you please.

Just decide on one convension and stick with it.

...
> 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?

Again, this has to do with conventions. The most common formulas
for the DFT/IDFT pair goes as (view with fixed-width font)

       N-1
X[k] = sum x[n] exp(-j2pi n k /N)
       n=0

        1  N-1
x[n] = --- sum X[k] exp( j2pi n k /N)
        N  k=0

The reason would presumably be that one on average saves some
computations by leaving out the 1/N term in the forward transform,
since there are far more DFTs computed than IDFTs.

The price to be paid is that the norm is not preserved,

N-1              N=1
sum |x[n]|^2 =/= sum |X[k]|
n=0              k=0

and that the spectrum shows amplitudes that are off by a
scale factor related to N.

Again, this is not a big deal as long as one knows what is
going on and inserts the appropriate scaling factors when
needed.

Rune

Hi Rune,

thank you for the quick answer.

On 18 Feb., 14:18, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> The 1/N term is placed according to preference. You need one
> 1/N scaling term somewhere in order to reproduce x[n] as IDFT{DFT{x
> [n]}}
> but you can assign it to the DFT, the IDFT or as two 1/sqrt(N) terms
> as you please.
>
> Just decide on one convension and stick with it.

You are right, now I remember that I have learned that some time in
the past. So this question was a stupid one, sorry.

> > 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?
>
> Again, this has to do with conventions.
...
> The price to be paid is that the norm is not preserved,
>
> N-1              N=1
> sum |x[n]|^2 =/= sum |X[k]|
> n=0              k=0
>
> and that the spectrum shows amplitudes that are off by a
> scale factor related to N.

This was not quite what I was talking about. You are right of course
that without scaling factor, the norm is not preserved. But assuming
that the scaling factor is so that the spectral response is normalized
(1/N in my case), I nevertheless would expect it to differ from unity
in cases where we do not have an integer multiple of the period of the
center frequency inside the (rectangular) window.

Example:

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 an output
of 1. With a 1900 points window, the response is 1.049, and with 2100
points, it is 0.956. 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.

Matthias

On 18 Feb, 15:05, Matthias Zenker <matthiaszen...@nurfuerspam.de>
wrote:
>  But assuming
> that the scaling factor is so that the spectral response is normalized
> (1/N in my case), I nevertheless would expect it to differ from unity
> in cases where we do not have an integer multiple of the period of the
> center frequency inside the (rectangular) window.

Ah. You are talking about spectral leakage. I have never seen
any closed-form expressions to reconstruct the amplitude from
the spectrum in such cases. Even if if such expressions exist
(I don't know if they do) they would be less than trustworthy
if the signal contains additive noise.

Rune

On 18 Feb., 15:17, Rune Allnor <all...@tele.ntnu.no> wrote:

> Ah. You are talking about spectral leakage. I have never seen
> any closed-form expressions to reconstruct the amplitude from
> the spectrum in such cases. Even if if such expressions exist
> (I don't know if they do) they would be less than trustworthy
> if the signal contains additive noise.

Hmm.
So the spectral response function does not contain those effects?
I will have to stick with my numerical experiments, though.

Thank you,

Matthias

On 18 Feb, 15:34, Matthias Zenker <matthiaszen...@nurfuerspam.de>
wrote:
> On 18 Feb., 15:17, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > Ah. You are talking about spectral leakage. I have never seen
> > any closed-form expressions to reconstruct the amplitude from
> > the spectrum in such cases. Even if if such expressions exist
> > (I don't know if they do) they would be less than trustworthy
> > if the signal contains additive noise.
>
> Hmm.
> So the spectral response function does not contain those effects?

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 have to stick with my numerical experiments, though.

That's the best way to proceed, to build intuition.

Rune

On Feb 18, 8:34&#4294967295;am, Matthias Zenker <matthiaszen...@nurfuerspam.de>
wrote:
> On 18 Feb., 15:17, Rune Allnor <all...@tele.ntnu.no> wrote:
>
> > Ah. You are talking about spectral leakage. I have never seen
> > any closed-form expressions to reconstruct the amplitude from
> > the spectrum in such cases. Even if if such expressions exist
> > (I don't know if they do) they would be less than trustworthy
> > if the signal contains additive noise.
>
> Hmm.
> So the spectral response function does not contain those effects?
> I will have to stick with my numerical experiments, though.
>
> Thank you,
>
> Matthias

I usually normalize.  Notice 2000/1900 = 1.0526 and 2000/2100 = .
9524.  I put the normalization in the recurrsion for the algorithm.
However, I use the reciprocal 1900/2000 = .95, 2100/2000 = 1.05

Maurice Givens

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.


Matthias

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)
     = 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 =  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;   % Amplitude ratio for N = 1900 and N = 2000
r2 = X2/X0;   % Amplitude ratio for N = 2100 and N=2000

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