|
| 1 | +/** |
| 2 | + * Copyright © 2018 8eo Inc. |
| 3 | + */ |
| 4 | +package za.co.monadic.scopus.dsp |
| 5 | + |
| 6 | +/** |
| 7 | + * Structure needed to define a filter |
| 8 | + * @param order The order of the filter |
| 9 | + * @param a The coefficients of the polynomial defining the poles in the Z-domain |
| 10 | + * @param b The coefficients of the polynomial defining the xeros in the Z-domain |
| 11 | + */ |
| 12 | +case class Filter(order: Int, a: Array[Float], b: Array[Float]) { |
| 13 | + require(order == a.length -1 && order == b.length - 1, "Order and coefficient array sizes must be equal") |
| 14 | +} |
| 15 | + |
| 16 | +/** |
| 17 | + * Supplies Elliptical filters that are suitable as interpolation and decimation filters |
| 18 | + * These were generated using scipy's filter design tools. A filter for a decimation/interpolation |
| 19 | + * factor of 4 has a cut off frequency of 90% of the Nyquist frequency of the original sample rate |
| 20 | + * divided by 4. |
| 21 | + */ |
| 22 | +object MultirateFilterFactory { |
| 23 | + |
| 24 | + /** |
| 25 | + * Select a filter for a cut-off frequency for the envisaged decimation/interpolation rate |
| 26 | + * @param factor Decimation/interpolation |
| 27 | + * @return A Filter |
| 28 | + */ |
| 29 | + def apply(factor: Int): Filter = factor match { |
| 30 | + case 2 ⇒ |
| 31 | + Filter( |
| 32 | + 6, |
| 33 | + Array[Float](1.0000000000e+00f, -2.2150939834e+00f, 3.6340884990e+00f, -3.6053900178e+00f, 2.5896791922e+00f, |
| 34 | + -1.1678192254e+00f, 2.9534451224e-01f), |
| 35 | + Array[Float](2.9743123894e-02f, 5.3683139827e-02f, 9.9997243561e-02f, 1.0623698403e-01f, 9.9997243561e-02f, |
| 36 | + 5.3683139827e-02f, 2.9743123894e-02f) |
| 37 | + ) |
| 38 | + case 3 ⇒ |
| 39 | + Filter( |
| 40 | + 6, |
| 41 | + Array[Float](1.0000000000e+00f, -3.9597283064e+00f, 7.5464419229e+00f, -8.4777972935e+00f, 5.8803344415e+00f, |
| 42 | + -2.3780608271e+00f, 4.4209252816e-01f), |
| 43 | + Array[Float](1.1659133162e-02f, -7.6466313128e-03f, 2.3681010412e-02f, -7.8989770916e-03f, 2.3681010412e-02f, |
| 44 | + -7.6466313128e-03f, 1.1659133162e-02f) |
| 45 | + ) |
| 46 | + case 4 ⇒ |
| 47 | + Filter( |
| 48 | + 6, |
| 49 | + Array[Float](1.0000000000e+00f, -4.6768497674e+00f, 9.7621423046e+00f, -1.1501934308e+01f, 8.0339245654e+00f, |
| 50 | + -3.1489284186e+00f, 5.4189391928e-01f), |
| 51 | + Array[Float](7.2935097847e-03f, -1.5409691786e-02f, 2.5558502830e-02f, -2.5750839286e-02f, 2.5558502830e-02f, |
| 52 | + -1.5409691786e-02f, 7.2935097847e-03f) |
| 53 | + ) |
| 54 | + case 5 ⇒ |
| 55 | + Filter( |
| 56 | + 6, |
| 57 | + Array[Float](1.0000000000e+00f, -5.0470082922e+00f, 1.1055203383e+01f, -1.3388074677e+01f, 9.4329631691e+00f, |
| 58 | + -3.6626915823e+00f, 6.1243214263e-01f), |
| 59 | + Array[Float](5.5699828667e-03f, -1.7047032929e-02f, 3.0001912362e-02f, -3.4532704919e-02f, 3.0001912362e-02f, |
| 60 | + -1.7047032929e-02f, 5.5699828667e-03f) |
| 61 | + ) |
| 62 | + case 6 ⇒ |
| 63 | + Filter( |
| 64 | + 6, |
| 65 | + Array[Float](1.0000000000e+00f, -5.2664800962e+00f, 1.1874502680e+01f, -1.4637401088e+01f, 1.0390542726e+01f, |
| 66 | + -4.0247175699e+00f, 6.6453255902e-01f), |
| 67 | + Array[Float](4.7105692599e-03f, -1.7496959962e-02f, 3.3343643097e-02f, -4.0241782497e-02f, 3.3343643097e-02f, |
| 68 | + -1.7496959962e-02f, 4.7105692599e-03f) |
| 69 | + ) |
| 70 | + case _ ⇒ throw new RuntimeException("Unsupported decimation factor") |
| 71 | + } |
| 72 | +} |
| 73 | + |
| 74 | +/** |
| 75 | + * Retain every n'th sample in the sequence. The input array length must be an |
| 76 | + * integer multiple of the decimation factor else the decimate method will throw |
| 77 | + * an exception. |
| 78 | + */ |
| 79 | +trait Decimator { |
| 80 | + val factor: Int |
| 81 | + def decimate(x: Array[Float]): Array[Float] = { |
| 82 | + require(x.length % factor == 0, "Input array length must be a multiple of the decimation rate") |
| 83 | + var n = 0 |
| 84 | + var m = 0 |
| 85 | + val y = new Array[Float](x.length / factor) |
| 86 | + while (m < x.length) { |
| 87 | + y(n) = x(m) |
| 88 | + n += 1 |
| 89 | + m += factor |
| 90 | + } |
| 91 | + y |
| 92 | + } |
| 93 | +} |
| 94 | + |
| 95 | +/** |
| 96 | + * Inserts N-1 zeros between samples provided, taking care to account for array boundaries. If the |
| 97 | + * interpolation factor is 3, then the sequence [1,2,3] is mapped to [1,0,0,2,0,0,3,0,0] |
| 98 | + */ |
| 99 | +trait Interpolator { |
| 100 | + val factor: Int |
| 101 | + def interpolate(x: Array[Float]): Array[Float] = { |
| 102 | + val l = x.length |
| 103 | + val y = new Array[Float](l * factor) |
| 104 | + var n = 0 |
| 105 | + while (n < l) { |
| 106 | + y(factor * n) = x(n) |
| 107 | + n += 1 |
| 108 | + } |
| 109 | + y |
| 110 | + } |
| 111 | +} |
| 112 | + |
| 113 | +/** |
| 114 | + * Perform an IIR filter operation using the filter configuration provided in the constructor. |
| 115 | + * This implementation is unoptimised. |
| 116 | + * @param f The input data signal to filter |
| 117 | + */ |
| 118 | +class FilterIIR(f: Filter) { |
| 119 | + |
| 120 | + val state = new Array[Float](f.order + 1) |
| 121 | + |
| 122 | + @inline |
| 123 | + def filterOne(x: Float): Float = { |
| 124 | + var sumA = x |
| 125 | + var sumB = 0.0f |
| 126 | + var i = f.order |
| 127 | + while (i > 0) { |
| 128 | + sumA -= state(i) * f.a(i) |
| 129 | + sumB += state(i) * f.b(i) |
| 130 | + state(i) = state(i - 1) |
| 131 | + i -= 1 |
| 132 | + } |
| 133 | + state(1) = sumA |
| 134 | + sumA * f.b(0) + sumB |
| 135 | + } |
| 136 | + |
| 137 | + /** |
| 138 | + * IIR filter. The multiplication factor is present to compensate for the loss in energy |
| 139 | + * caused by interpolation. |
| 140 | + * @param x Input sequence |
| 141 | + * @param mult A multiplication factor by which the output is multiplied. |
| 142 | + * @return Filtered sequence using the configured filter parameters |
| 143 | + */ |
| 144 | + def filter(x: Array[Float], mult: Float = 1.0f): Array[Float] = { |
| 145 | + val y = new Array[Float](x.length) |
| 146 | + var n = 0 |
| 147 | + while (n < x.length) { |
| 148 | + y(n) = filterOne(x(n)) * mult |
| 149 | + n += 1 |
| 150 | + } |
| 151 | + y |
| 152 | + } |
| 153 | +} |
| 154 | + |
| 155 | +trait Multirate { |
| 156 | + def process(x: Array[Float]): Array[Float] |
| 157 | +} |
| 158 | + |
| 159 | +/** |
| 160 | + * Upsample a signal by the factor specified. If a signal originally sampled at 8kHz is upsampled by a factor |
| 161 | + * of 6, the returned signal will have a sample frequency of 48kHz and will retain its original bandwidth |
| 162 | + * @param factor The interpolation factor to use |
| 163 | + */ |
| 164 | +case class Upsampler(factor: Int) extends FilterIIR(MultirateFilterFactory(factor)) with Interpolator with Multirate { |
| 165 | + |
| 166 | + /** |
| 167 | + * Process a signal, increasing its effective sample rate |
| 168 | + * @param x Signal to be upsampled |
| 169 | + */ |
| 170 | + def process(x: Array[Float]): Array[Float] = filter(interpolate(x), factor) |
| 171 | +} |
| 172 | + |
| 173 | +/** |
| 174 | + * Reduce a signal's sample rate by first filtering it to remove all potential alias frequencies then down-sampling |
| 175 | + * (decimating) the signal by the specified factor. If a signal is sampled at 48kHz and a decimation factor of |
| 176 | + * 6 is specified, the output signal will have a sample rate of 8 kHz |
| 177 | + * @param factor Decimation factor. |
| 178 | + */ |
| 179 | +case class Downsampler(factor: Int) extends FilterIIR(MultirateFilterFactory(factor)) with Decimator with Multirate { |
| 180 | + |
| 181 | + /** |
| 182 | + * Process the signal, effectively decreasing its sample rate. |
| 183 | + * @param x Input audio signal |
| 184 | + */ |
| 185 | + def process(x: Array[Float]): Array[Float] = decimate(filter(x)) |
| 186 | +} |
0 commit comments