Skip to content

Commit fde1592

Browse files
dthalerdthaler
authored
Update FFT algorithm (#283)
* Switch FFT library to FftSharp Signed-off-by: Dave Thaler <[email protected]> * Save to a file for debugging Signed-off-by: Dave Thaler <[email protected]> Change display scale to decibels But keep magnitude internally Add Hann window to FFT algorithm Signed-off-by: Dave Thaler <[email protected]> Debugging Signed-off-by: Dave Thaler <[email protected]> Fix WAV header parsing Signed-off-by: Dave Thaler <[email protected]> * DEBUG: read from a file Signed-off-by: Dave Thaler <[email protected]> * Fix line fill for decibels Signed-off-by: Dave Thaler <[email protected]> * Update definitions of intelligibility threshold configuration Signed-off-by: Dave Thaler <[email protected]> * Sort labels to fix spectral density display Signed-off-by: Dave Thaler <[email protected]> * Update thresholds Signed-off-by: Dave Thaler <[email protected]> --------- Signed-off-by: Dave Thaler <[email protected]>
1 parent 8931407 commit fde1592

File tree

6 files changed

+109
-74
lines changed

6 files changed

+109
-74
lines changed

‎OrcanodeMonitor/Core/Fetcher.cs‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -961,7 +961,11 @@ private static void AddHydrophoneStreamStatusEvent(OrcanodeMonitorContext contex
961961
try
962962
{
963963
using Stream stream = await _httpClient.GetStreamAsync(uri);
964+
#if true
964965
FrequencyInfo frequencyInfo = await FfmpegCoreAnalyzer.AnalyzeAudioStreamAsync(stream, oldStatus);
966+
#else
967+
FrequencyInfo frequencyInfo = await FfmpegCoreAnalyzer.AnalyzeFileAsync("output-2channels.wav", oldStatus);
968+
#endif
965969
frequencyInfo.AudioSampleUrl = uri.AbsoluteUri;
966970
return frequencyInfo;
967971
}

‎OrcanodeMonitor/Core/FrequencyInfo.cs‎

Lines changed: 54 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
// Copyright (c) Orcanode Monitor contributors
22
// SPDX-License-Identifier: MIT
3-
using MathNet.Numerics.IntegralTransforms;
3+
using FftSharp;
44
using OrcanodeMonitor.Models;
55
using System.Diagnostics;
66
using System.Numerics;
@@ -46,33 +46,33 @@ private void ComputeFrequencyMagnitudes(float[] data, uint sampleRate, int chann
4646
}
4747

4848
int n = data.Length / channelCount;
49+
int nextPowerOfTwo = (int)Math.Pow(2, Math.Ceiling(Math.Log(n, 2)));
4950

50-
// Create an array of complex data for each channel.
51-
Complex[][] complexData = new Complex[channelCount][];
51+
// Create frequency magnitudes for each channel.
5252
for (int ch = 0; ch < channelCount; ch++)
5353
{
54-
complexData[ch] = new Complex[n];
55-
}
56-
57-
// Populate the complex arrays with channel data.
58-
for (int i = 0; i < n; i++)
59-
{
60-
for (int ch = 0; ch < channelCount; ch++)
54+
// Extract channel-specific data.
55+
double[] channelData = new double[nextPowerOfTwo];
56+
for (int i = 0; i < n; i++)
6157
{
62-
complexData[ch][i] = new Complex(data[i * channelCount + ch], 0);
58+
channelData[i] = data[i * channelCount + ch];
6359
}
64-
}
6560

66-
// Perform Fourier transform for each channel.
67-
for (int ch = 0; ch < channelCount; ch++)
68-
{
69-
Fourier.Forward(complexData[ch], FourierOptions.Matlab);
61+
// Apply Hann window.
62+
var hannWindow = new FftSharp.Windows.Hanning();
63+
hannWindow.ApplyInPlace(channelData);
64+
65+
// Perform FFT.
66+
Complex[] fftResult = FFT.Forward(channelData);
67+
double[] magnitudes = FFT.Magnitude(fftResult);
68+
69+
// Store frequency magnitudes for this channel.
7070
FrequencyMagnitudesForChannel[ch] = new Dictionary<double, double>(n / 2);
71-
for (int i = 0; i < n / 2; i++)
71+
for (int i = 0; i < magnitudes.Length / 2; i++) // Use only the first half (positive frequencies).
7272
{
73-
double magnitude = complexData[ch][i].Magnitude;
74-
double frequency = (((double)i) * sampleRate) / n;
75-
FrequencyMagnitudesForChannel[ch][frequency] = magnitude;
73+
double frequency = (((double)i) * sampleRate) / nextPowerOfTwo;
74+
FrequencyMagnitudesForChannel[ch][frequency] = magnitudes[i];
75+
double decibels = magnitudes[i] > 0 ? 20 * Math.Log10(magnitudes[i]) : double.NegativeInfinity; // DEBUG
7676
}
7777
}
7878

@@ -99,30 +99,45 @@ private void ComputeFrequencyMagnitudes(float[] data, uint sampleRate, int chann
9999
}
100100
}
101101

102-
// We consider anything above this average magnitude as not silence.
103-
const double _defaultMaxSilenceMagnitude = 20.0;
104-
public static double MaxSilenceMagnitude
102+
private static double DecibelsToMagnitude(double decibels)
103+
{
104+
double magnitude = Math.Pow(10, decibels / 20);
105+
return magnitude;
106+
}
107+
108+
private static double MagnitudeToDecibels(double magnitude)
109+
{
110+
double dB = 20 * Math.Log10(magnitude);
111+
return dB;
112+
}
113+
114+
// We consider anything above this average decibels as not silence.
115+
const double _defaultMaxSilenceDecibels = -80;
116+
public static double MaxSilenceDecibels
105117
{
106118
get
107119
{
108-
string? maxSilenceMagnitudeString = Environment.GetEnvironmentVariable("ORCASOUND_MAX_SILENCE_MAGNITUDE");
109-
double maxSilenceMagnitude = double.TryParse(maxSilenceMagnitudeString, out var magnitude) ? magnitude : _defaultMaxSilenceMagnitude;
110-
return maxSilenceMagnitude;
120+
string? maxSilenceDecibelsString = Environment.GetEnvironmentVariable("ORCASOUND_MAX_SILENCE_DECIBELS");
121+
double maxSilenceDecibels = double.TryParse(maxSilenceDecibelsString, out var decibels) ? decibels : _defaultMaxSilenceDecibels;
122+
return maxSilenceDecibels;
111123
}
112124
}
113125

114-
// We consider anything below this average magnitude as silence.
115-
// The lowest normal value we have seen is 7.7.
116-
const double _defaultMinNoiseMagnitude = 7.0;
117-
public static double MinNoiseMagnitude
126+
public static double MaxSilenceMagnitude => DecibelsToMagnitude(MaxSilenceDecibels);
127+
128+
// We consider anything below this average decibels as silence.
129+
// The lowest normal value we have seen is -98.
130+
const double _defaultMinNoiseDecibels = -90;
131+
public static double MinNoiseDecibels
118132
{
119133
get
120134
{
121-
string? minNoiseMagnitudeString = Environment.GetEnvironmentVariable("ORCASOUND_MIN_NOISE_MAGNITUDE");
122-
double minNoiseMagnitude = double.TryParse(minNoiseMagnitudeString, out var magnitude) ? magnitude : _defaultMinNoiseMagnitude;
123-
return minNoiseMagnitude;
135+
string? minNoiseDecibelsString = Environment.GetEnvironmentVariable("ORCASOUND_MIN_NOISE_DECIBELS");
136+
double minNoiseDecibels = double.TryParse(minNoiseDecibelsString, out var decibels) ? decibels : _defaultMinNoiseDecibels;
137+
return minNoiseDecibels;
124138
}
125139
}
140+
public static double MinNoiseMagnitude => DecibelsToMagnitude(MinNoiseDecibels);
126141

127142
// Minimum ratio of magnitude outside the hum range to magnitude
128143
// within the hum range. So far the max in a known-unintelligible
@@ -280,22 +295,24 @@ public double GetTotalHumMagnitude(Dictionary<double, double> frequencyMagnitude
280295

281296
private OrcanodeOnlineStatus GetStatus(OrcanodeOnlineStatus oldStatus, int? channel = null)
282297
{
283-
double max = GetMaxMagnitude(channel);
284-
if (max < MinNoiseMagnitude)
298+
double maxMagnitude = GetMaxMagnitude(channel);
299+
double maxDecibels = MagnitudeToDecibels(maxMagnitude);
300+
if (maxDecibels < MinNoiseDecibels)
285301
{
286302
// File contains mostly silence across all frequencies.
287303
return OrcanodeOnlineStatus.Silent;
288304
}
289305

290-
if ((max <= MaxSilenceMagnitude) && (oldStatus == OrcanodeOnlineStatus.Silent))
306+
if (maxDecibels <= MaxSilenceDecibels)
291307
{
292308
// In between the min and max silence range, so keep previous status.
293309
return oldStatus;
294310
}
295311

296312
// Find the total magnitude outside the audio hum range.
297313
double maxNonHumMagnitude = GetMaxNonHumMagnitude(channel);
298-
if (maxNonHumMagnitude < MinNoiseMagnitude)
314+
double maxNonHumDecibels = MagnitudeToDecibels(maxNonHumMagnitude);
315+
if (maxNonHumDecibels < MinNoiseDecibels)
299316
{
300317
// Just silence outside the hum range, no signal.
301318
return OrcanodeOnlineStatus.Unintelligible;

‎OrcanodeMonitor/OrcanodeMonitor.csproj‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -13,7 +13,7 @@
1313
<ItemGroup>
1414
<PackageReference Include="FFMpegCore" Version="5.1.0" />
1515
<PackageReference Include="FFMpegInstaller.Windows.x64" Version="0.1.10" />
16-
<PackageReference Include="MathNet.Numerics" Version="5.0.0" />
16+
<PackageReference Include="FftSharp" Version="2.2.0" />
1717
<PackageReference Include="Microsoft.AspNetCore.Diagnostics.EntityFrameworkCore" Version="8.0.10" />
1818
<PackageReference Include="Microsoft.EntityFrameworkCore.Cosmos" Version="8.0.10" />
1919
<PackageReference Include="Microsoft.EntityFrameworkCore.Proxies" Version="8.0.10" />

‎OrcanodeMonitor/Pages/SpectralDensity.cshtml‎

Lines changed: 13 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -21,8 +21,8 @@
2121
// Client side function to initialize all charts.
2222
function initializeCharts2(labels, jsonChannelDatasets)
2323
{
24-
var maxSilenceMagnitude = @Model.MaxSilenceMagnitude;
25-
var minNoiseMagnitude = @Model.MinNoiseMagnitude;
24+
var maxSilenceDecibels = @Model.MaxSilenceDecibels;
25+
var minNoiseDecibels = @Model.MinNoiseDecibels;
2626
2727
// Channel Count = @Model.ChannelCount
2828
@@ -36,8 +36,8 @@
3636
@: '@id',
3737
@: labels,
3838
@: jsonChannelDatasets[@i],
39-
@: maxSilenceMagnitude,
40-
@: minNoiseMagnitude,
39+
@: maxSilenceDecibels,
40+
@: minNoiseDecibels,
4141
@: @i,
4242
@: '@backgroundColor',
4343
@: '@borderColor'
@@ -53,7 +53,7 @@
5353
}
5454
5555
// Client side function to initialize a chart.
56-
function initializeChart(elementId, labels, jsonChannelDataset, maxSilenceMagnitude, minNoiseMagnitude, i, backgroundColor, borderColor)
56+
function initializeChart(elementId, labels, jsonChannelDataset, maxSilenceDecibels, minNoiseDecibels, i, backgroundColor, borderColor)
5757
{
5858
var ctx = document.getElementById(elementId).getContext('2d');
5959
var myLineChart = new Chart(ctx, {
@@ -67,24 +67,24 @@
6767
backgroundColor: backgroundColor,
6868
borderColor: borderColor,
6969
borderWidth: 1,
70-
fill: true
70+
fill: "start"
7171
},
7272
{
73-
label: 'Max Noise Magnitude',
73+
label: 'Max Silence Decibels',
7474
data: [
75-
{ x: labels[0], y: maxSilenceMagnitude },
76-
{ x: labels[labels.length - 1], y: maxSilenceMagnitude }
75+
{ x: labels[0], y: maxSilenceDecibels },
76+
{ x: labels[labels.length - 1], y: maxSilenceDecibels }
7777
],
7878
borderColor: 'rgba(54, 162, 235, 1)',
7979
borderWidth: 2,
8080
fill: false,
8181
pointRadius: 0 // Hide points on this line
8282
},
8383
{
84-
label: 'Min Signal Magnitude',
84+
label: 'Min Noise Decibels',
8585
data: [
86-
{ x: labels[0], y: minNoiseMagnitude },
87-
{ x: labels[labels.length - 1], y: minNoiseMagnitude }
86+
{ x: labels[0], y: minNoiseDecibels },
87+
{ x: labels[labels.length - 1], y: minNoiseDecibels }
8888
],
8989
borderColor: 'rgba(255, 99, 132, 1)',
9090
borderWidth: 2,
@@ -102,10 +102,9 @@
102102
}
103103
},
104104
y: {
105-
beginAtZero: true,
106105
title: {
107106
display: true,
108-
text: 'Magnitude'
107+
text: 'dB'
109108
}
110109
}
111110
}

‎OrcanodeMonitor/Pages/SpectralDensity.cshtml.cs‎

Lines changed: 35 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -36,8 +36,13 @@ public class SpectralDensityModel : PageModel
3636
public int MaxNonHumMagnitude { get; private set; }
3737
public int SignalRatio { get; private set; }
3838
public string Status { get; private set; }
39-
public double MaxSilenceMagnitude => FrequencyInfo.MaxSilenceMagnitude;
40-
public double MinNoiseMagnitude => FrequencyInfo.MinNoiseMagnitude;
39+
private static double MagnitudeToDecibels(double magnitude)
40+
{
41+
double dB = 20 * Math.Log10(magnitude);
42+
return dB;
43+
}
44+
public double MaxSilenceDecibels => FrequencyInfo.MaxSilenceDecibels;
45+
public double MinNoiseDecibels => FrequencyInfo.MinNoiseDecibels;
4146
public string LastModified { get; private set; }
4247

4348
public SpectralDensityModel(OrcanodeMonitorContext context, ILogger<SpectralDensityModel> logger)
@@ -63,7 +68,7 @@ public SpectralDensityModel(OrcanodeMonitorContext context, ILogger<SpectralDens
6368
/// </summary>
6469
private const int POINT_COUNT = 1000;
6570

66-
private void FillInGraphPoints(List<string> labels, List<double> maxBucketMagnitudeList, int? channel = null)
71+
private void FillInGraphPoints(List<string> labels, List<double> maxBucketDecibelsList, int? channel = null)
6772
{
6873
if (_frequencyInfo == null)
6974
{
@@ -74,37 +79,41 @@ private void FillInGraphPoints(List<string> labels, List<double> maxBucketMagnit
7479
double b = Math.Pow(MAX_FREQUENCY, 1.0 / POINT_COUNT);
7580
double logb = Math.Log(b);
7681

77-
var maxBucketMagnitude = new double[POINT_COUNT];
82+
var maxBucketDecibels = new double[POINT_COUNT];
83+
for (int i = 0; i < POINT_COUNT; i++) {
84+
maxBucketDecibels[i] = double.NegativeInfinity;
85+
}
7886
var maxBucketFrequency = new int[POINT_COUNT];
7987
foreach (var pair in _frequencyInfo.GetFrequencyMagnitudes(channel))
8088
{
8189
double frequency = pair.Key;
8290
double magnitude = pair.Value;
91+
double dB = MagnitudeToDecibels(magnitude);
8392
int bucket = (frequency < 1) ? 0 : Math.Min(POINT_COUNT - 1, (int)(Math.Log(frequency) / logb));
84-
if (maxBucketMagnitude[bucket] < magnitude)
93+
if (maxBucketDecibels[bucket] < dB)
8594
{
86-
maxBucketMagnitude[bucket] = magnitude;
95+
maxBucketDecibels[bucket] = dB;
8796
maxBucketFrequency[bucket] = (int)Math.Round(frequency);
8897
}
8998
}
9099
for (int i = 0; i < POINT_COUNT; i++)
91100
{
92-
if (maxBucketMagnitude[i] > 0)
101+
if (maxBucketDecibels[i] > double.NegativeInfinity)
93102
{
94103
labels.Add(maxBucketFrequency[i].ToString());
95-
maxBucketMagnitudeList.Add(maxBucketMagnitude[i]);
104+
maxBucketDecibelsList.Add(maxBucketDecibels[i]);
96105
}
97106
}
98107
}
99108

100-
private double GetBucketMagnitude(string label, List<string> labels, List<double> magnitudes)
109+
private double GetBucketDecibels(string label, List<string> labels, List<double> decibels)
101110
{
102-
double max = 0;
111+
double max = double.NegativeInfinity;
103112
for (int i = 0; i < labels.Count; i++)
104113
{
105-
if (labels[i] == label && magnitudes[i] > max)
114+
if (labels[i] == label && decibels[i] > max)
106115
{
107-
max = magnitudes[i];
116+
max = decibels[i];
108117
}
109118
}
110119
return max;
@@ -119,15 +128,15 @@ private void UpdateFrequencyInfo()
119128

120129
// Compute graph points.
121130
var summaryLabels = new List<string>();
122-
var summaryMaxBucketMagnitude = new List<double>();
123-
FillInGraphPoints(summaryLabels, summaryMaxBucketMagnitude);
131+
var summaryMaxBucketDecibels = new List<double>();
132+
FillInGraphPoints(summaryLabels, summaryMaxBucketDecibels);
124133
var channelLabels = new List<string>[_frequencyInfo.ChannelCount];
125-
var channelMaxBucketMagnitude = new List<double>[_frequencyInfo.ChannelCount];
134+
var channelMaxBucketDecibels = new List<double>[_frequencyInfo.ChannelCount];
126135
for (int i = 0; i < _frequencyInfo.ChannelCount; i++)
127136
{
128137
channelLabels[i] = new List<string>();
129-
channelMaxBucketMagnitude[i] = new List<double>();
130-
FillInGraphPoints(channelLabels[i], channelMaxBucketMagnitude[i], i);
138+
channelMaxBucketDecibels[i] = new List<double>();
139+
FillInGraphPoints(channelLabels[i], channelMaxBucketDecibels[i], i);
131140
}
132141

133142
// Collect all labels.
@@ -136,21 +145,27 @@ private void UpdateFrequencyInfo()
136145
{
137146
mainLabels.UnionWith(channelLabels[i]);
138147
}
139-
_labels = mainLabels.ToList();
148+
149+
// Sort labels numerically.
150+
_labels = mainLabels
151+
.Select(label => int.Parse(label)) // Convert to integers for sorting.
152+
.OrderBy(num => num) // Sort in ascending order.
153+
.Select(num => num.ToString()) // Convert back to strings.
154+
.ToList();
140155

141156
// Align data.
142157
var summaryDataset = _labels.Select(label => new
143158
{
144159
x = label,
145-
y = summaryLabels.Contains(label) ? GetBucketMagnitude(label, summaryLabels, summaryMaxBucketMagnitude) : (double?)null
160+
y = summaryLabels.Contains(label) ? GetBucketDecibels(label, summaryLabels, summaryMaxBucketDecibels) : (double?)null
146161
}).ToList<object>();
147162
var channelDatasets = new List<List<object>>();
148163
for (int i = 0; i < _frequencyInfo.ChannelCount; i++)
149164
{
150165
var channelDataset = _labels.Select(label => new
151166
{
152167
x = label,
153-
y = channelLabels[i].Contains(label) ? GetBucketMagnitude(label, channelLabels[i], channelMaxBucketMagnitude[i]) : (double?)null
168+
y = channelLabels[i].Contains(label) ? GetBucketDecibels(label, channelLabels[i], channelMaxBucketDecibels[i]) : (double?)null
154169
}).ToList<object>();
155170
channelDatasets.Add(channelDataset);
156171
}

‎docs/Design.md‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -115,9 +115,9 @@ The following state will be stored per orcanode:
115115

116116
**ORCASOUND_MIN_INTELLIGIBLE_SIGNAL_PERCENT**: The minimum percentage of total magnitude across all frequencies outside the hum range vs magnitude in hum range (multiples of 60 Hz), needed to determine that an audio stream is intelligible. Default: 1400
117117

118-
**ORCASOUND_MAX_SILENCE_MAGNITUDE**: The maximum magnitude at which a stream stream might still be considered unintelligible due to silence. Default: 20
118+
**ORCASOUND_MAX_SILENCE_DECIBELS**: The maximum decibel level at which a stream stream might still be considered unintelligible due to silence. Default: -80
119119

120-
**ORCASOUND_MIN_NOISE_MAGNITUDE**: The minimum magnitude at which a stream stream might still be considered intelligible. Default: 7
120+
**ORCASOUND_MIN_NOISE_DECIBELS**: The minimum decibel level at which a stream stream might still be considered intelligible. Default: -90
121121

122122
These magnitude thresholds work together to implement hysteresis in the noise detection:
123123
- Magnitudes below ORCASOUND_MIN_NOISE_MAGNITUDE are always considered silent.

0 commit comments

Comments
 (0)