Add audio resampling util (#41)

* Add resampling util

* fmt

CHANGELOG.md

@@ -1,3 +1,7 @@
+## 1.5.0
+
+- Added an audio resampling util.
+
 ## 1.4.1
 
 - Special case FFTs of size 4 and 5, which are base cases of composite FFT. This

lib/fftea.dart

@@ -14,5 +14,6 @@
 
 export 'conv.dart';
 export 'impl.dart' show FFT;
+export 'resample.dart';
 export 'stft.dart';
 export 'util.dart' show ComplexArray;

lib/impl.dart

@@ -153,8 +153,9 @@ abstract class FFT {
   /// Performs an inverse FFT and discards the imaginary components of the
   /// result. Returns a newly allocated [Float64List].
   ///
-  /// This method expects the full result of [realFft], so don't use
-  /// [ComplexArray.discardConjugates] if you need to call [realInverseFft].
+  /// This method expects the full result of [realFft], so if you use
+  /// [ComplexArray.discardConjugates], remember to use
+  /// [ComplexArray.createConjugates] before calling [realInverseFft].
   ///
   /// WARINING: For efficiency reasons, this modifies [complexArray]. If you
   /// need the original values in [complexArray] to remain unmodified, make a

lib/resample.dart

@@ -0,0 +1,75 @@
+// Copyright 2023 The fftea authors
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+import 'dart:typed_data';
+
+import 'impl.dart' show FFT;
+import 'util.dart' show ComplexArray;
+
+/// Resamples the [input] audio to the given [outputLength].
+///
+/// Returns a new array of length [outputLength] containing the resampled
+/// [input]. Doesn't modify the input array.
+///
+/// This function FFTs the input, truncates or zero pads the frequencies to the
+/// output length, then IFFTs to get the output. This isn't the best way of
+/// resampling audio. It's intended to be simple and good enough for most
+/// purposes. A more typical approach is convolution with a windowed sinc
+/// function, which will often be more efficient and produce better results, but
+/// requires a bit more design work to match the parameters to the use case. If
+/// you just want something that works well enough, this function is a good
+/// starting point.
+Float64List resample(List<double> input, int outputLength) {
+  if (input.length == outputLength) {
+    return Float64List.fromList(input);
+  }
+  final inf = FFT(input.length).realFft(input).discardConjugates();
+  final outflen = ComplexArray.discardConjugatesLength(outputLength);
+  late Float64x2List outf;
+  if (outflen < inf.length) {
+    // Truncate.
+    outf = Float64x2List.sublistView(inf, 0, outflen);
+  } else {
+    // Zero pad.
+    outf = Float64x2List(outflen);
+    for (int i = 0; i < inf.length; ++i) {
+      outf[i] = inf[i];
+    }
+  }
+  final out =
+      FFT(outputLength).realInverseFft(outf.createConjugates(outputLength));
+  // Resampling like this changes the amplitude, so we need to fix that.
+  final ratio = outputLength.toDouble() / input.length;
+  for (int i = 0; i < out.length; ++i) {
+    out[i] *= ratio;
+  }
+  return out;
+}
+
+/// Resamples the [input] audio by the given sampling [ratio]. If [ratio] > 1
+/// the result will have more samples.
+///
+/// See [resample] for more information.
+Float64List resampleByRatio(List<double> input, double ratio) =>
+    resample(input, (input.length * ratio).round());
+
+/// Resamples the [input] audio from [inputSampleRate] to [outputSampleRate].
+///
+/// See [resample] for more information.
+Float64List resampleByRate(
+  List<double> input,
+  double inputSampleRate,
+  double outputSampleRate,
+) =>
+    resampleByRatio(input, outputSampleRate / inputSampleRate);

lib/util.dart

@@ -81,7 +81,9 @@ extension ComplexArray on Float64x2List {
     }
   }
 
-  static int _discConjLen(int length) {
+  /// Returns the length of the result of calling [discardConjugates] on a list
+  /// of the given [length]. See [discardConjugates] for more information.
+  static int discardConjugatesLength(int length) {
     return (length == 0) ? 0 : ((length >>> 1) + 1);
   }
 
@@ -100,7 +102,7 @@ extension ComplexArray on Float64x2List {
   /// This method returns a new array (which is a view into the same data). It
   /// does not modify this array, or make a copy of the data.
   Float64x2List discardConjugates() {
-    return Float64x2List.sublistView(this, 0, _discConjLen(length));
+    return Float64x2List.sublistView(this, 0, discardConjugatesLength(length));
   }
 
   /// Creates redundant conjugate terms. This is the inverse of
@@ -131,7 +133,7 @@ extension ComplexArray on Float64x2List {
   /// This method returns a totally new array containing a copy of this array,
   /// with the extra values appended at the end.
   Float64x2List createConjugates(int outputLength) {
-    if (_discConjLen(outputLength) != length) {
+    if (discardConjugatesLength(outputLength) != length) {
       throw ArgumentError(
         'Output length must be either (2 * length - 2) or (2 * length - 1).',
         'outputLength',

pubspec.yaml

@@ -13,7 +13,7 @@
 # limitations under the License.
 
 name: fftea
-version: 1.4.1
+version: 1.5.0
 description: >
   Fast Fourier Transform (FFT) library that can handle inputs of any size.
   Includes related tools such as STFT and convolution.

test/generate_test.py

@@ -278,6 +278,35 @@ def maker():
 
   write('}\n')
 
+def generateResamp(write, sizes):
+  write(kPreamble)
+
+  for insize, outsize in sizes:
+    matfile = 'test/data/resample_%d_%d.mat' % (insize, outsize)
+    def includeFreq(f, size):
+      return f * 2 <= size
+    def addFreq(a, f):
+      if not includeFreq(f, len(a)):
+        return
+      dt = f * 2 * math.pi / len(a)
+      for i in range(len(a)):
+        a[i] += math.sin(i * dt) / math.sqrt(f)
+    def maker():
+      a = [0] * insize
+      b = [0] * outsize
+      f = 1
+      while includeFreq(f, insize) and includeFreq(f, outsize):
+        addFreq(a, f)
+        addFreq(b, f)
+        f = math.ceil(f * 1.5)
+      return [a, b]
+    createDataset(matfile, maker)
+    write("  test('Resample %d %d', () async {" % (insize, outsize))
+    write("    await testResample('%s');" % (matfile))
+    write('  });\n')
+
+  write('}\n')
+
 def run(gen, testName, *args):
   outFile = os.path.normpath(os.path.join(
       os.path.dirname(__file__), testName + '_generated_test.dart'))
@@ -304,5 +333,7 @@ def run(gen, testName, *args):
     (1024, 1024, 1024), (2000, 3000, 1400), (123, 456, None), (456, 789, None),
     (1234, 1234, None)],
     [(1, 1), (4, 4), (5, 47), (91, 12), (127, 129), (337, 321), (1024, 1024)])
+run(generateResamp, 'resample',
+    [(1000, outsize) for outsize in [100, 300, 800, 1000, 1500, 2000, 5000]])
 run(generateMisc, 'misc')
 print('Done :)')

test/resample_generated_test.dart

@@ -0,0 +1,57 @@
+// Copyright 2022 The fftea authors
+//
+// Licensed under the Apache License, Version 2.0 (the "License");
+// you may not use this file except in compliance with the License.
+// You may obtain a copy of the License at
+//
+//     https://www.apache.org/licenses/LICENSE-2.0
+//
+// Unless required by applicable law or agreed to in writing, software
+// distributed under the License is distributed on an "AS IS" BASIS,
+// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+// See the License for the specific language governing permissions and
+// limitations under the License.
+
+// GENERATED FILE. DO NOT EDIT.
+
+// Test cases generated with numpy as a reference implementation, using:
+//   python3 test/generate_test.py && dart format .
+
+// ignore_for_file: unused_import
+// ignore_for_file: require_trailing_commas
+
+import 'package:fftea/fftea.dart';
+import 'package:fftea/impl.dart';
+import 'package:test/test.dart';
+
+import 'test_util.dart';
+
+void main() {
+  test('Resample 1000 100', () async {
+    await testResample('test/data/resample_1000_100.mat');
+  });
+
+  test('Resample 1000 300', () async {
+    await testResample('test/data/resample_1000_300.mat');
+  });
+
+  test('Resample 1000 800', () async {
+    await testResample('test/data/resample_1000_800.mat');
+  });
+
+  test('Resample 1000 1000', () async {
+    await testResample('test/data/resample_1000_1000.mat');
+  });
+
+  test('Resample 1000 1500', () async {
+    await testResample('test/data/resample_1000_1500.mat');
+  });
+
+  test('Resample 1000 2000', () async {
+    await testResample('test/data/resample_1000_2000.mat');
+  });
+
+  test('Resample 1000 5000', () async {
+    await testResample('test/data/resample_1000_5000.mat');
+  });
+}

test/test_util.dart

@@ -140,6 +140,18 @@ Future<void> testLinConv(String filename) async {
   expectClose(result2, raw[2]);
 }
 
+Future<void> testResample(String filename) async {
+  final raw = await readMatFile(filename);
+  expect(raw.length, 2);
+  final result1 = resample(raw[0], raw[1].length);
+  expectClose(result1, raw[1]);
+  final ratio = raw[1].length / raw[0].length;
+  final result2 = resampleByRatio(raw[0], ratio);
+  expectClose(result2, raw[1]);
+  final result3 = resampleByRate(raw[0], 44100, ratio * 44100);
+  expectClose(result3, raw[1]);
+}
+
 Future<List<List<double>>> readMatFile(String filename) async {
   final bytes = await File(filename).readAsBytes();
   int p = 0;