From 36bbfd0cae593bce3a446b8a7b431b869a134964 Mon Sep 17 00:00:00 2001 From: nicolas-f <1382241+nicolas-f@users.noreply.github.com> Date: Wed, 15 Nov 2023 10:48:07 +0100 Subject: [PATCH] unit test ok --- .../noisecapture/signal/realfft.kt | 44 ++++++++++++++----- .../noisecapture/signal/realifft.kt | 32 +++++++------- .../noisecapture/signal/TestFFT.kt | 27 +++++++++--- 3 files changed, 71 insertions(+), 32 deletions(-) diff --git a/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realfft.kt b/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realfft.kt index 57599c2..f9881a0 100644 --- a/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realfft.kt +++ b/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realfft.kt @@ -6,11 +6,23 @@ import kotlin.math.log2 import kotlin.math.pow import kotlin.math.sin -fun realFFT(length: Int, realArray: DoubleArray) { - val m = log2(length.toDouble()).toInt() +/** + * The bitwise AND operation of the number and its predecessor (number - 1) should result in 0. + * This is because powers of two in binary have only one bit set, and subtracting 1 from them flips + * that bit and sets all the lower bits to 1. So, the bitwise AND with the original number should + * be 0 if it is a power of two. + */ +fun isPowerOfTwo(number: Int): Boolean { + return number > 0 && (number and (number - 1)) == 0 +} + +fun realFFT(realArray: DoubleArray): DoubleArray { + require(isPowerOfTwo(realArray.size)) + val outArray = DoubleArray(realArray.size + 2) + realArray.copyInto(outArray) + val m = log2(realArray.size.toDouble()).toInt() val n = (2.0.pow(m) + 0.5).toInt() - require(n <= realArray.size) - fft(n/2, realArray) + fft(n/2, outArray) val a = DoubleArray(n) val b = DoubleArray(n) @@ -24,17 +36,25 @@ fun realFFT(length: Int, realArray: DoubleArray) { for (k in 1 until n / 4 + 1) { val k2 = 2 * k val xr = - realArray[k2] * a[k2] - realArray[k2 + 1] * a[k2 + 1] + realArray[n - k2] * b[k2] + realArray[n - k2 + 1] * b[k2 + 1] + outArray[k2] * a[k2] - outArray[k2 + 1] * a[k2 + 1] + outArray[n - k2] * b[k2] + outArray[n - k2 + 1] * b[k2 + 1] val xi = - realArray[k2] * a[k2 + 1] + realArray[k2 + 1] * a[k2] + realArray[n - k2] * b[k2 + 1] - realArray[n - k2 + 1] * b[k2] + outArray[k2] * a[k2 + 1] + outArray[k2 + 1] * a[k2] + outArray[n - k2] * b[k2 + 1] - outArray[n - k2 + 1] * b[k2] val xrN = - realArray[n - k2] * a[n - k2] - realArray[n - k2 + 1] * a[n - k2 + 1] + realArray[k2] * b[n - k2] + realArray[k2 + 1] * b[n - k2 + 1] + outArray[n - k2] * a[n - k2] - outArray[n - k2 + 1] * a[n - k2 + 1] + outArray[k2] * b[n - k2] + outArray[k2 + 1] * b[n - k2 + 1] val xiN = - realArray[n - k2] * a[n - k2 + 1] + realArray[n - k2 + 1] * a[n - k2] + realArray[k2] * b[n - k2 + 1] - realArray[k2 + 1] * b[n - k2] - realArray[k2] = xr - realArray[k2 + 1] = xi - realArray[n - k2] = xrN - realArray[n - k2 + 1] = xiN + outArray[n - k2] * a[n - k2 + 1] + outArray[n - k2 + 1] * a[n - k2] + outArray[k2] * b[n - k2 + 1] - outArray[k2 + 1] * b[n - k2] + outArray[k2] = xr + outArray[k2 + 1] = xi + outArray[n - k2] = xrN + outArray[n - k2 + 1] = xiN } + + val tmp = outArray[0] + outArray[n] = outArray[0] - outArray[1] + outArray[0] = tmp + outArray[1] + outArray[1] = 0.0 + outArray[n+1] = 0.0 + + return outArray } diff --git a/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realifft.kt b/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realifft.kt index f550f83..cd0acba 100644 --- a/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realifft.kt +++ b/NoiseCapture/src/commonMain/kotlin/org/noise_planet/noisecapture/signal/realifft.kt @@ -6,10 +6,11 @@ import kotlin.math.log2 import kotlin.math.pow import kotlin.math.sin -fun realIFFT(length: Int, realArray: DoubleArray) { - val m = log2(length.toDouble()).toInt() +fun realIFFT(realArray: DoubleArray) : DoubleArray { + require(isPowerOfTwo(realArray.size - 2)) + val outArray = realArray.copyOf() + val m = log2((realArray.size - 2).toDouble()).toInt() val n = (2.0.pow(m) + 0.5).toInt() - require(n <= realArray.size) val a = DoubleArray(n) val b = DoubleArray(n) @@ -23,21 +24,22 @@ fun realIFFT(length: Int, realArray: DoubleArray) { for (k in 1 until n / 4 + 1) { val k2 = 2 * k val xr = - realArray[k2] * a[k2] + realArray[k2 + 1] * a[k2 + 1] + realArray[n - k2] * b[k2] - realArray[n - k2 + 1] * b[k2 + 1] + outArray[k2] * a[k2] + outArray[k2 + 1] * a[k2 + 1] + outArray[n - k2] * b[k2] - outArray[n - k2 + 1] * b[k2 + 1] val xi = - -realArray[k2] * a[k2 + 1] + realArray[k2 + 1] * a[k2] - realArray[n - k2] * b[k2 + 1] - realArray[n - k2 + 1] * b[k2] + -outArray[k2] * a[k2 + 1] + outArray[k2 + 1] * a[k2] - outArray[n - k2] * b[k2 + 1] - outArray[n - k2 + 1] * b[k2] val xrN = - realArray[n - k2] * a[n - k2] + realArray[n - k2 + 1] * a[n - k2 + 1] + realArray[k2] * b[n - k2] - realArray[k2 + 1] * b[n - k2 + 1] + outArray[n - k2] * a[n - k2] + outArray[n - k2 + 1] * a[n - k2 + 1] + outArray[k2] * b[n - k2] - outArray[k2 + 1] * b[n - k2 + 1] val xiN = - -realArray[n - k2] * a[n - k2 + 1] + realArray[n - k2 + 1] * a[n - k2] - realArray[k2] * b[n - k2 + 1] - realArray[k2 + 1] * b[n - k2] - realArray[k2] = xr - realArray[k2 + 1] = xi - realArray[n - k2] = xrN - realArray[n - k2 + 1] = xiN + -outArray[n - k2] * a[n - k2 + 1] + outArray[n - k2 + 1] * a[n - k2] - outArray[k2] * b[n - k2 + 1] - outArray[k2 + 1] * b[n - k2] + outArray[k2] = xr + outArray[k2 + 1] = xi + outArray[n - k2] = xrN + outArray[n - k2 + 1] = xiN } - val temp = realArray[0] - realArray[0] = 0.5 * realArray[0] + 0.5 * realArray[n-1] - realArray[1] = 0.5 * temp - 0.5 * realArray[n-1] - iFFT(n/2, realArray) + val temp = outArray[0] + outArray[0] = 0.5 * outArray[0] + 0.5 * outArray[n] + outArray[1] = 0.5 * temp - 0.5 * outArray[n] + iFFT(n/2, outArray) + return outArray.copyOf(outArray.size - 2) } diff --git a/NoiseCapture/src/commonTest/kotlin/org/noise_planet/noisecapture/signal/TestFFT.kt b/NoiseCapture/src/commonTest/kotlin/org/noise_planet/noisecapture/signal/TestFFT.kt index 90373b2..ae3efab 100644 --- a/NoiseCapture/src/commonTest/kotlin/org/noise_planet/noisecapture/signal/TestFFT.kt +++ b/NoiseCapture/src/commonTest/kotlin/org/noise_planet/noisecapture/signal/TestFFT.kt @@ -24,11 +24,28 @@ class TestFFT { val sampleRate = 64.0 // Hz val duration = 2.0.pow(ceil(log2(sampleRate))) / sampleRate val samples = generateSinusoidalSignal(frequency, sampleRate, duration) - val fftResult = samples.copyOf() - realFFT(samples.size, fftResult) - realIFFT(samples.size, fftResult) - fftResult.forEachIndexed { - index, value -> assertEquals(value, samples[index], 1e-8) + val spectrum = realFFT(samples) + val result = realIFFT(spectrum) + assertEquals(samples.size, result.size) + samples.forEachIndexed { + index, value -> assertEquals(value, result[index], 1e-8) + } + } + + @Test + fun testRFFTIncremental() { + val samples = doubleArrayOf(1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0) + val expected = + doubleArrayOf(36.0, 0.0, -4.0, 9.65685425, -4.0, 4.0, -4.0, 1.65685425, -4.0, 0.0) + val result = realFFT(samples) + assertEquals(expected.size, result.size) + expected.forEachIndexed { index, value -> + assertEquals(value, result[index], 1e-8) + } + val origin = realIFFT(result) + assertEquals(samples.size, origin.size) + samples.forEachIndexed { index, value -> + assertEquals(value, origin[index], 1e-8) } } } \ No newline at end of file