|
1 | 1 | using System.Data; |
2 | | -using System.Diagnostics; |
3 | 2 |
|
4 | 3 | namespace NSD |
5 | 4 | { |
@@ -84,52 +83,93 @@ public static Spectrum StackedLinear(Memory<double> input, double sampleRate, in |
84 | 83 | return new Spectrum() { Frequencies = outputFrequencies.ToArray(), Values = outputValues.ToArray(), Averages = averages, Stacking = widths.Count }; |
85 | 84 | } |
86 | 85 |
|
87 | | - public static Spectrum Log(Memory<double> input, double sampleRateHz, double freqMin, double freqMax, int pointsPerDecade, int minimumAverages) |
| 86 | + private record WelchGoertzelJob(double Frequency, int SpectrumLength, int CalculatedAverages); |
| 87 | + public static Spectrum Log(Memory<double> input, double sampleRateHz, double freqMin, double freqMax, int pointsPerDecade, int minimumAverages, int minimumLength) |
88 | 88 | { |
89 | 89 | if (freqMax <= freqMin) |
90 | 90 | throw new ArgumentException("freqMax must be greater than freqMin"); |
91 | | - if (pointsPerDecade <= 0 || minimumAverages <= 0) |
92 | | - throw new ArgumentException("pointsPerDecade, and minimumAverages must be positive"); |
| 91 | + if (pointsPerDecade <= 0 || minimumAverages <= 0 || minimumLength <= 0) |
| 92 | + throw new ArgumentException("pointsPerDecade, minimumAverages, and minimumLength must be positive"); |
93 | 93 | if (sampleRateHz <= 0) |
94 | 94 | throw new ArgumentException("sampleRateHz must be positive"); |
95 | | - if (freqMin < sampleRateHz / input.Length) |
96 | | - freqMin = sampleRateHz / input.Length; |
97 | | - if (freqMax > sampleRateHz / 2) |
98 | | - freqMax = sampleRateHz / 2; |
99 | 95 |
|
100 | 96 | Windows.FTNI(1, out double optimumOverlap, out double NENBW); |
101 | 97 | int firstUsableBinForWindow = (int)Math.Ceiling(NENBW); |
102 | | - double decades = Math.Log10(freqMax / freqMin); |
| 98 | + |
| 99 | + // To do: |
| 100 | + // For the purposes of the frequencies calculation, round freqMax/freqMin to nearest major decade line. |
| 101 | + // This ensures consistency of X-coordinate over various view widths. |
| 102 | + double decadeMin = RoundToDecade(freqMin, RoundingMode.Down); |
| 103 | + double decadeMax = RoundToDecade(freqMax, RoundingMode.Up); |
| 104 | + double decades = Math.Log10(decadeMax / decadeMin); |
103 | 105 | int desiredNumberOfPoints = (int)(decades * pointsPerDecade) + 1; // + 1 to get points on the decade grid lines |
104 | 106 |
|
105 | | - double g = Math.Log(freqMax) - Math.Log(freqMin); |
106 | | - double[] frequencies = Enumerable.Range(0, desiredNumberOfPoints).Select(j => freqMin * Math.Exp(j * g / (desiredNumberOfPoints - 1))).ToArray(); |
| 107 | + double g = Math.Log(decadeMax) - Math.Log(decadeMin); |
| 108 | + double[] frequencies = Enumerable.Range(0, desiredNumberOfPoints).Select(j => decadeMin * Math.Exp(j * g / (desiredNumberOfPoints - 1))).ToArray(); |
107 | 109 | double[] spectrumResolution = frequencies.Select(freq => freq / firstUsableBinForWindow).ToArray(); |
108 | | - // spectrumResolution contains the 'desired resolutions' for each frequency bin, given the rule that we want the first usuable bin in the flat-top'd data. |
| 110 | + // spectrumResolution contains the 'desired resolutions' for each frequency bin, respecting the rule that we want the first usuable bin for the given window. |
| 111 | + |
| 112 | + int[] spectrumLengths = spectrumResolution.Select(resolution => (int)Math.Round(sampleRateHz / resolution)).ToArray(); |
109 | 113 |
|
110 | | - int[] spectrumLength = spectrumResolution.Select(val => (int)Math.Round(sampleRateHz / val)).ToArray(); // Segment lengths |
111 | | - //double[] actualSpectrumResolution = spectrumLength.Select(val => sampleRateHz / val).ToArray(); // Actual resolution |
112 | | - //double[] binNumber = frequencies.Select((val, index) => val / actualSpectrumResolution[index]).ToArray(); // Fourier tranform bin number (maybe validate that it doesn't deviate beyond +/-10%?) |
113 | | - int[] estimatedAverages = spectrumLength.Select(val => (int)((input.Length - val) / (val * (1.0 - optimumOverlap)))).ToArray(); |
| 114 | + // Create a job list of valid points to calculate |
| 115 | + double nyquistMax = sampleRateHz / 2; |
| 116 | + List<WelchGoertzelJob> jobs = []; |
| 117 | + for (int i = 0; i < frequencies.Length; i++) |
| 118 | + { |
| 119 | + if (frequencies[i] > nyquistMax) |
| 120 | + continue; |
| 121 | + if (TryCalculateAverages(input.Length, spectrumLengths[i], optimumOverlap, out var averages)) |
| 122 | + { |
| 123 | + if (averages >= minimumAverages) |
| 124 | + { |
| 125 | + // Increase spectrum length until minimumLength is met, or averages drops below minimumAverages. |
| 126 | + // This increases the spectral resolution at the top end of the chart, allowing 50Hz spikes (& similar) to be more visible |
| 127 | + var spectrumLength = spectrumLengths[i]; |
| 128 | + bool continueLoop = true; |
| 129 | + while (continueLoop) |
| 130 | + { |
| 131 | + if (spectrumLength < minimumLength && averages > minimumAverages) |
| 132 | + { |
| 133 | + var success = TryCalculateAverages(input.Length, spectrumLength * 2, optimumOverlap, out var newAverages); |
| 134 | + if (!success) |
| 135 | + break; |
| 136 | + if (averages > minimumAverages) |
| 137 | + { |
| 138 | + spectrumLength *= 2; |
| 139 | + averages = newAverages; |
| 140 | + continueLoop = true; |
| 141 | + } |
| 142 | + else |
| 143 | + { |
| 144 | + continueLoop = false; |
| 145 | + break; |
| 146 | + } |
| 147 | + |
| 148 | + } |
| 149 | + else |
| 150 | + { |
| 151 | + continueLoop = false; |
| 152 | + } |
| 153 | + } |
| 154 | + jobs.Add(new WelchGoertzelJob(frequencies[i], spectrumLength, averages)); |
| 155 | + } |
| 156 | + } |
| 157 | + } |
114 | 158 |
|
115 | 159 | var spectrum = new Dictionary<double, double>(); |
116 | | - var indices = Enumerable.Range(0, desiredNumberOfPoints).ToArray(); |
117 | | - for(int i = 0; i < frequencies.Length; i++) |
| 160 | + for (int i = 0; i < jobs.Count; i++) |
118 | 161 | { |
119 | | - spectrum[frequencies[i]] = double.NaN; |
| 162 | + spectrum[jobs[i].Frequency] = double.NaN; |
120 | 163 | } |
121 | 164 | object averageLock = new(); |
122 | 165 | int cumulativeAverage = 0; |
123 | | - //for (int i = 0; i < desiredNumberOfPoints; i++) |
124 | | - //foreach(var i in indices) |
125 | | - Parallel.ForEach(indices, new ParallelOptions { MaxDegreeOfParallelism = 8 }, i => |
| 166 | + //foreach(var job in jobs) |
| 167 | + Parallel.ForEach(jobs, new ParallelOptions { MaxDegreeOfParallelism = 8 }, job => |
126 | 168 | { |
127 | | - if (estimatedAverages[i] < minimumAverages) |
128 | | - return; |
129 | | - var result = RunWelchGoertzel(input, spectrumLength[i], frequencies[i], sampleRateHz, out var actualAverages); |
130 | | - if (estimatedAverages[i] != actualAverages) |
131 | | - Debug.WriteLine($"{estimatedAverages[i]} {actualAverages}"); |
132 | | - spectrum[frequencies[i]] = result; |
| 169 | + var result = RunWelchGoertzel(input, job.SpectrumLength, job.Frequency, sampleRateHz, out var actualAverages); |
| 170 | + if (job.CalculatedAverages != actualAverages) |
| 171 | + throw new Exception("Actual averages does not match calculated averages"); |
| 172 | + spectrum[job.Frequency] = result; |
133 | 173 | lock (averageLock) |
134 | 174 | { |
135 | 175 | cumulativeAverage += actualAverages; |
@@ -285,6 +325,39 @@ private static double S2(ReadOnlySpan<double> window) |
285 | 325 | } |
286 | 326 | return sumSquared; |
287 | 327 | } |
| 328 | + |
| 329 | + enum RoundingMode { Nearest, Up, Down } |
| 330 | + private static double RoundToDecade(double value, RoundingMode mode) |
| 331 | + { |
| 332 | + if (value <= 0) |
| 333 | + throw new ArgumentOutOfRangeException(nameof(value), "Value must be positive."); |
| 334 | + |
| 335 | + double log10 = Math.Log10(value); |
| 336 | + double exponent = mode switch |
| 337 | + { |
| 338 | + RoundingMode.Nearest => Math.Round(log10), |
| 339 | + RoundingMode.Up => Math.Ceiling(log10), |
| 340 | + RoundingMode.Down => Math.Floor(log10), |
| 341 | + _ => throw new ArgumentOutOfRangeException(nameof(mode), "Invalid rounding mode.") |
| 342 | + }; |
| 343 | + |
| 344 | + return Math.Pow(10, exponent); |
| 345 | + } |
| 346 | + |
| 347 | + private static bool TryCalculateAverages(int dataLength, int spectrumLength, double optimumOverlap, out int averages) |
| 348 | + { |
| 349 | + averages = 0; |
| 350 | + int overlap = (int)(spectrumLength * (1.0 - optimumOverlap)); |
| 351 | + if (overlap < 1) |
| 352 | + return false; |
| 353 | + int endIndex = spectrumLength; |
| 354 | + while (endIndex < dataLength) |
| 355 | + { |
| 356 | + averages++; |
| 357 | + endIndex += overlap; |
| 358 | + } |
| 359 | + return true; |
| 360 | + } |
288 | 361 | } |
289 | 362 | } |
290 | 363 |
|
0 commit comments