Fix errors in octave frequency scale

Issue #332

src/AudioAnalysisTools/DSP/OctaveFreqScale.cs

@@ -90,7 +90,7 @@ public static FrequencyScale GetStandardOctaveScale(FrequencyScale scale)
             var powerSpectrogram = ConvertAmplitudeToPowerSpectrogram(octaveScaleM, windowPower, sampleRate);
 
             // Convert the power values to log using: dB = 10*log(power)
-            var powerEpsilon = epsilon* epsilon / windowPower / sampleRate;
+            var powerEpsilon = epsilon * epsilon / windowPower / sampleRate;
             var decibelSpectrogram = MatrixTools.SpectrogramPower2DeciBels(powerSpectrogram, powerEpsilon, out var min, out var max);
             return decibelSpectrogram;
         }
@@ -137,6 +137,20 @@ public static FrequencyScale GetStandardOctaveScale(FrequencyScale scale)
             return octaveSpectrogram;
         }
 
+        /// <summary>
+        /// Converts Amplitude Spectrogram to Power Spectrogram.
+        /// Square the amplitude values to get power.
+        /// Power values must be adjusted for the power in the FFT window and also for the sample rate.
+        /// Must divide by the window power to remove its contribution to amplitude values.
+        /// Must divide by sample rate to get average power per signal sample.
+        /// NOTE: Multiply by 2 to accomodate two spectral components, ie positive and neg freq.
+        ///       BUT the last nyquist bin does not require a factor of two.
+        ///       However this method is called only by octave reduced frequency scales where the nyquist bin is just one of several.
+        /// </summary>
+        /// <param name="amplitudeM">The frequency scaled amplitude spectrogram.</param>
+        /// <param name="windowPower">Power of the FFT window.</param>
+        /// <param name="sampleRate">The sample rate of the original recording.</param>
+        /// <returns>The Power Spectrogram as matrix. Each spectrum is a matrix row.</returns>
         public static double[,] ConvertAmplitudeToPowerSpectrogram(double[,] amplitudeM, double windowPower, int sampleRate)
         {
             int frameCount = amplitudeM.GetLength(0);
@@ -147,7 +161,7 @@ public static FrequencyScale GetStandardOctaveScale(FrequencyScale scale)
 
             // Square the values to calculate power.
             // Must multiply by 2 to accomodate two spectral components, ie positive and neg freq.
-            for (int j = 0; j < binCount - 1; j++)
+            for (int j = 0; j < binCount; j++)
             {
                 //foreach time step or frame
                 for (int i = 0; i < frameCount; i++)
@@ -161,15 +175,7 @@ public static FrequencyScale GetStandardOctaveScale(FrequencyScale scale)
             //foreach time step or frame
             for (int i = 0; i < frameCount; i++)
             {
-                //calculate power of the DC value
-                if (amplitudeM[i, binCount - 1] < epsilon)
-                {
-                    powerSpectra[i, binCount - 1] = minPow;
-                }
-                else
-                {
-                    powerSpectra[i, binCount - 1] = amplitudeM[i, binCount - 1] * amplitudeM[i, binCount - 1] / windowPower / sampleRate;
-                }
+                powerSpectra[i, binCount - 1] = amplitudeM[i, binCount - 1] * amplitudeM[i, binCount - 1] / windowPower / sampleRate;
             }
             */
 
@@ -328,6 +334,7 @@ public static FrequencyScale GetDataReductionScale(FrequencyScale scale)
             int sr = 22050;
             int frameSize = 512;
             scale.Nyquist = sr / 2;
+            scale.WindowSize = frameSize;
 
             // linear reduction of the lower spectrum from 0 - 2 kHz.
             scale.LinearBound = 2000;

src/TowseyLibrary/MatrixTools.cs

@@ -519,7 +519,7 @@ public static double[] Matrix2Array(double[,] matrix)
         }
 
         /// <summary>
-        /// Convert the power values in a matrix of spectrogram values to Decibels.
+        /// Convert the power values in a matrix of spectrogram values to Decibels using: dB = 10*log10(power).
         /// Assume that all matrix values are positive i.e. due to prior noise removal.
         /// NOTE: This method also returns the min and max decibel values in the passed matrix.
         /// NOTE: A decibel value should be a ratio.
@@ -528,48 +528,52 @@ public static double[] Matrix2Array(double[,] matrix)
         /// <param name="m">matrix of positive power values.</param>
         /// <param name="min">min value to be return by out.</param>
         /// <param name="max">max value to be return by out.</param>
-        public static double[,] Power2DeciBels(double[,] m, double epsilon, out double min, out double max)
+        public static double[,] SpectrogramPower2DeciBels(double[,] m, double powerEpsilon, out double min, out double max)
         {
+            //convert epsilon power to decibels
+            double minDecibels = 10 * Math.Log10(powerEpsilon);
+
             min = double.MaxValue;
             max = -double.MaxValue;
 
             int rows = m.GetLength(0);
             int cols = m.GetLength(1);
-            double[,] ret = new double[rows, cols];
+            double[,] returnM = new double[rows, cols];
 
             for (int i = 0; i < rows; i++)
             {
                 for (int j = 0; j < cols; j++)
                 {
-                    double dBels = 10 * Math.Log10(epsilon); //convert power to decibels
-                    if (m[i, j] > epsilon)
+                    if (m[i, j] <= powerEpsilon)
                     {
-                        dBels = 10 * Math.Log10(m[i, j]); //convert power to decibels
+                        returnM[i, j] = minDecibels;
+                    }
+                    else
+                    {
+                        returnM[i, j] = 10 * Math.Log10(m[i, j]);
                     }
 
-                    if (dBels < min)
+                    if (returnM[i, j] < min)
                     {
-                        min = dBels;
+                        min = returnM[i, j];
                     }
                     else
-                    if (dBels > max)
+                    if (returnM[i, j] > max)
                     {
-                        max = dBels;
+                        max = returnM[i, j];
                     }
-
-                    ret[i, j] = dBels;
                 }
             }
 
-            return ret;
+            return returnM;
         }
 
         /// <summary>
-        /// Convert the power values in a matrix of spectrogram values to Decibels.
+        /// Convert the Decibels values in a matrix of spectrogram values to power values.
         /// Assume that all matrix values are positive due to prior noise removal.
         /// </summary>
         /// <param name="m">matrix of positive Decibel values.</param>
-        public static double[,] Decibels2Power(double[,] m)
+        public static double[,] SpectrogramDecibels2Power(double[,] m)
         {
             int rows = m.GetLength(0);
             int cols = m.GetLength(1);