|
| 1 | +//----------------------------------------------------------------------------- |
| 2 | +// Filename: AudioScope.cs |
| 3 | +// |
| 4 | +// Description: Implementation of a Hilbert filter to visualise audio input. |
| 5 | +// Originally based on https://github.com/conundrumer/visual-music-workshop. |
| 6 | + |
| 7 | +// Author(s): |
| 8 | +// Aaron Clauson (aaron@sipsorcery.com) |
| 9 | +// |
| 10 | +// History: |
| 11 | +// 29 Feb 2020 Aaron Clauson Created, Dublin, Ireland. |
| 12 | +// |
| 13 | +// License: |
| 14 | +// BSD 3-Clause "New" or "Revised" License, see included LICENSE.md file. |
| 15 | +//----------------------------------------------------------------------------- |
| 16 | + |
| 17 | +using System; |
| 18 | +using System.Linq; |
| 19 | +using System.Numerics; |
| 20 | +using MathNet.Numerics; |
| 21 | +using MathNet.Numerics.IntegralTransforms; |
| 22 | + |
| 23 | +namespace AudioScope |
| 24 | +{ |
| 25 | + public class LowPassFilter |
| 26 | + { |
| 27 | + private readonly float _k; |
| 28 | + private readonly float _norm; |
| 29 | + private readonly float _a0; |
| 30 | + private readonly float _a1; |
| 31 | + private readonly float _a2; |
| 32 | + private readonly float _b1; |
| 33 | + private readonly float _b2; |
| 34 | + |
| 35 | + private float _w1 = 0.0f; |
| 36 | + private float _w2 = 0.0f; |
| 37 | + |
| 38 | + public LowPassFilter(float n, float q) |
| 39 | + { |
| 40 | + _k = (float)Math.Tan((0.5 * n * Math.PI)); |
| 41 | + _norm = 1.0f / (1.0f + _k / q + _k * _k); |
| 42 | + _a0 = _k * _k * _norm; |
| 43 | + _a1 = 2.0f * _a0; |
| 44 | + _a2 = _a0; |
| 45 | + _b1 = 2.0f * (_k * _k - 1.0f) * _norm; |
| 46 | + _b2 = (1.0f - _k / q + _k * _k) * _norm; |
| 47 | + } |
| 48 | + |
| 49 | + public float Apply(float x) |
| 50 | + { |
| 51 | + float w0 = x - _b1 * _w1 - _b2 * _w2; |
| 52 | + float y = _a0 * w0 + _a1 * _w1 + _a2 * _w2; |
| 53 | + _w2 = _w1; |
| 54 | + _w1 = w0; |
| 55 | + |
| 56 | + return y; |
| 57 | + } |
| 58 | + } |
| 59 | + |
| 60 | + public class AudioScope |
| 61 | + { |
| 62 | + public const int NUM_CHANNELS = 1; |
| 63 | + public const int SAMPLE_RATE = 44100; |
| 64 | + public const float maxAmplitude = 4.0F; |
| 65 | + public const int B = (1 << 16) - 1; |
| 66 | + public const int M = 4; |
| 67 | + public const int FFT_SIZE = 1024; |
| 68 | + public const int MID = (FFT_SIZE - 1) / 2; |
| 69 | + public const float DELAY_TIME = MID / SAMPLE_RATE; |
| 70 | + public const float GAIN = 1.0f; |
| 71 | + public const int BUFFER_SIZE = 256; |
| 72 | + public const int CIRCULAR_BUFFER_SAMPLES = 3; |
| 73 | + public const float CUTOFF_FREQ = 0.5f; |
| 74 | + |
| 75 | + private const int DISPLAY_ARRAY_STRIDE = 4; // Each element sent to the display function needs to have 4 floats. |
| 76 | + private const int PREVIOUS_SAMPLES_LENGTH = 3 * DISPLAY_ARRAY_STRIDE; |
| 77 | + |
| 78 | + private Complex[] _analytic; |
| 79 | + private LowPassFilter _angleLowPass; |
| 80 | + private LowPassFilter _noiseLowPass; |
| 81 | + |
| 82 | + private Complex[] _timeRingBuffer = new Complex[2 * FFT_SIZE]; |
| 83 | + private int _timeIndex = 0; |
| 84 | + private float[] _previousResults = new float[3 * 4]; |
| 85 | + private Complex _prevInput = new Complex(0.0f, 0.0f); |
| 86 | + private Complex _prevDiff = new Complex(0.0f, 0.0f); |
| 87 | + private float[] _lastSample = []; |
| 88 | + |
| 89 | + public AudioScope() |
| 90 | + { |
| 91 | + uint n = FFT_SIZE; |
| 92 | + if (n % 2 == 0) |
| 93 | + { |
| 94 | + n -= 1; |
| 95 | + } |
| 96 | + |
| 97 | + _analytic = MakeAnalytic(n, FFT_SIZE); |
| 98 | + _angleLowPass = new LowPassFilter(0.01f, 0.5f); |
| 99 | + _noiseLowPass = new LowPassFilter(0.5f, 0.7f); |
| 100 | + } |
| 101 | + |
| 102 | + public float[] GetSample() |
| 103 | + { |
| 104 | + return _lastSample; |
| 105 | + } |
| 106 | + |
| 107 | + /// <summary> |
| 108 | + /// Called to process the audio input once the required number of samples are available. |
| 109 | + /// </summary> |
| 110 | + public void ProcessSample(Complex[] samples) |
| 111 | + { |
| 112 | + Array.Copy(samples, 0, _timeRingBuffer, _timeIndex, samples.Length > FFT_SIZE ? FFT_SIZE : samples.Length); |
| 113 | + Array.Copy(samples, 0, _timeRingBuffer, _timeIndex + FFT_SIZE, samples.Length > (_timeRingBuffer.Length/2 - _timeIndex) ? _timeRingBuffer.Length / 2 - _timeIndex : samples.Length); |
| 114 | + |
| 115 | + _timeIndex = (_timeIndex + samples.Length) % FFT_SIZE; |
| 116 | + |
| 117 | + var freqBuffer = _timeRingBuffer.Skip(_timeIndex).Take(FFT_SIZE).ToArray(); |
| 118 | + |
| 119 | + Fourier.Forward(freqBuffer, FourierOptions.NoScaling); |
| 120 | + |
| 121 | + for (int j = 0; j < freqBuffer.Length; j++) |
| 122 | + { |
| 123 | + freqBuffer[j] = freqBuffer[j] * _analytic[j]; |
| 124 | + } |
| 125 | + |
| 126 | + Fourier.Inverse(freqBuffer, FourierOptions.NoScaling); |
| 127 | + |
| 128 | + float scale = (float)FFT_SIZE; |
| 129 | + |
| 130 | + var complexAnalyticBuffer = freqBuffer.Skip(FFT_SIZE - BUFFER_SIZE).Take(BUFFER_SIZE).ToArray(); |
| 131 | + var data = new float[BUFFER_SIZE * DISPLAY_ARRAY_STRIDE + PREVIOUS_SAMPLES_LENGTH]; |
| 132 | + |
| 133 | + for (int k = 0; k < complexAnalyticBuffer.Length; k++) |
| 134 | + { |
| 135 | + var diff = complexAnalyticBuffer[k] - _prevInput; |
| 136 | + _prevInput = complexAnalyticBuffer[k]; |
| 137 | + |
| 138 | + var angle = (float)Math.Max(Math.Log(Math.Abs(GetAngle(diff, _prevDiff)), 2.0f), -1.0e12); |
| 139 | + _prevDiff = diff; |
| 140 | + var output = _angleLowPass.Apply(angle); |
| 141 | + |
| 142 | + data[k * DISPLAY_ARRAY_STRIDE] = (float)(complexAnalyticBuffer[k].Real / scale); |
| 143 | + data[k * DISPLAY_ARRAY_STRIDE + 1] = (float)(complexAnalyticBuffer[k].Imaginary / scale); |
| 144 | + data[k * DISPLAY_ARRAY_STRIDE + 2] = (float)Math.Pow(2, output); // Smoothed angular velocity. |
| 145 | + data[k * DISPLAY_ARRAY_STRIDE + 3] = _noiseLowPass.Apply((float)Math.Abs(angle - output)); // Average angular noise. |
| 146 | + } |
| 147 | + |
| 148 | + Array.Copy(_previousResults, 0, data, 0, PREVIOUS_SAMPLES_LENGTH); |
| 149 | + _lastSample = data; |
| 150 | + |
| 151 | + _previousResults = data.Skip(data.Length - PREVIOUS_SAMPLES_LENGTH).ToArray(); |
| 152 | + } |
| 153 | + |
| 154 | + public static float GetAngle(Complex v, Complex u) |
| 155 | + { |
| 156 | + var len_v_mul_u = v.Norm() * u; |
| 157 | + var len_u_mul_v = u.Norm() * v; |
| 158 | + var left = (len_v_mul_u - len_u_mul_v).Norm(); |
| 159 | + var right = (len_v_mul_u + len_u_mul_v).Norm(); |
| 160 | + |
| 161 | + return (float)(Math.Atan2(left, right) / Math.PI); |
| 162 | + } |
| 163 | + |
| 164 | + private static Complex[] MakeAnalytic(uint n, uint m) |
| 165 | + { |
| 166 | + var impulse = new Complex[m]; |
| 167 | + |
| 168 | + var mid = (n - 1) / 2; |
| 169 | + |
| 170 | + impulse[mid] = new Complex(1.0f, 0.0f); |
| 171 | + float re = -1.0f / (mid - 1); |
| 172 | + for (int i = 1; i < mid + 1; i++) |
| 173 | + { |
| 174 | + if (i % 2 == 0) |
| 175 | + { |
| 176 | + impulse[mid + i] = new Complex(re, impulse[mid + i].Imaginary); |
| 177 | + impulse[mid - i] = new Complex(re, impulse[mid - i].Imaginary); |
| 178 | + } |
| 179 | + else |
| 180 | + { |
| 181 | + float im = (float)(2.0 / Math.PI / i); |
| 182 | + impulse[mid + i] = new Complex(impulse[mid + i].Real, im); |
| 183 | + impulse[mid - i] = new Complex(impulse[mid - i].Real, -im); |
| 184 | + } |
| 185 | + // hamming window |
| 186 | + var k = 0.53836 + 0.46164 * Math.Cos(i * Math.PI / (mid + 1)); |
| 187 | + impulse[mid + i] = new Complex((float)(impulse[mid + i].Real * k), (float)(impulse[mid + i].Imaginary * k)); |
| 188 | + impulse[mid - i] = new Complex((float)(impulse[mid - i].Real * k), (float)(impulse[mid - i].Imaginary * k)); |
| 189 | + } |
| 190 | + |
| 191 | + Fourier.Forward(impulse, FourierOptions.NoScaling); |
| 192 | + |
| 193 | + return impulse; |
| 194 | + } |
| 195 | + } |
| 196 | +} |
0 commit comments