From 0664803e1d30e53b76064e9d5d7e9f9cd1d224b7 Mon Sep 17 00:00:00 2001 From: Edwin Jakobs Date: Sun, 17 Mar 2024 16:01:03 +0100 Subject: [PATCH] [orx-fft] Add orx-fft for a simple fast fourier transform routine --- orx-fft/README.md | 5 + orx-fft/build.gradle.kts | 23 ++ orx-fft/src/commonMain/kotlin/FFT.kt | 225 ++++++++++++++++++ orx-fft/src/commonMain/kotlin/HannWindow.kt | 9 + .../src/commonMain/kotlin/IdentityWindow.kt | 6 + .../src/commonMain/kotlin/WindowFunction.kt | 39 +++ orx-fft/src/jvmDemo/kotlin/DemoFFTShape01.kt | 109 +++++++++ 7 files changed, 416 insertions(+) create mode 100644 orx-fft/README.md create mode 100644 orx-fft/build.gradle.kts create mode 100644 orx-fft/src/commonMain/kotlin/FFT.kt create mode 100644 orx-fft/src/commonMain/kotlin/HannWindow.kt create mode 100644 orx-fft/src/commonMain/kotlin/IdentityWindow.kt create mode 100644 orx-fft/src/commonMain/kotlin/WindowFunction.kt create mode 100644 orx-fft/src/jvmDemo/kotlin/DemoFFTShape01.kt diff --git a/orx-fft/README.md b/orx-fft/README.md new file mode 100644 index 00000000..38ab0299 --- /dev/null +++ b/orx-fft/README.md @@ -0,0 +1,5 @@ +# orx-fft + +Simple forward and inverse FFT routine + +The FFT routine found in `orx-fft` is a Kotlin port of Minim's FFT routine. diff --git a/orx-fft/build.gradle.kts b/orx-fft/build.gradle.kts new file mode 100644 index 00000000..5c6a74f8 --- /dev/null +++ b/orx-fft/build.gradle.kts @@ -0,0 +1,23 @@ +plugins { + org.openrndr.extra.convention.`kotlin-multiplatform` +} + +kotlin { + sourceSets { + @Suppress("UNUSED_VARIABLE") + val commonMain by getting { + dependencies { + + } + } + + @Suppress("UNUSED_VARIABLE") + val jvmDemo by getting { + dependencies { + implementation(project(":orx-shapes")) + implementation(project(":orx-noise")) + + } + } + } +} \ No newline at end of file diff --git a/orx-fft/src/commonMain/kotlin/FFT.kt b/orx-fft/src/commonMain/kotlin/FFT.kt new file mode 100644 index 00000000..6bda6de5 --- /dev/null +++ b/orx-fft/src/commonMain/kotlin/FFT.kt @@ -0,0 +1,225 @@ +package org.openrndr.extra.fft + +import kotlin.math.* + +/* +Based on https://github.com/ddf/Minim/blob/e294e2881a20340603ee0156cb9188c15b5915c2/src/main/java/ddf/minim/analysis/FFT.java +I (EJ) stripped away spectrum and averages. + +This is the original license (GPLv2). I am not sure if my low-effort Kotlin port falls under the same license. + + * Copyright (c) 2007 - 2008 by Damien Di Fede + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU Library General Public License as published + * by the Free Software Foundation; either version 2 of the License, or + * (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU Library General Public License for more details. + * + * You should have received a copy of the GNU Library General Public + * License along with this program; if not, write to the Free Software + * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. + */ + +class FFT(val size: Int, private val windowFunction: WindowFunction = IdentityWindow()) { + + var real = FloatArray(size) + var imag = FloatArray(size) + + private fun setComplex(r: FloatArray, i: FloatArray) { + if (real.size != r.size && imag.size != i.size) { + error("FourierTransform.setComplex: the two arrays must be the same length as their member counterparts.") + } else { + r.copyInto(real) + i.copyInto(imag) + } + } + + fun magnitudeSum(includeDC: Boolean = false): Double { + var sum = 0.0 + for (i in (if (includeDC) 0 else 1)..size / 2) { + sum += magnitude(i) + } + return sum + } + + fun scaleAll(sr: Float, includeDC: Boolean = false) { + for (i in (if (includeDC) 0 else 1)..size / 2) { + scaleBand(i, sr) + } + } + + fun magnitude(i: Int): Float { + return sqrt(real[i] * real[i] + imag[i] * imag[i]) + } + + fun phase(i: Int): Float { + return atan2(imag[i], real[i]) + } + + fun shiftPhase(i: Int, shift: Double) { + val m = magnitude(i) + val phase = phase(i) + real[i] = (cos(phase + shift) * m).toFloat() + imag[i] = (sin(phase + shift) * m).toFloat() + + if (i != 0 && i != size / 2) { + real[(size - i)] = real[i] + imag[(size - i)] = -imag[i] + } + } + + fun scaleBand(i: Int, sr: Float) { + if (sr < 0) { + error("Can't scale a frequency band by a negative value.") + } + + real[i] *= sr + imag[i] *= sr + + if (i != 0 && i != size / 2) { + real[(size - i)] = real[i] + imag[(size - i)] = -imag[i] + } + } + + + // performs an in-place fft on the data in the real and imag arrays + // bit reversing is not necessary as the data will already be bit reversed + private fun fft() { + var halfSize = 1 + while (halfSize < real.size) { + // float k = -(float)Math.PI/halfSize; + // phase shift step + // float phaseShiftStepR = (float)Math.cos(k); + // float phaseShiftStepI = (float)Math.sin(k); + // using lookup table + val phaseShiftStepR = cos[halfSize] + val phaseShiftStepI = sin[halfSize] + // current phase shift + var currentPhaseShiftR = 1.0f + var currentPhaseShiftI = 0.0f + for (fftStep in 0 until halfSize) { + var i = fftStep + while (i < real.size) { + val off = i + halfSize + val tr = (currentPhaseShiftR * real[off]) - (currentPhaseShiftI * imag[off]) + val ti = (currentPhaseShiftR * imag[off]) + (currentPhaseShiftI * real[off]) + real[off] = real[i] - tr + imag[off] = imag[i] - ti + real[i] += tr + imag[i] += ti + i += 2 * halfSize + } + val tmpR = currentPhaseShiftR + currentPhaseShiftR = (tmpR * phaseShiftStepR) - (currentPhaseShiftI * phaseShiftStepI) + currentPhaseShiftI = (tmpR * phaseShiftStepI) + (currentPhaseShiftI * phaseShiftStepR) + } + halfSize *= 2 + } + } + + private fun doWindow(samples: FloatArray) { + windowFunction.apply(samples) + } + + + fun forward(buffer: FloatArray) { + if (buffer.size != size) { + error("FFT.forward: The length of the passed sample buffer must be equal to timeSize().") + } + doWindow(buffer) + // copy samples to real/imag in bit-reversed order + bitReverseSamples(buffer, 0) + // perform the fft + fft() + } + + fun forward(buffer: FloatArray, startAt: Int) { + if (buffer.size - startAt < size) { + error( + "FourierTransform.forward: not enough samples in the buffer between " + + startAt + " and " + buffer.size + " to perform a transform." + ) + } + windowFunction.apply(buffer, startAt, size) + bitReverseSamples(buffer, startAt) + fft() + } + + /** + * Performs a forward transform on the passed buffers. + * + * @param buffReal the real part of the time domain signal to transform + * @param buffImag the imaginary part of the time domain signal to transform + */ + fun forward(buffReal: FloatArray, buffImag: FloatArray) { + if (buffReal.size != size || buffImag.size != size) { + error("FFT.forward: The length of the passed buffers must be equal to timeSize().") + } + setComplex(buffReal, buffImag) + bitReverseComplex() + fft() + } + + fun inverse(buffer: FloatArray) { + if (buffer.size > real.size) { + error("FFT.inverse: the passed array's length must equal FFT.timeSize().") + } + // conjugate + for (i in 0 until size) { + imag[i] *= -1.0f + } + bitReverseComplex() + fft() + // copy the result in real into buffer, scaling as we do + for (i in buffer.indices) { + buffer[i] = real[i] / real.size + } + } + + private val reverse by lazy { buildReverseTable() } + + private fun buildReverseTable(): IntArray { + val reverse = IntArray(size) + // set up the bit reversing table + reverse[0] = 0 + var limit = 1 + var bit = size / 2 + while (limit < size) { + for (i in 0 until limit) reverse[i + limit] = reverse[i] + bit + limit = limit shl 1 + bit = bit shr 1 + } + return reverse + } + + // copies the values in the samples array into the real array + // in bit reversed order. the imag array is filled with zeros. + private fun bitReverseSamples(samples: FloatArray, startAt: Int) { + for (i in 0 until size) { + real[i] = samples[startAt + reverse[i]] + imag[i] = 0.0f + } + } + + // bit reverse real[] and imag[] + private fun bitReverseComplex() { + val revReal = FloatArray(real.size) + val revImag = FloatArray(imag.size) + for (i in real.indices) { + revReal[i] = real[reverse[i]] + revImag[i] = imag[reverse[i]] + } + real = revReal + imag = revImag + } + + // lookup tables + private val sin by lazy { FloatArray(size) { i -> sin((-PI.toFloat() / i).toDouble()).toFloat() } } + private val cos by lazy { FloatArray(size) { i -> cos((-PI.toFloat() / i).toDouble()).toFloat() } } +} \ No newline at end of file diff --git a/orx-fft/src/commonMain/kotlin/HannWindow.kt b/orx-fft/src/commonMain/kotlin/HannWindow.kt new file mode 100644 index 00000000..87a2c429 --- /dev/null +++ b/orx-fft/src/commonMain/kotlin/HannWindow.kt @@ -0,0 +1,9 @@ +package org.openrndr.extra.fft + +import kotlin.math.PI +import kotlin.math.cos + +class HannWindow : WindowFunction() { + override fun value(length: Int, index: Int): Float = 0.5f * (1f - cos((PI * 2.0 * index / (length - 1f))) + .toFloat()) +} \ No newline at end of file diff --git a/orx-fft/src/commonMain/kotlin/IdentityWindow.kt b/orx-fft/src/commonMain/kotlin/IdentityWindow.kt new file mode 100644 index 00000000..a5949796 --- /dev/null +++ b/orx-fft/src/commonMain/kotlin/IdentityWindow.kt @@ -0,0 +1,6 @@ +package org.openrndr.extra.fft + +class IdentityWindow + : WindowFunction() { + override fun value(length: Int, index: Int): Float = 1.0f +} \ No newline at end of file diff --git a/orx-fft/src/commonMain/kotlin/WindowFunction.kt b/orx-fft/src/commonMain/kotlin/WindowFunction.kt new file mode 100644 index 00000000..71939641 --- /dev/null +++ b/orx-fft/src/commonMain/kotlin/WindowFunction.kt @@ -0,0 +1,39 @@ +package org.openrndr.extra.fft + +abstract class WindowFunction { + private var length: Int = 0 + + /** + * Apply the window function to a sample buffer. + * + * @param samples a sample buffer + */ + fun apply(samples: FloatArray) { + this.length = samples.size + + for (n in samples.indices) { + samples[n] *= value(samples.size, n) + } + } + + /** + * Apply the window to a portion of this sample buffer, + * given an offset from the beginning of the buffer + * and the number of samples to be windowed. + * + * @param samples + * float[]: the array of samples to apply the window to + * @param offset + * int: the index in the array to begin windowing + * @param length + * int: how many samples to apply the window to + */ + fun apply(samples: FloatArray, offset: Int, length: Int) { + this.length = length + + for (n in offset until offset + length) { + samples[n] *= value(length, n - offset) + } + } + protected abstract fun value(length: Int, index: Int): Float +} \ No newline at end of file diff --git a/orx-fft/src/jvmDemo/kotlin/DemoFFTShape01.kt b/orx-fft/src/jvmDemo/kotlin/DemoFFTShape01.kt new file mode 100644 index 00000000..74beb4e4 --- /dev/null +++ b/orx-fft/src/jvmDemo/kotlin/DemoFFTShape01.kt @@ -0,0 +1,109 @@ +import org.openrndr.application + +import org.openrndr.color.ColorRGBa +import org.openrndr.extra.fft.FFT +import org.openrndr.extra.noise.scatter +import org.openrndr.extra.shapes.hobbycurve.hobbyCurve +import org.openrndr.math.Vector2 +import org.openrndr.extra.shapes.splines.catmullRom +import org.openrndr.extra.shapes.splines.toContour +import org.openrndr.math.smoothstep +import org.openrndr.math.transforms.buildTransform +import kotlin.math.* +import kotlin.random.Random + +/** + * Demonstration of using FFT to filter a two-dimensional shape. Mouse xy-position is mapped + * to lowpass and highpass settings of the filter. + */ +fun main() { + application { + configure { + width = 720 + height = 720 + } + program { + val fftSize = 512 + val fft = FFT(fftSize) + fun List.toFloatArrays(x: FloatArray, y: FloatArray) { + for ((index, segment) in this.withIndex()) { + x[index] = segment.x.toFloat() + y[index] = segment.y.toFloat() + } + } + + fun vectorsFromFloatArrays(x: FloatArray, y: FloatArray): List { + val n = x.size + val result = mutableListOf() + for (i in 0 until n) { + result.add(Vector2(x[i].toDouble(), y[i].toDouble())) + } + return result + } + + fun lp(t: Double, c: Double): Double { + return smoothstep(c, c - 0.1, t) + } + + fun hp(t: Double, c: Double): Double { + return smoothstep(c, c + 0.1, t) + } + + val c = hobbyCurve( + drawer.bounds.scatter(40.0, distanceToEdge = 100.0, random = Random(0)), + true + ).transform(buildTransform { translate(-drawer.bounds.center) }) + + val x = FloatArray(fftSize) + val y = FloatArray(fftSize) + + val xFiltered = FloatArray(fftSize) + val yFiltered = FloatArray(fftSize) + + extend { + c.equidistantPositions(fftSize).take(fftSize).toFloatArrays(x, y) + + // process x-component + fft.forward(x) + val xpower = fft.magnitudeSum() + + val lpc = mouse.position.x / width + val hpc = mouse.position.y / height + + for (i in 1..fftSize / 2) { + val t = i.toDouble() / (fftSize / 2 - 1) + val f = max(lp(t, lpc), hp(t, hpc)) + fft.scaleBand(i, f.toFloat()) + } + val xfpower = fft.magnitudeSum().coerceAtLeast(1.0) + + fft.scaleAll((xpower / xfpower).toFloat()) + fft.inverse(xFiltered) + + // process y-component + fft.forward(y) + val ypower = fft.magnitudeSum() + + for (i in 1..fftSize / 2) { + val t = i.toDouble() / (fftSize / 2 - 1) + val f = max(lp(t, lpc), hp(t, hpc)) + fft.scaleBand(i, f.toFloat()) + } + val yfpower = fft.magnitudeSum().coerceAtLeast(1.0) + + fft.scaleAll((ypower / yfpower).toFloat()) + fft.inverse(yFiltered) + + val cr = vectorsFromFloatArrays(xFiltered, yFiltered).catmullRom(closed = true).toContour() + //val cr = ShapeContour.fromPoints(vectorsFromFloatArrays(xr, yr), closed=true) + + val recenteredShape = cr.transform(buildTransform { + translate(drawer.bounds.center) + }) + drawer.fill = null + drawer.stroke = ColorRGBa.WHITE + drawer.contour(recenteredShape) + } + } + } +} \ No newline at end of file