diff --git a/app/src/main/java/xdsopl/robot36/FastFourierTransform.java b/app/src/main/java/xdsopl/robot36/FastFourierTransform.java new file mode 100644 index 0000000..f446fc1 --- /dev/null +++ b/app/src/main/java/xdsopl/robot36/FastFourierTransform.java @@ -0,0 +1,353 @@ +/* +Fast Fourier Transform + +Copyright 2025 Ahmet Inan +*/ + +package xdsopl.robot36; + +public class FastFourierTransform { + private final Complex[] tf; + private final Complex tmpA, tmpB, tmpC, tmpD, tmpE, tmpF, tmpG, tmpH, tmpI, tmpJ, tmpK, tmpL, tmpM; + private final Complex tin0, tin1, tin2, tin3, tin4, tin5, tin6; + + FastFourierTransform(int length) { + int rest = length; + while (rest > 1) { + if (rest % 2 == 0) + rest /= 2; + else if (rest % 3 == 0) + rest /= 3; + else if (rest % 5 == 0) + rest /= 5; + else if (rest % 7 == 0) + rest /= 7; + else + break; + } + if (rest != 1) + throw new IllegalArgumentException( + "Transform length must be a composite of 2, 3, 5 and 7, but was: " + + length); + tf = new Complex[length]; + for (int i = 0; i < length; ++i) { + double x = -(2.0 * Math.PI * i) / length; + float a = (float) Math.cos(x); + float b = (float) Math.sin(x); + tf[i] = new Complex(a, b); + } + tmpA = new Complex(); + tmpB = new Complex(); + tmpC = new Complex(); + tmpD = new Complex(); + tmpE = new Complex(); + tmpF = new Complex(); + tmpG = new Complex(); + tmpH = new Complex(); + tmpI = new Complex(); + tmpJ = new Complex(); + tmpK = new Complex(); + tmpL = new Complex(); + tmpM = new Complex(); + tin0 = new Complex(); + tin1 = new Complex(); + tin2 = new Complex(); + tin3 = new Complex(); + tin4 = new Complex(); + tin5 = new Complex(); + tin6 = new Complex(); + } + + private boolean isPowerOfTwo(int n) { + return n > 0 && (n & (n - 1)) == 0; + } + + private boolean isPowerOfFour(int n) { + return isPowerOfTwo(n) && (n & 0x55555555) != 0; + } + + private float cos(int n, int N) { + return (float) Math.cos(n * 2.0 * Math.PI / N); + } + + private float sin(int n, int N) { + return (float) Math.sin(n * 2.0 * Math.PI / N); + } + + private void dft2(Complex out0, Complex out1, Complex in0, Complex in1) { + out0.set(in0).add(in1); + out1.set(in0).sub(in1); + } + + private void radix2(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 2) { + dft2(out[O], out[O + 1], in[I], in[I + S]); + } else { + int Q = N / 2; + dit(out, in, O, I, Q, 2 * S, F); + dit(out, in, O + Q, I + S, Q, 2 * S, F); + for (int k0 = O, k1 = O + Q, l1 = 0; k0 < O + Q; ++k0, ++k1, l1 += S) { + tin1.set(tf[l1]); + if (!F) + tin1.conj(); + tin0.set(out[k0]); + tin1.mul(out[k1]); + dft2(out[k0], out[k1], tin0, tin1); + } + } + } + + private void fwd3(Complex out0, Complex out1, Complex out2, Complex in0, Complex in1, Complex in2) { + tmpA.set(in1).add(in2); + tmpB.set(in1.imag - in2.imag, in2.real - in1.real); + tmpC.set(tmpA).mul(cos(1, 3)); + tmpD.set(tmpB).mul(sin(1, 3)); + out0.set(in0).add(tmpA); + out1.set(in0).add(tmpC).add(tmpD); + out2.set(in0).add(tmpC).sub(tmpD); + } + + private void radix3(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 3) { + if (F) + fwd3(out[O], out[O + 1], out[O + 2], + in[I], in[I + S], in[I + 2 * S]); + else + fwd3(out[O], out[O + 2], out[O + 1], + in[I], in[I + S], in[I + 2 * S]); + } else { + int Q = N / 3; + dit(out, in, O, I, Q, 3 * S, F); + dit(out, in, O + Q, I + S, Q, 3 * S, F); + dit(out, in, O + 2 * Q, I + 2 * S, Q, 3 * S, F); + for (int k0 = O, k1 = O + Q, k2 = O + 2 * Q, l1 = 0, l2 = 0; + k0 < O + Q; ++k0, ++k1, ++k2, l1 += S, l2 += 2 * S) { + tin1.set(tf[l1]); + tin2.set(tf[l2]); + if (!F) { + tin1.conj(); + tin2.conj(); + } + tin0.set(out[k0]); + tin1.mul(out[k1]); + tin2.mul(out[k2]); + if (F) + fwd3(out[k0], out[k1], out[k2], tin0, tin1, tin2); + else + fwd3(out[k0], out[k2], out[k1], tin0, tin1, tin2); + } + } + } + + private void fwd4(Complex out0, Complex out1, Complex out2, Complex out3, + Complex in0, Complex in1, Complex in2, Complex in3) { + tmpA.set(in0).add(in2); + tmpB.set(in0).sub(in2); + tmpC.set(in1).add(in3); + tmpD.set(in1.imag - in3.imag, in3.real - in1.real); + out0.set(tmpA).add(tmpC); + out1.set(tmpB).add(tmpD); + out2.set(tmpA).sub(tmpC); + out3.set(tmpB).sub(tmpD); + } + + private void radix4(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 4) { + if (F) + fwd4(out[O], out[O + 1], out[O + 2], out[O + 3], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S]); + else + fwd4(out[O], out[O + 3], out[O + 2], out[O + 1], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S]); + } else { + int Q = N / 4; + radix4(out, in, O, I, Q, 4 * S, F); + radix4(out, in, O + Q, I + S, Q, 4 * S, F); + radix4(out, in, O + 2 * Q, I + 2 * S, Q, 4 * S, F); + radix4(out, in, O + 3 * Q, I + 3 * S, Q, 4 * S, F); + for (int k0 = O, k1 = O + Q, k2 = O + 2 * Q, k3 = O + 3 * Q, l1 = 0, l2 = 0, l3 = 0; + k0 < O + Q; ++k0, ++k1, ++k2, ++k3, l1 += S, l2 += 2 * S, l3 += 3 * S) { + tin1.set(tf[l1]); + tin2.set(tf[l2]); + tin3.set(tf[l3]); + if (!F) { + tin1.conj(); + tin2.conj(); + tin3.conj(); + } + tin0.set(out[k0]); + tin1.mul(out[k1]); + tin2.mul(out[k2]); + tin3.mul(out[k3]); + if (F) + fwd4(out[k0], out[k1], out[k2], out[k3], tin0, tin1, tin2, tin3); + else + fwd4(out[k0], out[k3], out[k2], out[k1], tin0, tin1, tin2, tin3); + } + } + } + + private void fwd5(Complex out0, Complex out1, Complex out2, Complex out3, Complex out4, + Complex in0, Complex in1, Complex in2, Complex in3, Complex in4) { + tmpA.set(in1).add(in4); + tmpB.set(in2).add(in3); + tmpC.set(in1.imag - in4.imag, in4.real - in1.real); + tmpD.set(in2.imag - in3.imag, in3.real - in2.real); + tmpF.set(tmpA).mul(cos(1, 5)).add(tmpE.set(tmpB).mul(cos(2, 5))); + tmpG.set(tmpC).mul(sin(1, 5)).add(tmpE.set(tmpD).mul(sin(2, 5))); + tmpH.set(tmpA).mul(cos(2, 5)).add(tmpE.set(tmpB).mul(cos(1, 5))); + tmpI.set(tmpC).mul(sin(2, 5)).sub(tmpE.set(tmpD).mul(sin(1, 5))); + out0.set(in0).add(tmpA).add(tmpB); + out1.set(in0).add(tmpF).add(tmpG); + out2.set(in0).add(tmpH).add(tmpI); + out3.set(in0).add(tmpH).sub(tmpI); + out4.set(in0).add(tmpF).sub(tmpG); + } + + private void radix5(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 5) { + if (F) + fwd5(out[O], out[O + 1], out[O + 2], out[O + 3], out[O + 4], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S], in[I + 4 * S]); + else + fwd5(out[O], out[O + 4], out[O + 3], out[O + 2], out[O + 1], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S], in[I + 4 * S]); + } else { + int Q = N / 5; + dit(out, in, O, I, Q, 5 * S, F); + dit(out, in, O + Q, I + S, Q, 5 * S, F); + dit(out, in, O + 2 * Q, I + 2 * S, Q, 5 * S, F); + dit(out, in, O + 3 * Q, I + 3 * S, Q, 5 * S, F); + dit(out, in, O + 4 * Q, I + 4 * S, Q, 5 * S, F); + for (int k0 = O, k1 = O + Q, k2 = O + 2 * Q, k3 = O + 3 * Q, k4 = O + 4 * Q, l1 = 0, l2 = 0, l3 = 0, l4 = 0; + k0 < O + Q; ++k0, ++k1, ++k2, ++k3, ++k4, l1 += S, l2 += 2 * S, l3 += 3 * S, l4 += 4 * S) { + tin1.set(tf[l1]); + tin2.set(tf[l2]); + tin3.set(tf[l3]); + tin4.set(tf[l4]); + if (!F) { + tin1.conj(); + tin2.conj(); + tin3.conj(); + tin4.conj(); + } + tin0.set(out[k0]); + tin1.mul(out[k1]); + tin2.mul(out[k2]); + tin3.mul(out[k3]); + tin4.mul(out[k4]); + if (F) + fwd5(out[k0], out[k1], out[k2], out[k3], out[k4], tin0, tin1, tin2, tin3, tin4); + else + fwd5(out[k0], out[k4], out[k3], out[k2], out[k1], tin0, tin1, tin2, tin3, tin4); + } + } + } + + private void fwd7(Complex out0, Complex out1, Complex out2, Complex out3, Complex out4, Complex out5, Complex out6, + Complex in0, Complex in1, Complex in2, Complex in3, Complex in4, Complex in5, Complex in6) { + tmpA.set(in1).add(in6); + tmpB.set(in2).add(in5); + tmpC.set(in3).add(in4); + tmpD.set(in1.imag - in6.imag, in6.real - in1.real); + tmpE.set(in2.imag - in5.imag, in5.real - in2.real); + tmpF.set(in3.imag - in4.imag, in4.real - in3.real); + tmpH.set(tmpA).mul(cos(1, 7)).add(tmpG.set(tmpB).mul(cos(2, 7))).add(tmpG.set(tmpC).mul(cos(3, 7))); + tmpI.set(tmpD).mul(sin(1, 7)).add(tmpG.set(tmpE).mul(sin(2, 7))).add(tmpG.set(tmpF).mul(sin(3, 7))); + tmpJ.set(tmpA).mul(cos(2, 7)).add(tmpG.set(tmpB).mul(cos(3, 7))).add(tmpG.set(tmpC).mul(cos(1, 7))); + tmpK.set(tmpD).mul(sin(2, 7)).sub(tmpG.set(tmpE).mul(sin(3, 7))).sub(tmpG.set(tmpF).mul(sin(1, 7))); + tmpL.set(tmpA).mul(cos(3, 7)).add(tmpG.set(tmpB).mul(cos(1, 7))).add(tmpG.set(tmpC).mul(cos(2, 7))); + tmpM.set(tmpD).mul(sin(3, 7)).sub(tmpG.set(tmpE).mul(sin(1, 7))).add(tmpG.set(tmpF).mul(sin(2, 7))); + out0.set(in0).add(tmpA).add(tmpB).add(tmpC); + out1.set(in0).add(tmpH).add(tmpI); + out2.set(in0).add(tmpJ).add(tmpK); + out3.set(in0).add(tmpL).add(tmpM); + out4.set(in0).add(tmpL).sub(tmpM); + out5.set(in0).add(tmpJ).sub(tmpK); + out6.set(in0).add(tmpH).sub(tmpI); + } + + private void radix7(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 7) { + if (F) + fwd7(out[O], out[O + 1], out[O + 2], out[O + 3], out[O + 4], out[O + 5], out[O + 6], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S], in[I + 4 * S], in[I + 5 * S], in[I + 6 * S]); + else + fwd7(out[O], out[O + 6], out[O + 5], out[O + 4], out[O + 3], out[O + 2], out[O + 1], + in[I], in[I + S], in[I + 2 * S], in[I + 3 * S], in[I + 4 * S], in[I + 5 * S], in[I + 6 * S]); + } else { + int Q = N / 7; + dit(out, in, O, I, Q, 7 * S, F); + dit(out, in, O + Q, I + S, Q, 7 * S, F); + dit(out, in, O + 2 * Q, I + 2 * S, Q, 7 * S, F); + dit(out, in, O + 3 * Q, I + 3 * S, Q, 7 * S, F); + dit(out, in, O + 4 * Q, I + 4 * S, Q, 7 * S, F); + dit(out, in, O + 5 * Q, I + 5 * S, Q, 7 * S, F); + dit(out, in, O + 6 * Q, I + 6 * S, Q, 7 * S, F); + for (int k0 = O, k1 = O + Q, k2 = O + 2 * Q, k3 = O + 3 * Q, k4 = O + 4 * Q, k5 = O + 5 * Q, k6 = O + 6 * Q, l1 = 0, l2 = 0, l3 = 0, l4 = 0, l5 = 0, l6 = 0; + k0 < O + Q; ++k0, ++k1, ++k2, ++k3, ++k4, ++k5, ++k6, l1 += S, l2 += 2 * S, l3 += 3 * S, l4 += 4 * S, l5 += 5 * S, l6 += 6 * S) { + tin1.set(tf[l1]); + tin2.set(tf[l2]); + tin3.set(tf[l3]); + tin4.set(tf[l4]); + tin5.set(tf[l5]); + tin6.set(tf[l6]); + if (!F) { + tin1.conj(); + tin2.conj(); + tin3.conj(); + tin4.conj(); + tin5.conj(); + tin6.conj(); + } + tin0.set(out[k0]); + tin1.mul(out[k1]); + tin2.mul(out[k2]); + tin3.mul(out[k3]); + tin4.mul(out[k4]); + tin5.mul(out[k5]); + tin6.mul(out[k6]); + if (F) + fwd7(out[k0], out[k1], out[k2], out[k3], out[k4], out[k5], out[k6], tin0, tin1, tin2, tin3, tin4, tin5, tin6); + else + fwd7(out[k0], out[k6], out[k5], out[k4], out[k3], out[k2], out[k1], tin0, tin1, tin2, tin3, tin4, tin5, tin6); + } + } + } + + private void dit(Complex[] out, Complex[] in, int O, int I, int N, int S, boolean F) { + if (N == 1) + out[O].set(in[I]); + else if (isPowerOfFour(N)) + radix4(out, in, O, I, N, S, F); + else if (N % 7 == 0) + radix7(out, in, O, I, N, S, F); + else if (N % 5 == 0) + radix5(out, in, O, I, N, S, F); + else if (N % 3 == 0) + radix3(out, in, O, I, N, S, F); + else if (N % 2 == 0) + radix2(out, in, O, I, N, S, F); + } + + void forward(Complex[] out, Complex[] in) { + if (in.length != tf.length) + throw new IllegalArgumentException("Input array length (" + in.length + + ") must be equal to Transform length (" + tf.length + ")"); + if (out.length != tf.length) + throw new IllegalArgumentException("Output array length (" + out.length + + ") must be equal to Transform length (" + tf.length + ")"); + dit(out, in, 0, 0, tf.length, 1, true); + } + + @SuppressWarnings("unused") + void backward(Complex[] out, Complex[] in) { + if (in.length != tf.length) + throw new IllegalArgumentException("Input array length (" + in.length + + ") must be equal to Transform length (" + tf.length + ")"); + if (out.length != tf.length) + throw new IllegalArgumentException("Output array length (" + out.length + + ") must be equal to Transform length (" + tf.length + ")"); + dit(out, in, 0, 0, tf.length, 1, false); + } +} diff --git a/app/src/main/java/xdsopl/robot36/Hann.java b/app/src/main/java/xdsopl/robot36/Hann.java new file mode 100644 index 0000000..892046a --- /dev/null +++ b/app/src/main/java/xdsopl/robot36/Hann.java @@ -0,0 +1,13 @@ +/* +Hann window + +Copyright 2025 Ahmet Inan +*/ + +package xdsopl.robot36; + +public class Hann { + static double window(int n, int N) { + return 0.5 * (1.0 - Math.cos((2.0 * Math.PI * n) / (N - 1))); + } +} \ No newline at end of file diff --git a/app/src/main/java/xdsopl/robot36/MainActivity.java b/app/src/main/java/xdsopl/robot36/MainActivity.java index b412578..96608e1 100644 --- a/app/src/main/java/xdsopl/robot36/MainActivity.java +++ b/app/src/main/java/xdsopl/robot36/MainActivity.java @@ -71,6 +71,7 @@ public class MainActivity extends AppCompatActivity { private PixelBuffer peakMeterBuffer; private ImageView peakMeterView; private PixelBuffer imageBuffer; + private ShortTimeFourierTransform stft; private short[] shortBuffer; private float[] recordBuffer; private AudioRecord audioRecord; @@ -78,6 +79,7 @@ public class MainActivity extends AppCompatActivity { private Menu menu; private String currentMode; private String language; + private Complex input; private int recordRate; private int recordChannel; private int audioSource; @@ -137,8 +139,9 @@ public class MainActivity extends AppCompatActivity { recordBuffer[i] = .000030517578125f * shortBuffer[i]; } processPeakMeter(); + processSpectrogram(); boolean newLines = decoder.process(recordBuffer, recordChannel); - processFreqPlot(); + //processFreqPlot(); if (newLines) { processScope(); processImage(); @@ -161,6 +164,77 @@ public class MainActivity extends AppCompatActivity { peakMeterView.invalidate(); } + private double clamp(double x) { + return Math.min(Math.max(x, 0), 1); + } + + private int argb(double a, double r, double g, double b) { + a = clamp(a); + r = clamp(r); + g = clamp(g); + b = clamp(b); + r *= a; + g *= a; + b *= a; + r = Math.sqrt(r); + g = Math.sqrt(g); + b = Math.sqrt(b); + int A = (int) Math.rint(255 * a); + int R = (int) Math.rint(255 * r); + int G = (int) Math.rint(255 * g); + int B = (int) Math.rint(255 * b); + return (A << 24) | (R << 16) | (G << 8) | B; + } + + private int rainbow(double v) { + v = clamp(v); + double t = 4 * v - 2; + return argb(4 * v, t, 1 - Math.abs(t), -t); + } + + private void processSpectrogram() { + boolean process = false; + int channels = recordChannel > 0 ? 2 : 1; + for (int j = 0; j < recordBuffer.length / channels; ++j) { + switch (recordChannel) { + case 1: + input.set(recordBuffer[2 * j]); + break; + case 2: + input.set(recordBuffer[2 * j + 1]); + break; + case 3: + input.set(recordBuffer[2 * j] + recordBuffer[2 * j + 1]); + break; + case 4: + input.set(recordBuffer[2 * j], recordBuffer[2 * j + 1]); + break; + default: + input.set(recordBuffer[j]); + } + if (stft.push(input)) { + process = true; + int stride = freqPlotBuffer.width; + int line = stride * freqPlotBuffer.line; + double lowest = Math.log(0.000001); + double highest = Math.log(1); + double range = highest - lowest; + for (int i = 0; i < stride; ++i) + freqPlotBuffer.pixels[line + i] = rainbow((Math.log(stft.power[i]) - lowest) / range); + System.arraycopy(freqPlotBuffer.pixels, line, freqPlotBuffer.pixels, line + stride * (freqPlotBuffer.height / 2), stride); + freqPlotBuffer.line = (freqPlotBuffer.line + 1) % (freqPlotBuffer.height / 2); + } + } + if (process) { + int width = freqPlotBitmap.getWidth(); + int height = freqPlotBitmap.getHeight(); + int stride = freqPlotBuffer.width; + int offset = stride * (freqPlotBuffer.line + freqPlotBuffer.height / 2 - height); + freqPlotBitmap.setPixels(freqPlotBuffer.pixels, offset, stride, 0, 0, width, height); + freqPlotView.invalidate(); + } + } + private void processFreqPlot() { int width = freqPlotBitmap.getWidth(); int height = freqPlotBitmap.getHeight(); @@ -238,6 +312,7 @@ public class MainActivity extends AppCompatActivity { if (rateChanged) { decoder = new Decoder(scopeBuffer, imageBuffer, getString(R.string.raw_mode), recordRate); decoder.setMode(currentMode); + stft = new ShortTimeFourierTransform(recordRate / 10, 3); } startListening(); } else { @@ -458,6 +533,7 @@ public class MainActivity extends AppCompatActivity { freqPlotBuffer = new PixelBuffer(256, 2 * 256); peakMeterBuffer = new PixelBuffer(1, 16); imageBuffer = new PixelBuffer(800, 616); + input = new Complex(); createScope(config); createFreqPlot(config); createPeakMeter(); diff --git a/app/src/main/java/xdsopl/robot36/ShortTimeFourierTransform.java b/app/src/main/java/xdsopl/robot36/ShortTimeFourierTransform.java new file mode 100644 index 0000000..7216944 --- /dev/null +++ b/app/src/main/java/xdsopl/robot36/ShortTimeFourierTransform.java @@ -0,0 +1,50 @@ +/* +Short Time Fourier Transform + +Copyright 2025 Ahmet Inan +*/ + +package xdsopl.robot36; + +public class ShortTimeFourierTransform { + private final FastFourierTransform fft; + private final Complex[] prev, fold, freq; + private final float[] weight; + private final Complex temp; + private int index; + + public final float[] power; + + ShortTimeFourierTransform(int length, int overlap) { + fft = new FastFourierTransform(length); + prev = new Complex[length * overlap]; + for (int i = 0; i < length * overlap; ++i) + prev[i] = new Complex(); + fold = new Complex[length]; + for (int i = 0; i < length; ++i) + fold[i] = new Complex(); + freq = new Complex[length]; + for (int i = 0; i < length; ++i) + freq[i] = new Complex(); + temp = new Complex(); + power = new float[length]; + weight = new float[length * overlap]; + for (int i = 0; i < length * overlap; ++i) + weight[i] = (float)(Filter.lowPass(1, length, i, length * overlap) * Hann.window(i, length * overlap)); + } + + boolean push(Complex input) { + prev[index].set(input); + index = (index + 1) % prev.length; + if (index % fold.length != 0) + return false; + for (int i = 0; i < fold.length; ++i, index = (index + 1) % prev.length) + fold[i].set(prev[index]).mul(weight[i]); + for (int i = fold.length; i < prev.length; ++i, index = (index + 1) % prev.length) + fold[i % fold.length].add(temp.set(prev[index]).mul(weight[i])); + fft.forward(freq, fold); + for (int i = 0; i < power.length; ++i) + power[i] = freq[i].norm(); + return true; + } +} \ No newline at end of file