Forums

Goertzel spectral response

Started by Unknown February 18, 2009
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&#2013266080;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&#2013266080;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) > &#2013266080; &#2013266080; &#2013266080;= 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 = &#2013266080;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; &#2013266080; % Amplitude ratio for N = 1900 and N = 2000 > r2 = X2/X0; &#2013266080; % 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