Introduction
The current analytical digital processing methods of signal frequency analysis uses Correlation, Fast Fourier Transform (FFT), and/or the Goertzel Algorithm to convert signal data to frequency data. This is a discussion of the algorithms and a comparison of the response of the FFT, Goertzel, and Frequency Convolution, both for center frequency and off center frequency, which shows the signal response and frequency results to be equivalent to a high degree of accuracy. Frequency measurement result scaling, as related to the input signal level, is also discussed.
Convolution, Correlation
Convolution, sometimes referred to as correlation, compares input data with a known data pattern. This can be done in single dimension such as a sound signal, done in 2 dimensions such as an image, or even done with 3 dimensions such as multiple viewpoints used to correlate multiple images to physical object such as an airplane viewed from multiple positions. An example of looking for an item on a 2 dimensional image, the target shape can be sized to an image item of interest to compare and is correlated by taking each pixel of the sized target shape and multiplying the image pixel and summing the results in the target image area to get a correlation result. The higher the result, the better correlation. There may also be other considerations such as rotation of target shape, etc.
Frequency Convolution
For frequency convolution analysis, the data is convoluted with a unit amplitude sine waveform of a fixed target frequency to determine the target frequency content of the data. The target convolution frequency can be any arbitrary frequency between 0 and the Nyquist frequency (half the sample rate frequency). Doing multiple independent convolutions will allow multiple target frequencies to be measured.
Since the data signal is probably not in exact phase with the target convolution frequency, the data must also be convoluted with another sine wave form, 90 degrees out of phase with the target sine wave, which is the Cosine wave. Correlations with these two relationships where the data has some arbitrary input phase, are often referred to as the real and imaginary components of the data. From tradition, the polar complex form notation for data z is z = r(cos φ + i sin φ), where r is the circular polar radius, φ is the position phase and i stands for the imaginary component. This is particular definition is unfortunate since this implies that "real" signals start at full amplitude at 0 phase. Natural data always starts at a 0 amplitude at the signal's start and proceeds to cycle to its maximum amplitude over time while the definition implies to start evaluation with a Cosine wave. Starting at full amplitude, introduces an impulse into the data, which an impulse has frequency components that spread across the frequency spectrum, which results in the addition of noise at the other frequencies when trying to measure the signal's frequency components. For these discussions the polar complex form definition will be ignored and the primary signal of interest for the primary correlation will be the Sine wave since it starts at 0 amplitude at a signal at 0 phase and has direct correlation with the unit Sine wave.
The start of a detailed discussion of the frequency convolution begins with a plot of the input data shown in Figure 1: Convolution Input Data. In this case, the input data is exactly the same 500 Hz frequency as the target convolution frequency and starting with a Sinusoidal wave at zero phase and with an arbitrary amplitude of 100. Data sample n is d_{n} = A*sin(2π(F_{d}/F_{sr})) * n ), where n is the sample number starting at 0, amplitude A = 100.0, data frequency F_{d} = 500 Hz, sample rate frequency F_{sr} = 8000 Hz.
Figure 1: Convolution Input Data
For the Sine wave portion of the convolution, we start with the unit Sine wave at the target frequency of 500 Hz shown in Figure 2: Convolution, Unit Sine Wave at Target Frequency. The Unit Sine Wave sample n is S_{n} = sin(2π(F_{t}/F_{sr}) * n ), where target frequency F_{t} = 500 Hz. Each sample of the input waveform is multiplied by the corresponding waveform of Unit Sine Wave, s_{n} = d_{n}*S_{n}, and plotted in Figure 3: Convolution, Input Data Multiplied by Unit Sine Wave. Each of these results are added into the running sine wave sum total t_{sn} = t_{sn} + s_{n} and the total result accumulates increasingly as more data samples are added as shown in Figure 4: Convolution, Running Sum of Sine Convolution.
Figure 2: Convolution, Unit Sine Wave at Target Frequency
Figure 3: Convolution, Input Data Multiplied by Unit Sine Wave
Figure 4: Convolution, Running Sum of Sine Convolution
The Cosine Convolution occurs in the same way using the unit cosine samples C_{n} = cos(2π(F_{t}/F_{sr}) * n ), plotted in Figure 5: Unit Cosine Wave at Target Frequency. The input data is multiplied by the unit cosine, c_{n} = d_{n}*C_{n} and plotted in Figure 6: Convolution, Input Data Multiplied by Unit Cosine Wave. The running cosine wave sum total t_{cn} = t_{cn} + c_{n} is plotted in Figure 7: Convolution, Running Sum of Cosine Convolution. Since input data is in phase with Unit Sine wave but 90 degrees out of phase with the cosine wave, the correlated function is a wave with both positive and negative areas that results in the convolution accumulation sum oscillating between 0 and in this cases 120.7107.
Figure 5: Unit Cosine Wave at Target Frequency
Figure 6: Convolution, Input Data Multiplied by Unit Cosine Wave
Figure 7: Convolution, Running Sum of Cosine Convolution
The complete convolution is the trigonometric sum of the accumulated sin and cosine convolution results, which is done by summing the squares of each of the results, m_{t}^{2} = t_{sn}^{2} + t_{cn}^{2} from the sine and cosine wave convolutions to get the magnitude squared is shown in Figure 8: Running Convolution of Relative Power. The result is the square of the amplitude (or magnitude) which is directly related to power without the resistive constant seen in a circuit (or resistive constant of 1.0) where power P = V^{2}/R and V is the instantaneous voltage amplitude. Since signal power is generally the desired information and converting back to amplitude level requires calculating the square root, which is a very expensive operation, the square root operation is generally not done in high performance signal processing systems, if it can be avoided. However the accumulation of the signal during the convolution is linear (excluding the ripple) in the amplitude scale, which must be used to understand how the correlation sum result is accumulating.
Figure 8: Running Convolution of Relative Power
Frequency Convolution with Phase Shifted Input Data
When the input signal is not in phase with the unit sine wave, the sine and cosine convolutions have different results. For the case of having the input signal shifted forward by 36 degrees, the result sine convolution t_{sn} is shown in Figure 9: Convolution, Running Sum of Sine Convolution, 36 Degree Phase Shift and cosine convolution t_{cn} is shown in Figure 10: Convolution, Running Sum of Cosine Convolution, 36 Degree Phase Shift. Both of these have an upward ramp and the squares sum up to a waveform very similar to Figure 8: Running Convolution of Relative Power, however there are some small result differences which will be covered later in this document.
Figure 9: Convolution, Running Sum of Sine Convolution, 36 Degree Phase Shift
Figure 10: Convolution, Running Sum of Cosine Convolution, 36 Degree Phase Shift
Magnitude and Accuracy Increase with Convolution Size
The convolutions plots shown so far were done with very few samples so that the effects that is occurring from point to point could be easily seen across the plot for a thorough understanding of the frequency convolution process. Although this is showing all of the sample points of the process, except for very complex high accuracy processing, it is generally the last data point at the end of the data block, which is the last point on the graph, that is the target quantitative frequency power result of the processing. The convolution result of adding more processing data points from doing just 41 to 200 data points is plotted in Figure 11: Running Convolution of Relative Power, 200 Samples. We see that the curve is a similar exponential shape but with smaller ripples.
Figure 11: Running Convolution of Relative Power, 200 Samples
The magnitude, which is square root of the power, is plotted in Figure 12: Running Convolution of Magnitude, 200 Samples. This is almost a straight line with a little ripple in it, showing that the process has linear amplitude relationship. This strait line indicates that this process can go on to an indefinitely larger data set sizes and get valid center frequency results. However the ripple can be thought of as process error which affects the end result measurement accuracy depending on where the last point of the output is calculated. The relative ripple remains the same, but as more samples are added, the ripple error becomes less significant compared with end magnitude result and end measurement accuracy goes up. One way to improve accuracy is to determine the best fit line through the data and use that last point intersection of that line as the result measurement. Another is to just process more data points for a higher result data to the ripple ratio value decrease, for more center frequency accuracy. However the number of samples processed greatly affects the off-center frequency results, which will be discussed later.
Figure 12: Running Convolution Magnitude, 200 Samples
In some cases there are data points that have the exact result. At the π point of each frequency Sine and Cosine convolution, where Sine waveform crosses the 0 Axis, the results are exact, however there are only a small number of frequencies that these points will occur. The example sample rate and frequency was chosen so that these points would be present in the data set. For 500 Hz sine wave at a 8000 Hz sample rate, the π point of Sine waveform is at every 500/8000/2 = 8 data points. As can be seen in Figure 2: Convolution, Unit Sine Wave at Target Frequency the Sine wave starts at data point 0 and π sine wave points occur at data point 8, 16, 24, 32 and data point 40. At data point 200, 8 * 25 = 200, the exact result value present and will be used in scaling discussions. In reality, the convolution ending in an exact π position point result point for an arbitrary frequency is not common, unless the frequencies are chosen by design, resulting in most convolution data measurements having a ripple related error.
Normalization
Normalization is not being applied to help make it clear what is happening and to show how the data relates across algorithm types. Normalization is simply the end result divided by the number of samples processed (1/n) and is applied to the amplitude or (1/n)^{2} is applied to the relative power result. This allows the same result, outside the numerical process errors, phase related impulse errors, and process variations for different convolution lengths, for exact input signal frequency measurements of uniform center frequency data.
Since division is a complex algorithmic process that takes a series of shifts and subtractions, the result normalization divides always taking considerable number of CPU clock cycles to complete. However most processors can do a full multiply operation in 1 CPU clock cycle. The 1/n or (1/n)^{2} is calculated once during initialization and normalization is applied by using a multiply operation.
Convolution Result Scaling
When we normalize the amplitude of the exact result data point 200 of Figure 12: Running Convolution Magnitude, 200 Samples we get 10,000/200 = 50.0, which is exactly half of the input wave amplitude of 100.0. Where does this result scaling come from?
The convolution is combining the signal's amplitude with the unit Sine wave amplitude. However the signals amplitude changes through time for a sinusoidal wave has an average amplitude of 1/(square root of 2) = (1/(2^{1/2})) = 0.7071... called the RMS (Root Means Square) value. Multiplying the two average amplitudes together = (1/(2^{1/2}))* (1/(2^{1/2})) = 1/2 = 0.5. To convert the convolution result to the input signals RMS value, multiple the result by square root of 2 = (2^{1/2}) = 1.414... . To get the input signals amplitude value, multiple by 2. The RMS value is the closest value related to actual time averaged signal power.
Alternatively the unit convolution Sine and Cosine waves can be scaled up by the RMS value or by 2 to get an end amplitude RMS or amplitude result at the desired signal related result. However the unit amplitude response is more easily related to the Goertzel and FFT response results.
Process Ripple Error
The ripple observed in Figure 12: Running Convolution Magnitude, 200 Samples, is related to the area of the multiplied/convoluted sine wave curve as the sequential convolution result is being added in. As seen the Sine convolution plotted in Figure 3: Convolution, Input Data Multiplied by Unit Sine Wave, the convolution input varies across time like a sinusoidal wave (but is not a true sinusoidal wave) and the addition of the peak area is accumulating at a higher rate while the trough areas is accumulating at a lower rate, resulting in a rippling effect in the convolution result. Viewed from just the data input of power accumulation, when the input signal is close to zero the signal power is close to zero which will accumulate little power for the result as these samples are added in. When the signal is at the peak amplitude, the per sample accumulation of the signal power is at its maximum. The bottom line is that the ripple is inherent in the frequency measurement of a sinusoidal signal over time and not some general algorithm processing error. For general and lower quality measurements, this ripple can be ignored, while high precision measurements, the ripple may need to be considered.
Subtracting out the exact response from the convolution we get the ripple plot in Figure 13: Ripple of Running Convolution Amplitude. The ripple passes through the zero point at each π radian location, which is at every 8 sample result (0, 8, 16 ... ) location. The ripple is highly biased on the positive side compared to the negative side. Eliminating the ripple error by calculating the signal based on a line through the average would result in an error, which might indicate that the exact result at π radian locations might not be correct, however other factors such as the math has an exact matches for scaling and the line needing to meet at the 0,0 axis locations indicates that the 0 plot position is correct. Looking at the plot, the first portion of the ripple is notably different than the later portions. A difference plot between successive ripple amplitude is shown in Figure 14: Change in Amplitude Ripple Difference for Running Convolution. The biggest change is at the beginning and the ripple quickly becomes more uniform. however the ripple takes considerable time to dampen out as can be seen in the up-scaled plot in Figure 15: Up-Scaled Amplitude Ripple Difference for Running Convolution. Correcting for result ripple error is not strait forward.
Figure 13: Ripple of Running Convolution Amplitude
Figure 14: Change in Amplitude Ripple Difference for Running Convolution
Figure 15: Up-Scaled Amplitude Ripple Difference for Running Convolution
The example input data is data has no data noise, other than the very small numerical noise from limited binary size that will cause some small digital related error in the process results. In the same way, processing more samples would make the end result less affected by ripple noise, so processing more data samples might be more desirable when the input frequency and target measurement frequency are exactly the same. When the data contains data noise or other large data components not at the target frequency range, along with the target center frequency signal, the measurement result of the correlated signal to non-signal noise improves with the longer convolutions. This convolution process can be done with an arbitrary number of data samples, the real question comes down to how many samples should be processed to get an adequate result for the data measurement goals?
Input Signal Phase Shift Result Difference
There are result differences of a convolution output when the data input data phase is different than the unit sine convolution wave, which is substantially much more probable for actual data to not be at a 0 degree phase, than to be at exactly 0 degree phase. Going back to the result differences for 36 degree case used to show the change in the sine and cosine convolutions discussed earlier. The difference in relative power output results are plotted in Figure 16: Convolution Difference of Magnitude^{2} for Phase Change, which is ramping up as the convolution proceeds. However the amplitude plot in Figure 17: Convolution Difference of Magnitude for Phase Change plot shows that the error remains in fixed range as the net convolution amplitude result in Figure 12: Running Convolution Magnitude, 200 Samples is ramping up. The relative phase related difference becoming less significant as more data samples are processed. The 500 Hz target convolution frequency is an integer divisor to the 8000 Hz sample rate, resulting in data points at each half cycle which is the π radian points. It turns out that the difference in processing between 0 and 36 degree phase data start, shows no difference at interval half wavelengths which can be seen in the difference plots where they go down to zero. All other partial wave forms points has some potential measurement difference associated with them. There are only a few frequencies that have an integer number of samples in a complete cycle that would have these zero difference levels. From a different prospective, we can say that there is generally a phase measurement error associated with a non zero phase difference at the start of the convolution, with the bad point being that this error changes depending on where in the partial cycle that the data measurement starts, assuming non π radian multiple process span.
Figure 16: Convolution Difference of Magnitude^{2} for Phase Change
Figure 17: Convolution Difference of Magnitude for Phase Change
All other processing methods with equivalent processing characteristics will have the same phase offset processing errors, not that "it is not a phase caused error", since they will all have the same frequency measurement results. In the non-center frequencies, this can be considered at phase offset step at the first sample which may have an impulse and/or step type effects of spreading a noise response across all the convolution frequencies.
Frequency Convolution with Off Frequency Data
There may be some cases that the exact convolution target frequency is present, but more probably contains data with a frequency close to but not centered on the exact convolution frequency. In fact there is a many types of data that changes frequency across the data samples rather than have a mono tone across the signal time, such as voice. This is to give a feel for this effect and also happens to be a big source of potential data measurement errors, and can be the primary source of the error in some data measurement cases.
The off frequency response effects can be observed by convolution of input data with a frequency that is different than the target convolution frequency. The ongoing convolution process of 600 Hz input sine wave is fed into a 500 Hz frequency convolution is shown in Figure 18: 500Hz Frequency Convolution with 600 Hz Input Data. The result power level starts to increase as power is integrated in time, until the phase of the input signal input phase changes enough and goes out of phase with the convolution frequency and starts to subtract. However since we are convoluting with both a sine and cosine wave, the out of phase change for these occurs at different sample points resulting in both addition and subtraction for the same input data points, which results in a both up going and down going fluxuations seen across the waveform. Eventually the energy subtracts and the results goes back to zero (not always zero if waveform is not an integer sample count and/or also contains accumulated processing error) and the data starts adding again for another in cycle. This particular frequency produced 2 1/2 cycles during the 200 sample convolution, but would continue cycling until the input data stopped or the convolution ended.
Figure 18: 500Hz Frequency Convolution with 600 Hz Input Signal
A different arbitrary input frequency will have a different number of samples before going back to zero. Figure 19: 500Hz Frequency Convolution with 539.6666 Hz Input Signal shows the result for an arbitrary 539.6666 Hz frequency closer to the convolution frequency of 500 Hz, which approaches zero at 200 samples rather than about 80 samples. As the data input frequency moves away from the target center frequency, these peaks occur in shorter and shorter number of samples and at a lower amplitudes.
Figure 19: 500Hz Frequency Convolution with 539.6666 Hz Input Signal
Measuring data signal power off exact center frequency becomes problematic since the end measurement result will depend on where on the curve that the process ends. Anything pass the peak is problematic since data starts to subtract, and at that point there becomes multiple logical locations of the data measurement to relate back to the actual input frequency. Is there a low power signal at a close frequency or a higher power signal at frequency further from the center frequency? There is frequency information in the way the curve changes over time, but convolution (and FFT) generally are done to arrive at the end point representing the "frequency content" of the data at the particular convolution frequency points.
An important effect of this oscillating result indicates there is limited frequency range around the center frequency that is useful, which changes with the number of data samples processed. The number of processed samples and to the overall frequency response needs to be considered. Often, mistakes are published that indicate that the number of processed samples can be doubled and get the same result with just a scaling change, which except at the center frequency, all other frequency measurement results will have a different end result. This also requires that all data for measurements be blocked into a limited sample block length size rather than processing an arbitrary long data block.
The effect of the convolution going all the way down to zero relates an important property of this method. The output wave point going to zero indicate that the convolution has a zero frequency width measurement. This also indicates that mathematical relationships using change of phase relationship between input frequency and the convolution frequency across samples (time) can be used to derive the expected response. However, since in many cases real data has varying frequencies, this creates some issues on trying to do accurate measurements using zero frequency width functions on real world data.
Goertzel Algorithm
For single or small frequency count frequency measurement, or specific exact center frequencies not part of the FFT center frequency set, the Goertzel Algorithm is commonly used.
The Oscillating Function
There exists a simple oscillating function, which appears to be oscillating around a circle on the z-plane (perpendicular to the x-y plane) at a set frequency, that uses the previous output two values to generate the next sample value in the series.
B = 2cos(2πF_{c}/F_{sr}) | Constant for oscillation Center Frequency F_{c} |
φ = 2πF_{c}/F_{sr} | Sample step angle, F_{sr} is Frequency Sample Rate |
ѱ = 2πS_{a} | Starting angle in radians in S_{a} portion of 2π, 0 for 0 phase start |
n_{-1} = Asin( -φ + ѱ) | Initial Value for previous sample n_{-1}, scaled by Amplitude A |
n_{-2} = Asin( -2φ+ ѱ) | Initial Value for 2 samples before: n_{-2}, scaled by Amplitude A |
n = Bn_{-1} - n_{-2} | Calculate next data point n, repeated for series of points |
At start of next data point, n becomes n_{-1} and n_{-1} becomes n_{-2} |
The function is initialized for frequency, amplitude, and the two previous initial points are calculated. The phase ѱ is set to 0.0 for 0.0 phase sine wave, however the starting points can be adjusted for phase, where 90 degree (0.25*2π) phase increase will result in a cosine wave. Consecutive points can calculated continuously and indefinitely, but real systems have limited binary accuracy and errors will creep in over a very large sequence count and output will deviate from ideal oscillation. However is generally not an issue in short time frames with the accuracy of most current CPUs, with some compute architecture exceptions.
With some care, this oscillation function as well as the full Goertzel can done using a fixed point processor assigning a decimal point at some well chosen bit location, ideally 32 bit CPU, however the rate of result and error accumulation and the offset of the center frequency due to the truncated B value needs to be carefully considered. Ideally full 64 bit floating point (double) math should be used when possible but fair results can be obtained with 32 bit floating point (float) math.
The F_{c} center frequency and F_{sr} sample rate frequencies can be any arbitrary decimal number and does not need to be an integer value. For example F_{c} is 539.6666 Hz is used in some examples.
A plot of using the formula is shown in Figure 20: Simple Oscillating Function, which is a 100 amplitude, 500 Hz sine wave at a 8000 Hz sample rate starting with a 0 degree phase. This generates a sinusoidal result, but have not been able to find the mathematical proof relating this simple math operation to a circular function.
Figure 20: Simple Oscillation Function
The Data Correlation Operation
The next step is to match the input signal into the oscillating function. This is done by just adding in the data sample for each calculation, and the formula becomes:
n = Bn_{-1} - n_{-2} + d | Calculate next data point n, also add in the next sample data d |
B = 2cos(2πF_{c}/F_{sr}) | Goertzel Constant for Center Frequency: F_{c} and at Sample Rate Frequency: F_{sr} |
n_{-1} = 0.0 | Initial 0.0 value at start of Goertzel Algorithm operation |
n_{-2} = 0.0 | Initial 0.0 value at start of Goertzel Algorithm operation |
At start of next data point, n_{-1} becomes n_{-2} and n becomes n_{-1} | |
mag^{2} = n^{2} + n_{-1}^{2} - Bnn_{-1} | Calculates the Magnitude^{2} result of the Goertzel response |
The n position values are before the start of the next data point |
This simple addition of the input data to the oscillation function, is the Goertzel Algorithm.
Additional Goertzel supporting formulas:
imag = n_{-1}sin(2πF_{c}/F_{sr}) | Imaginary Axis value (x axis) of the circular function |
real = n - n_{-1}cos(2πF_{c}/F_{sr}) | Real Axis value (y axis) of the circular function |
The n position values are before the start of next data point | |
phase = atan2(real, imag) | Radian angle in range of -π to π |
Rotating plot is clockwise, 0 phase is Real Axis (positive y axis) | |
Starts at 0 phase for Goertzel doing first non-zero data sample | |
Positive increase in phase for increasing sample count | |
a = Scale*Norm*sqrt(mag^{2}) | Amplitude with Norm = 1/N Samples processed |
a = Scale*Norm*sqrt(real^{2} + imag^{2)} | Alternative Amplitude Calculation |
Scale | Scaling of 2^{1/2}: RMS value, 2.0: Signal Amplitude, and 1.0: Un-Scaled |
Increasing phase for a plot rotating in a counter clockwise can be achieved by using atan2(imag,real) which 0 phase starts on x axis with the real value plotted on the x axis and the imaginary part plotted on the y axis. I think of the circle rolling forward on the positive x axis, but is a matter of prospective.
Goertzel Algorithm Center Frequency Data Processing
A 500 Hz, 100 Amplitude, 0 start phase sine wave was generated using the described oscillating function and is plotted in Figure 20: Simple Oscillation Function. This was fed into the Goertzel Algorithm with the same center frequency, with the real and imaginary components calculated at every processed sample and plotted in Figure 21: Goertzel Process of Same Frequency, 100 Amplitude Sine Wave. The oscillation starts at the center point and the oscillation grows in radius (amplitude) as the a sine wave with the same frequency is processed by the Goertzel Algorithm. Calculating the Magnitude^{2} and plotting in Figure 22: Goertzel Magnitude^{2}, Same Frequency, 100 Amplitude Sine Wave shows an exponentially increasing value, but also includes ripple as described in the convolution section.
Figure 21: Goertzel Process of Same Frequency, 100 Amplitude Sine Wave
Figure 22: Goertzel Magnitude2, Same Frequency, 100 Amplitude Sine Wave
Calculation of the clockwise rotation Goertzel result phase in the range of -π to π of each data point is plotted in Figure 23: Goertzel Center Frequency Processing Phase. The input signal is a sine wave at equally spaced phase angles. There is noticeable non-linearity in the phase angle processing of the first sine wave of 0 to 16 samples as the Goertzel is shifting from a 0 phase at first sample to a 90 degree lag at sample 16 + 16/4 = 20 as the processing proceeds. Later phase changes deviations are much lower and the equal phase change error goes down as more center frequency samples are processed and creates visually strait lines in the -π to π plot regions at lower phase errors. The primary error is at the beginning as was the case for amplitude ripple in Figure 14: Change in Amplitude Ripple Difference for Running Convolution.
Figure 23: Goertzel Center Frequency Processing Phase
The input signal data, the Goertzel response pattern, and the accumulating signal are plotted together in Figure 24: Plot of Data, Goertzel Response, and Result Amplitude. The amplitude of the Goertzel increases with every data input sample and very quickly grows in amplitude, quickly dwarfing the input signal. The result measurement amplitude grows at a rate less than half the Goertzel amplitude increase. Although the number is large, there is an upper bound to the number samples that can be processed by a Goertzel algorithm, where the input data value relative to the Goertzel amplitude becomes equal to the digital number's binary accuracy and the data will no longer add in. For the 32 bit floating point number, the mantissa has a much lower number of bits and this results in a much smaller number of samples to reach this limitation.
Figure 24: Plot of Data, Goertzel Response, and Result Amplitude
Increasing the processing to 200 samples, the circular Goertzel plot is shown in Figure 25: Goertzel Process 200 Samples, 100 Amplitude Sine Wave which shows the Goertzel response amplitude for input of the exact frequency, continuing to grow as more data is processed. The phase change is constant and the plot points radiate outward in a equally phase spaced strait lines from the center. Note: The calculated values on the real y axis are very close to zero (one example: -2.86491E-11), it is the Excel spreadsheet plot software that is producing the visual offset from the y axis.
Figure 25: Goertzel Process 200 Samples, 100 Amplitude Sine Wave
The result Magnitude^{2} plotted in Figure 26: Goertzel Magnitude2, 200 Samples, 100 Amplitude Sine Wave shows the result growing as more data is added. When comparing to this to the frequency convolution plotted in Figure 11: Running Convolution of Relative Power, 200 Samples it shows that the plots are the same. Used a Excel spreadsheet to do the algorithms and generate plots and found that all values up to the 10 digit values in the spread sheet were exactly the same. Doing a difference to the numbers, the majority had a zero difference, but the numbers that had a difference got larger as numbers grew with largest and last non zero difference occurring at data point 191 which has a value of 92,160,000.0 with a difference of 1.341e-7, which would be about 14 digits of equivalent accuracy, which is acceptable process error range for 64 bit floating point operations. The convolutions was done doing using sine and cosine functions for the unit waveforms with input data generated using a sine function with a vector addition to get the result as described in the convolution section. While the Goertzel processing was using the oscillating function described above to generate the same input, and then used the Goertzel Algorithm to process the inputs and the Goertzel magnitude^{2} formula to generate the power response.
Figure 26: Goertzel Magnitude^{2}, 200 Samples, 100 Amplitude Sine Wave
This shows the unit Frequency Convolution and the Goertzel Algorithm, have equivalent response for the same input data starting with 0 phase sine wave having the same frequency as the measurement center frequency. The un-scaled Goertzel amplitude plot is shown in Figure 27: Goertzel Magnitude, 200 Samples, 100 Amplitude Sine Wave, with an amplitude of 10,000 and data point 200. Normalizing the end result gives 10,000/200 = 50.0, which is 1/2 the input signal amplitude. Multiplying by 2^{1/2} will give the RMS value, multiplying by 2.0 will gives amplitude of the input signal with a uniform amplitude. The result ripples are present and results match exactly to the unit frequency convolution results.
Figure 27: Goertzel Magnitude, 200 Samples, 100 Amplitude Sine Wave
A important effect for the Goertzel is that if the input data goes to zero during the processing, the Goertzel output level will remain constant and continue to oscillate at the center frequency and maintain the frequency's phase relationship to the original signal portion as it continues to oscillate, which would plot a circle on the Goertzel circular z-plane plot. From a data measurement prospective, the first portion of the data is measured and the measurement is maintained until the end of the data block being measured.
Goertzel Algorithm Processing Off Frequency Data
In natural data and in most cases the off center frequency processing is going to be more common than the exact center frequency data processing. Having an idea of how the Goertzel Algorithm responds to the off center frequency data is very desirable.
Figure 28: 500 Hz Goertzel Response to 600 Hz, 100 Amplitude Sine Wave shows the Goertzel algorithm beginning to accumulate the signal at the start. It accumulates and forms a peak then the starts to go back down as signals phase goes out of phase with the oscillating phase of the Goertzel algorithm at its center frequency. The response goes all the way down to zero (0.0) to fully cancel the response to the out of phase signal before starting to accumulate the signal to the new Goertzel algorithm response to the data, and the response pattern repeats.
Figure 28: 500 Hz Goertzel Response to 600 Hz, 100 Amplitude Sine Wave
This example has approximately 2 1/2 cycles of the repeating response. The z-plane plot shown in Figure 29: 500 Hz Goertzel Z-Plane Response Plot of 600 Hz, 100 Amplitude Sine Wave also contains the 2 1/2 cycles of the response pattern. In this case since the response goes to zero, the 2 1/2 cycles overlays one another and forms one pattern.
Figure 29: 500 Hz Goertzel Z-Plane Response Plot of 600 Hz, 100 Amplitude Sine Wave
Moving closer to the Goertzel algorithm center frequency, Figure 30: 500 Hz Goertzel Response to 539.6666 Hz, 100 Amplitude Sine Wave we see a response that peaks in the middle and goes down to zero close to the end of the 200 sample process. The zero point is not present in the data points, but a zero data point would lay somewhere between 2 samples points when the processing passes through this transition, and in this case is between the last 2 data points, which have roughly equal values.
Figure 30: 500 Hz Goertzel Response to 539.6666 Hz, 100 Amplitude Sine Wave
The un-scaled non-normalized amplitude is shown in Figure 31: 500 Hz Goertzel Amplitude of 539.6666 Hz, 100 Amplitude Sine Wave. The amplitude response curve starts responding, peaks out, and goes back to zero. The first portion is can be thought as somewhat linear but quickly goes non-linear as we move along the process curve. At the center frequency, the Goertzel Algorithm has a linear amplitude response (disregarding the sinusoidal introduced ripples), however all off frequency response is affected by this signal to center frequency phase relationship that varies in length depending how far off frequency the signal is. The off frequency measurement is a non-linear process and as you move away from the center frequency it becomes more non-linear.
Figure 31: 500 Hz Goertzel Amplitude of 539.6666 Hz, 100 Amplitude Sine Wave
The response shown in Figure 28: 500 Hz Goertzel Response to 600 Hz, 100 Amplitude Sine Wave which is a 100 Hz off frequency from the center frequency, has a magnitude^{2} response of 1,907,280.8 at the end of the 200 sample measurement. In Figure 30: 500 Hz Goertzel Response to 539.6666 Hz, 100 Amplitude Sine Wave the input Sine wave is off frequency by 39.6666 Hz and at the end of the 200 sample measurement has a end result of 6557.9, much smaller than the measurement result more than twice the frequency distance from the center frequency.
This effect creates frequency power measurement results that can cancel at some frequencies and add at other frequencies depending on the data content, which usually is the case for natural data that the content is not strictly controlled as the case of a natural sound source picked up by a microphone. Increasing the number of samples processed to get to the end result only narrows the off frequency distance that these patterns occur and off frequency addition and cancelation exist, but the problem persists. This effect is often observed in the FFT where many frequencies are being measured, and some patterns sometimes appear giving the impression that there are more low signal level data frequency components in the data content than is actually present in the data since a lower value response is often between the response and the response closest to the actual data frequency. From a strict point of view, this can be viewed as process measurement noise. However since the patterns have algorithmic type relationship, these correlated off frequencies result patterns could be used to measure the signals content more accurately in some cases, as well is important in an inverse frequency transformation. Just because the off frequency results appear and inverse frequency transformation works does not mean that all these frequency components actually exist and must be at the same levels on a real linear signal power frequency scale. Analyzing data with zero frequency width algorithms, has it's draw backs.
The zero Goertzel response value will only be present in the response data when the sample rate, Goertzel center frequency, and data frequency has a sampling relationship. In this case, the last 2 points are 6581.8 and 6557.9 and the zero response point is approximately halfway between these 2.
Goertzel Algorithm with Off Frequency Data Relationship to Frequency Convolution
The same off frequency data input cases for the Goertzel and Frequency Convolution along with the same process length was used to allow a direct comparison between the Frequency Convolution and the Goertzel Algorithm responses. The 2 off frequency responses for the Frequency Convolution is shown in Figure 18: 500Hz Frequency Convolution with 600 Hz Input and in Figure 19: 500Hz Frequency Convolution with 539.6666 Hz Input Signal. The 2 off frequency responses for the Goertzel Algorithm are shown in Figure 28: 500 Hz Goertzel Response to 600 Hz, 100 Amplitude Sine Wave and in Figure 30: 500 Hz Goertzel Response to 539.6666 Hz, 100 Amplitude Sine Wave. If you look at these 2 sets closely, you will see that the response curves look almost identical in shape and amplitude.
All of these operations was done in an Excel Spreadsheet and can be numerically compared. Even though the frequency convolution was done using all Sine and Cosine functions to generated the data and unit waveforms and Goertzel algorithm which used the oscillation function above, the results agreed to 13 bits of accuracy, except for the zero response points that had similar noise level had an accumulated process noise close to 0.0 of a 1e-23 noise range which noise has no related digits of accuracy. This data implies that the frequency convolution and Goertzel algorithm has the same frequency measurement characteristics for both at the center frequency and for off frequency data. This implies that there might be a mathematical way to prove this response equivalence, but is not a strong interest or on my agenda to do so.
At this time, I have no supporting proof for this theory, however due to the very complex nature of the off frequency responses being equal to a high degree for accuracy, this provides a very strong indication that the Frequency Convolution and Goertzel Algorithm have identical frequency response characteristics, even though their mathematics are very different and one has the static waveforms and the other has a dynamic oscillating function.
FFT (Fast Fourier Transform)
The FFT is the most used algorithm for doing frequency analysis due to its process efficiency. It does a frequency analysis of equally spaced frequencies with a frequency count of half of the 2n sample count block size.
The following code generates the FFT frequencies for nFreq frequencies which is half the FFT length for the frequency array pFreq. The first frequency is 0.0 Hz, which is the DC value and nFreq = 512 for a 1024 point FFT.
double Bin_Width = ((double) Sample_Rate)/( 2 * nFreq );
for ( i = 0; i < nFreq; i++ ) *pFreq++ = Bin_Width * i;
The FFT came about from needing to do many sine and cosine wave convolutions as efficiently as possible when computers were very expensive and very slow. It was noted that formulas repeated in a power of 2 data size counts allowed sine and cosine formulas to be easily combined together and take advantage of symmetry and the repeating nature of sine and cosine waveforms. It was also found that in some cases the multiplies and sum operation could be combined. The work on the combination of the mathematics to get the minimum calculations possible resulted in the creation of the FFT.
A sine and cosine convolution at the FFT bin center frequency can be used to get the same result as an FFT for the same data sample count, excluding the accumulated numerical error from different count and order of the mathematical operations. Variations to the FFT, maybe the scaling of the internally generated sine and cosine waves has resulted in some different final result scaling values being used. Two FFT variations were used to make sure that an FFT was working correctly, but the results only were different by their final scaling value.
Although other algorithms like the frequency convolution and Goertzel Algorithm can process arbitrary block lengths, the nature of the algorithmic process used to combine the convolutions, and using symmetry of the property of the sine and cosine waves, the FFT is restricted to only 2n sample count block size and that only generates 2n-1 evenly spaced frequency bins. Due to the internal convoluted manor for combination of equations, it is very difficult to observe what is happening to the data during the process, but is the same process as the frequency convolution, and is why the frequency convolution was covered in such detail when it is rarely used. Whatever discussion applies to the frequency convolution, also applies to the FFT.
Rather than create program to do frequency convolutions to show the equivalence to the FFT an indirect approach is being taken. The frequency response equivalence has been shown between the Goertzel algorithm and the frequency convolution. The approach is to compare the response of the Goertzel response at the same frequencies and sample count as the FFT. To do this code was added to a frequency analysis program to support a 1024 point frequency analysis using the Goertzel Algorithm and then processed the same complex audio block with each to compare the result differences.
Shown in Figure 32: 1024 Point FFT Processing Arbitrary Complex Audio Data Block is an analysis program which is stopped and displaying the FFT results of a particular arbitrary 1024 sample data block in the bottom plot. The top plot is displaying the results of a series of FFTs where each point in the x direction is 1 1024 point FFT and the vertical direction is frequency amplitude with 40 frequency plots displayed. The scroll bar on the left side of the upper plot is selecting the range of frequencies being plotted, which is currently showing frequencies 7020 to 8899 Hz. The black vertical line is showing which data block is being displayed in the lower plot. In the lower frequency plot Red is showing the second or Right channel and green is showing data for the first or Left channel. In the upper plot black is used in place of green since with the higher detail, the pen width is narrow.
Figure 32: 1024 Point FFT Processing Arbitrary Complex Audio Data Block
Shown in Figure 33: 1024 Point Goertzel Processing Arbitrary Complex Audio Data Block as FFT is the same data being processed using the Goertzel algorithm at same frequencies and data block size as processed by the FFT shown in Figure 32: 1024 Point FFT Processing Arbitrary Complex Audio Data Block. The upper plot is showing a slightly different frequency range of 6503 to 8183 Hz. The Frequency plot scale of the FFT is at 928 and the plot scale of the Goertzel algorithm is at 589 since the FFT has been scaled and the Goertzel is un-scaled. The FFT in this case was scaled by 0.5π = 1.57079 which is close to the RMS value of (2^{1/2}) = 1.41421, which is more directly related to signal power measurement.
Figure 33: 1024 Point Goertzel Processing Arbitrary Complex Audio Data Block as FFT
The data for each frequency spectrum at the was copied using the clipboard to a spreadsheet. The 1024 data point processing difference after eliminating scaling difference for the 512 frequencies across the frequency range is plotted in Figure 34: Response Difference between FFT and Goertzel Removing Scaling. This plot agrees with better than 10 digits of accuracy with the worst difference is on the first channel and is located in the top frequency range, while the majority of the results has a higher number of matching digits. All frequencies for the second channel are all much better than 10 digit accuracy.
Figure 34: Response Difference between FFT and Goertzel Removing Scaling
An alternative FFT version with different algorithm coding was compiled and plots are shown in Figure 35: Alternate 1024 Point FFT Processing Arbitrary Complex Audio Data Block. The lower plot is scaled at 589 and upper plots frequency range is 7063 to 8742 Hz. Comparison the second FFT with the with Goertzel results shown in Figure 36: Response Difference between FFT Version 2 and Goertzel Algorithm. The lower frequencies has good comparison up to the higher frequencies where the most error occurred. The higher data only compares with 8 digits with most of the lower data having about 12 digits of comparison (looking at numbers). The FFT in this case had no scaling, so compared directly to the Goertzel values.
35: Alternate 1024 Point FFT Processing Arbitrary Complex Audio Data Block
Figure 36: Response Difference between FFT Version 2 and Goertzel Algorithm
Picking a different audio block that had several major frequency peaks, the FFT (first version) frequency spectrum response is shown in Figure 37: 1024 FFT of 2nd Arbitrary Complex Audio Block. Comparison of the FFT response to the Goertzel Response (not shown) is shown in Figure 38: 2nd Audio Block Difference Between FFT and Goertzel. The results compare better than 11 digits, with this time the worst difference is in the low frequency portion.
Figure 37: 1024 FFT of 2nd Arbitrary Complex Audio Block
Figure 38: 2nd Audio Block Difference Between FFT and Goertzel
Considering the highly complex frequency and unique response nature of the algorithms and high complexity of data within these particular data blocks, for the FFT and Goertzel to have numerical match to such a high degree of accuracy can only mean that the FFT and Goertzel must have identical frequency response and amplitude response levels!
Processing Frequency Response Across the 1024 Sample Data Block
The process of frequency accumulation within the FFT cannot be seen due to the intermixed nature of the FFT algorithm, however progressive sample by sample frequency response in the equivalent Goertzel Algorithm response can be seen in Figure 39: Progressive Frequency Response of 1024 Sample Data Block for Right Channel and Figure 40: Progressive Frequency Response of 1024 Sample Data Block for Left Channel, however are slightly scaled differently. The frequency range scroll bar has been adjusted to show frequency range 6374 to 8053 Hz (these number are rounded to be able to list the frequencies and are not integer frequencies). On the left side, the Goertzel has been cleared and data is processed proceeding to the right. Toward the right side, the Goertzel ends and next 1024 block is started at its zero initial starting values. Looking at these you can see signals accumulating and some data cancelling at other locations.
Figure 39: Progressive Frequency Response of 1024 Sample Data Block for Right Channel
The non-linear nature can be observed in Figure 39: Progressive Frequency Response of 1024 Sample Data Block for Right Channel. In the middle of the plot, the top purple 7967 plot is reaching its maximum, but drops by about 50% in the second half since the signal data has become out of phase with the center frequency data of Goertzel Algorithm. About half way through the second half the data ends and the Goertzel remains somewhat steady other then more minor interactions from other data components. The 1024 end result for this would be about 30% the current scale. However, if 512 point FFT or Goertzel was performed, the first result would be at about 60%. Doing a second 512 point Goertzel, the data in that frequency range has not been accumulated, and the data would add until the signal drops ends about half way though result in an output at somewhere around 30%. The result of the 1024 point FFT would be 30% while the 512 point FFT would get 60% on the first block and about 30% on the second block for a total of 90%.
Figure 40: Progressive Frequency Response of 1024 Sample Data Block for Left Channel
The bottom line is that the power measurement processing is non-linear for non-center frequencies and the data does not always accurately reflect the actual frequency contents of the data blocks and some care must be exercised in interpreting the results.