Help with understanding the magnitude of a Goertzel function output

Hello all,

I'm writing some code to break down the data recorded from a microphone into specific frequencies. Initially I was trying to use fft(), but that made it too difficult to select the frequencies. Using the goertzel() function I have the data expressed in the frequency bins correctly, but I don't understand the magnitude that I am getting.

The raw data is in pa (Pascal) as a sound pressure level. It's fine to stay in that format for now. Below is a sample of the raw data, these magnitudes are typical.

-0.8256
0.3371
1.6520
1.6950
0.9749
1.0670
0.0441
-0.5052
0.2395
0.4317
0.3432
-0.0444
-0.1115
-0.3343

The full dataset has 65536 points collected over 5.12 seconds, so I wont include it here. My code is below. It outputs a stem plot

SPL = importspl('SPLsample.xlsx','Sheet1','B1:B65536');

N = length(SPL)
Fs = 5000  %5kHz sample freq

f = [12.5 16 20 25 31.5 40 50 63 80 100 125 160 200 250 315 400 500 630 800 1000 1250 1600 2000 2500 3150 4000];
freq_indices = round(f*N/Fs) + 1;
dft_data = goertzel(SPL, freq_indices);

% plot dft output on stem plot with logarithmic x-axis
stem(f, abs(dft_data))
set(gca, 'xscal', 'log')
ax = gca;
ax.XTick = f;
xlabel('Frequency [Hz]')
title('DFT Magnitude')

I'm not sure what these magnitudes mean, they don't seem logical based on my input.

I'm not sure if I need to divide by the sampling frequency Fs or the number of data points N, but I don't see how that would make the data more meaningful. I have the signal processing toolbox available to me, but I'm not proficient with it.

Any guidance is appreciated!

Answer 1

[Star Strider]

I haven’t used goertzel in a while, and at the time I was interested in the relative amplitudes, not the ‘correct’ amplitudes. Looking through the Proakis reference and the Wikipedia article, it seems that the output is the unscaled result of the recursive filter.

As with the fft (link), it is likely necessary to normalise the raw goertzel amplitude results by the length of your input vector, specifically here:

dft_data_scaled = 2*dft_data/numel(SPL);

From my bit of research, that should give you the correct approximate amplitudes.

[Shaun Gair]

Thanks much for the answer. When I put that code in the magnitudes certainly SEEM correctly scaled. I'll have to do some more analysis to confirm!

[Star Strider]

My pleasure!

The factor of 2 exists because complex transforms have complex-conjugate representations (positive and negative frequencies), while the plot displays only the positive half of the spectrum. The 2 corrects for it having only half the signal energy in half the plot. Dividing by the length of th signal vector has to do with the algorithm summing over the entire vector, so it is necessary to divide by the length of the vector to normalise for that. (I used numel here assuming ‘SPL’ is a vector. If it is a two-channel (Nx2) matrix, use the length function instead.)

The first example in the goertzel function documentation supports this. The amplitude of the frequencies presented to the function are 1 with a signal length of 205. The magnitude (absolute value) of the goertzel output at those frequencies are [99.18,103.6]. Multiplying by 2 and dividing by 205 yields a ‘corrected’ amplitude of [0.9676,1.011], acceptably close to 1 for both.

Answer 2

[KABOURI MOUAD]

1) Implémenter une fonction Matlab, goertzel.m, qui lors de la réception d'un vecteur d'échantillons,

renvoie l'harmonique fondamentale desdits échantillons calculée par l'algorithme de Goertzel.

La fonction doit s'adapter à n'importe quelle taille du vecteur d'entrée.

2) Implémentez un programme, hw3.m, dans lequel vous comparez que, étant donné un vecteur d'entrée de

10 échantillons entiers dans l'intervalle de 0 à 100 choisi par vous, la fonction de la section précédente donne

le même résultat que le composant correspondant renvoyé par La fonction fft de Matlab.

[Star Strider]

Non!