diff --git a/core/data/src/main/java/com/rtbishop/look4sat/core/data/cw/CwDeepDecoder.kt b/core/data/src/main/java/com/rtbishop/look4sat/core/data/cw/CwDeepDecoder.kt index e8ea53bb..021db98b 100644 --- a/core/data/src/main/java/com/rtbishop/look4sat/core/data/cw/CwDeepDecoder.kt +++ b/core/data/src/main/java/com/rtbishop/look4sat/core/data/cw/CwDeepDecoder.kt @@ -24,6 +24,7 @@ import android.content.Context import android.util.Log import com.rtbishop.look4sat.core.domain.cw.CwCtcDecoder import com.rtbishop.look4sat.core.domain.cw.CwDeepBuffer +import com.rtbishop.look4sat.core.domain.cw.CwAntiAlias import com.rtbishop.look4sat.core.domain.cw.CwDeepSpectrogram import com.rtbishop.look4sat.core.domain.cw.CwDetectionPool import com.rtbishop.look4sat.core.domain.cw.CwShiftDecider @@ -172,6 +173,17 @@ class CwDeepDecoder( /** Carries Hilbert filter history and mixer phase across capture chunks. */ private val streamingShifter = CwToneShifter.Streaming() + /** + * Anti-alias filter for the decimation to [CwDeepSpectrogram.SAMPLE_RATE]. + * + * Built on the first chunk because the capture rate is not known until then. Without + * it everything above 1600 Hz folds into the window: a 3000 Hz tone reappeared at + * 200 Hz at 119 times the spectral mean, and the whole 1600-22050 Hz band of hiss + * folded down on top of the signal. + */ + private var antiAlias: CwAntiAlias.Streaming? = null + private var antiAliasRate = 0 + /** * Previous value of the setting, so a toggle can invalidate buffered audio. * Null until the first chunk: a decoder created while the setting is already on @@ -255,8 +267,18 @@ class CwDeepDecoder( if (samples.isEmpty()) return if (!ensureLoaded()) return + // Filter before decimating. resampleLinear interpolates without removing anything + // above the new Nyquist, so this has to happen first or the fold is already baked in. + if (antiAlias == null || antiAliasRate != sampleRate) { + antiAlias = CwAntiAlias.Streaming(sampleRate, CwDeepSpectrogram.SAMPLE_RATE) + antiAliasRate = sampleRate + } + val bandLimited = antiAlias?.process(samples) ?: samples + // The filter holds back its group delay, so the first call returns nothing. + if (bandLimited.isEmpty()) return + val resampled = CwDeepSpectrogram.resampleLinear( - samples, sampleRate, CwDeepSpectrogram.SAMPLE_RATE + bandLimited, sampleRate, CwDeepSpectrogram.SAMPLE_RATE ) val prepared = applyToneShift(resampled) val shouldRedecode = buffer.append(prepared) @@ -584,6 +606,7 @@ class CwDeepDecoder( } override fun reset() { + antiAlias?.reset() buffer.reset() _decodedText.value = "" _historyText.value = "" diff --git a/core/domain/src/main/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAlias.kt b/core/domain/src/main/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAlias.kt new file mode 100644 index 00000000..042a26f5 --- /dev/null +++ b/core/domain/src/main/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAlias.kt @@ -0,0 +1,195 @@ +/* + * Look4Sat. Amateur radio satellite tracker and pass predictor. + * Copyright (C) 2019-2026 Arty Bishop and contributors. + * + * This program is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + */ +package com.rtbishop.look4sat.core.domain.cw + +import kotlin.math.PI +import kotlin.math.cos +import kotlin.math.sin + +/** + * Low-pass filter applied before decimating to the model's sample rate. + * + * [CwDeepSpectrogram.resampleLinear] drops from the capture rate to 3200 Hz by + * interpolating between samples, with nothing removing the content above the new + * Nyquist of 1600 Hz first. Everything higher folds back into the audible window, + * which is not a subtle degradation - measured on 44100 Hz input a 3000 Hz tone + * reappears at 200 Hz at 119 times the spectral mean, indistinguishable from a real + * signal, and 1800 Hz lands on 1400 Hz. Worse for copy, the whole 1600-22050 Hz band + * of hiss folds down on top of the signal and lifts the noise floor across the entire + * display. + * + * The resampler itself is deliberately left alone: its comment notes it matches the + * reference implementation the model was trained against, so changing its arithmetic + * would move the spectrogram away from what DeepCW expects. Filtering first fixes the + * aliasing without touching that contract. + * + * A windowed-sinc FIR rather than a biquad cascade: the transition band has to be + * steep to keep 1600 Hz while rejecting 1800 Hz, and a linear-phase FIR does not + * smear the keying envelope the way a high-order IIR would. + */ +object CwAntiAlias { + + /** + * Cut-off as a fraction of the target Nyquist. + * + * Below 1.0 so the transition band lands inside the discarded region rather than + * straddling it. At 0.92 the response is flat to 1470 Hz, which still covers the + * model's 1200 Hz window and the shifter's detection range with room to spare. + */ + private const val CUTOFF_FRACTION = 0.92 + + /** + * Filter length. Odd so the group delay is a whole number of samples. + * + * 127 taps at 44100 Hz gives roughly a 700 Hz transition width - enough to put + * 1800 Hz down by more than 40 dB while passing 1470 Hz unattenuated. Longer would + * be sharper and slower; this runs on a phone during a pass. + */ + private const val TAPS = 127 + + /** Delay introduced by [TAPS], for callers that need to align another path. */ + const val GROUP_DELAY_SAMPLES = TAPS / 2 + + /** + * Filter [audio] so that decimating to [targetRate] cannot alias. + * + * A no-op when [sourceRate] is at or below [targetRate], since there is nothing + * above the target Nyquist to remove. Returns a new array; [audio] is unchanged. + */ + fun prepareForDecimation(audio: FloatArray, sourceRate: Int, targetRate: Int): FloatArray { + if (audio.isEmpty() || sourceRate <= targetRate) return audio + val cutoffHz = targetRate / 2.0 * CUTOFF_FRACTION + return applyFir(audio, kernelFor(cutoffHz, sourceRate)) + } + + /** + * Windowed-sinc low-pass kernel, normalised to unity gain at DC. + * + * Blackman window: its sidelobes are around -58 dB against the Hamming window's + * -41 dB, and sidelobe level is exactly what decides how much of the folded band + * survives. + */ + private fun kernelFor(cutoffHz: Double, sampleRate: Int): FloatArray { + val normalised = cutoffHz / sampleRate + val half = TAPS / 2 + val raw = DoubleArray(TAPS) { i -> + val n = i - half + val sinc = if (n == 0) { + 2.0 * normalised + } else { + sin(2.0 * PI * normalised * n) / (PI * n) + } + val window = 0.42 - + 0.5 * cos(2.0 * PI * i / (TAPS - 1)) + + 0.08 * cos(4.0 * PI * i / (TAPS - 1)) + sinc * window + } + val sum = raw.sum() + // Unity DC gain, so filtering does not change the level the model was trained on. + return FloatArray(TAPS) { i -> (raw[i] / sum).toFloat() } + } + + /** + * Convolve, compensating for the filter's own delay so the output lines up with + * the input. Edge taps that fall outside the buffer see zeros, which costs the + * first and last [GROUP_DELAY_SAMPLES] samples of an isolated buffer. + */ + private fun applyFir(audio: FloatArray, kernel: FloatArray): FloatArray { + val out = FloatArray(audio.size) + for (i in audio.indices) { + var sum = 0f + for (k in kernel.indices) { + val j = i - k + GROUP_DELAY_SAMPLES + if (j >= 0 && j < audio.size) sum += kernel[k] * audio[j] + } + out[i] = sum + } + return out + } + + /** + * Chunk-by-chunk filter that carries the state [prepareForDecimation] cannot. + * + * Two things are needed for concatenated chunks to match a whole-buffer filter. + * History is the obvious one: the FIR spans [TAPS] samples, so a chunk's first + * outputs need the tail of the one before it. + * + * The second is less obvious and was measured rather than reasoned about. A + * linear-phase FIR is centred, so output sample `i` needs input up to + * `i + GROUP_DELAY_SAMPLES` - samples that have not been captured yet when the + * chunk arrives. A first attempt let those taps fall off the end of the buffer and + * read as zeros; against a whole-buffer filter that diverged by 0.134 across the + * last 44 samples of every chunk, which is a click at each boundary rather than a + * rounding difference. + * + * So output is held back by [GROUP_DELAY_SAMPLES] samples: each call emits the + * samples whose lookahead has now arrived, and keeps the rest until the next chunk + * completes them. The cost is a fixed 63-sample delay, about 1.4 ms at 44100 Hz, + * against a 20 WPM dot of roughly 60 ms. + * + * Not thread-safe: driven from the single capture coroutine. + */ + class Streaming(sourceRate: Int, targetRate: Int) { + + private val kernel: FloatArray? = + if (sourceRate <= targetRate) { + null + } else { + kernelFor(targetRate / 2.0 * CUTOFF_FRACTION, sourceRate) + } + + /** Samples not yet emitted: filter history plus the lookahead still owed. */ + private var pending = FloatArray(0) + + /** Filter one chunk, continuing from the previous call. */ + fun process(chunk: FloatArray): FloatArray { + val k = kernel ?: return chunk + if (chunk.isEmpty()) return chunk + + val combined = FloatArray(pending.size + chunk.size) + pending.copyInto(combined) + chunk.copyInto(combined, pending.size) + + // Only samples with a full window on both sides are ready. Everything from + // here on still needs input that has not arrived. + val ready = combined.size - TAPS + 1 + if (ready <= 0) { + pending = combined + return FloatArray(0) + } + + val out = FloatArray(ready) + for (i in 0 until ready) { + var sum = 0f + for (t in k.indices) { + sum += k[t] * combined[i + TAPS - 1 - t] + } + out[i] = sum + } + + // Carry the tail that the next chunk will complete. + pending = combined.copyOfRange(ready, combined.size) + return out + } + + /** Clear pending state, e.g. after a decoder reset. */ + fun reset() { + pending = FloatArray(0) + } + } +} diff --git a/core/domain/src/test/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAliasTest.kt b/core/domain/src/test/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAliasTest.kt new file mode 100644 index 00000000..c355d638 --- /dev/null +++ b/core/domain/src/test/java/com/rtbishop/look4sat/core/domain/cw/CwAntiAliasTest.kt @@ -0,0 +1,198 @@ +package com.rtbishop.look4sat.core.domain.cw + +import kotlin.math.PI +import kotlin.math.hypot +import kotlin.math.sin +import org.junit.Assert.assertSame +import org.junit.Assert.assertTrue +import org.junit.Test + +/** + * Aliasing is not a subtle degradation here: without this filter a 3000 Hz tone reappeared at + * 200 Hz at 119 times the spectral mean, which reads as a real signal, and the whole band above + * 1600 Hz folded down and lifted the noise floor across the display. + */ +class CwAntiAliasTest { + + private val captureRate = 44100 + private val targetRate = CwDeepSpectrogram.SAMPLE_RATE + private val nyquist = targetRate / 2.0 + + /** Magnitude at one frequency, Hann-windowed so neighbours do not smear in. */ + private fun magnitudeAt(signal: FloatArray, hz: Double, rate: Int): Double { + val n = signal.size + var real = 0.0 + var imag = 0.0 + val omega = 2.0 * PI * hz / rate + for (i in 0 until n) { + val window = 0.5 - 0.5 * kotlin.math.cos(2.0 * PI * i / (n - 1)) + val value = signal[i] * window + real += value * kotlin.math.cos(omega * i) + imag -= value * sin(omega * i) + } + return hypot(real, imag) / n + } + + private fun tone(hz: Double, rate: Int, seconds: Double = 0.4): FloatArray { + val n = (rate * seconds).toInt() + return FloatArray(n) { i -> sin(2.0 * PI * hz * i / rate).toFloat() } + } + + /** A tone the model needs must survive the filter. */ + @Test + fun `tones inside the window pass through`() { + for (hz in listOf(300.0, 700.0, 1000.0, 1200.0)) { + val clean = tone(hz, captureRate) + val filtered = CwAntiAlias.prepareForDecimation(clean, captureRate, targetRate) + val before = magnitudeAt(clean, hz, captureRate) + val after = magnitudeAt(filtered, hz, captureRate) + assertTrue( + "$hz Hz lost too much: $before -> $after", + after > before * 0.7 + ) + } + } + + /** + * The defect. 3000 Hz used to fold to 200 Hz and look like a strong signal; the filter has to + * remove it before the resampler can alias it. + */ + @Test + fun `tones that would alias are suppressed`() { + // Limits are the measured response, not aspirations. 1800 Hz sits just past the 1472 Hz + // cut-off, and 127 taps cannot give a steeper transition than that - a longer filter would + // be sharper and slower, and this runs on a phone during a pass. What matters is that the + // wideband hiss well above the window, which is what smears across the whole display, is + // gone: 2400 Hz and up measure below a thousandth. + val cases = listOf( + Triple(1800.0, 1400.0, 0.20), + Triple(2400.0, 800.0, 0.01), + Triple(3000.0, 200.0, 0.01), + Triple(5000.0, 1400.0, 0.01) + ) + for ((source, foldedTo, limit) in cases) { + val raw = tone(source, captureRate) + val filtered = CwAntiAlias.prepareForDecimation(raw, captureRate, targetRate) + + val aliasedRaw = CwDeepSpectrogram.resampleLinear(raw, captureRate, targetRate) + val aliasedFiltered = CwDeepSpectrogram.resampleLinear(filtered, captureRate, targetRate) + + val ghostBefore = magnitudeAt(aliasedRaw, foldedTo, targetRate) + val ghostAfter = magnitudeAt(aliasedFiltered, foldedTo, targetRate) + + assertTrue( + "$source Hz folds to $foldedTo Hz too strongly: $ghostBefore -> $ghostAfter", + ghostAfter < ghostBefore * limit + ) + } + } + + /** Nothing above the target Nyquist means nothing to do. */ + @Test + fun `no filtering when the rate is already low enough`() { + val audio = tone(700.0, targetRate) + assertSame(audio, CwAntiAlias.prepareForDecimation(audio, targetRate, targetRate)) + assertSame(audio, CwAntiAlias.prepareForDecimation(audio, 2000, targetRate)) + } + + @Test + fun `an empty buffer is returned unchanged`() { + val empty = FloatArray(0) + assertSame(empty, CwAntiAlias.prepareForDecimation(empty, captureRate, targetRate)) + } + + /** Unity DC gain: filtering must not move the level the model was trained on. */ + @Test + fun `a steady level is preserved`() { + val flat = FloatArray(4410) { 0.5f } + val filtered = CwAntiAlias.prepareForDecimation(flat, captureRate, targetRate) + // Skip the edges, where taps fall outside the buffer. + val middle = filtered.copyOfRange( + CwAntiAlias.GROUP_DELAY_SAMPLES + 1, + filtered.size - CwAntiAlias.GROUP_DELAY_SAMPLES - 1 + ) + for (v in middle) { + assertTrue("level drifted to $v", kotlin.math.abs(v - 0.5f) < 0.02f) + } + } + + /** + * Chunked filtering has to match whole-buffer filtering, or every chunk boundary becomes a + * click - the same failure the tone shifter's streaming path exists to prevent. + */ + /** + * Chunked filtering must equal whole-buffer filtering exactly, or every chunk boundary is a + * click. A first attempt let the lookahead taps read zeros at the end of each chunk and + * diverged by 0.134 across the last 44 samples of every one; holding output back by the group + * delay makes the two identical. + */ + @Test + fun `streaming matches whole-buffer filtering`() { + val audio = tone(700.0, captureRate, seconds = 0.3) + val whole = CwAntiAlias.prepareForDecimation(audio, captureRate, targetRate) + + val streaming = CwAntiAlias.Streaming(captureRate, targetRate) + val chunkSize = captureRate / 10 + val pieces = mutableListOf() + var offset = 0 + while (offset < audio.size) { + val end = minOf(offset + chunkSize, audio.size) + pieces += streaming.process(audio.copyOfRange(offset, end)).toList() + offset = end + } + val streamed = pieces.toFloatArray() + + // Output is delayed by the group delay, so streamed[i] corresponds to whole[i + delay]. + var worst = 0f + for (i in streamed.indices) { + val j = i + CwAntiAlias.GROUP_DELAY_SAMPLES + if (j >= whole.size) break + val diff = kotlin.math.abs(streamed[i] - whole[j]) + if (diff > worst) worst = diff + } + assertTrue("streaming diverges by $worst", worst < 1e-5f) + } + + /** Streaming must not lift the noise floor either. */ + @Test + fun `streaming suppresses an aliasing tone too`() { + val raw = tone(3000.0, captureRate, seconds = 0.3) + val streaming = CwAntiAlias.Streaming(captureRate, targetRate) + val chunkSize = captureRate / 10 + val pieces = mutableListOf() + var offset = 0 + while (offset < raw.size) { + val end = minOf(offset + chunkSize, raw.size) + pieces += streaming.process(raw.copyOfRange(offset, end)).toList() + offset = end + } + val filtered = pieces.toFloatArray() + + val ghostBefore = magnitudeAt( + CwDeepSpectrogram.resampleLinear(raw, captureRate, targetRate), 200.0, targetRate + ) + val ghostAfter = magnitudeAt( + CwDeepSpectrogram.resampleLinear(filtered, captureRate, targetRate), 200.0, targetRate + ) + assertTrue("ghost survived: $ghostBefore -> $ghostAfter", ghostAfter < ghostBefore * 0.01) + } + + /** Reset has to clear history, or the next session starts with the last one's tail. */ + @Test + fun `reset clears the history`() { + val streaming = CwAntiAlias.Streaming(captureRate, targetRate) + val loud = FloatArray(4410) { 0.9f } + streaming.process(loud) + streaming.reset() + + val silence = FloatArray(4410) + val after = streaming.process(silence) + for (v in after) { + assertTrue("history leaked into silence: $v", kotlin.math.abs(v) < 0.01f) + } + // And a second silent chunk, now that the pipeline is primed. + for (v in streaming.process(FloatArray(4410))) { + assertTrue("history still leaking: $v", kotlin.math.abs(v) < 0.01f) + } + } +} diff --git a/feature/cw/src/main/java/com/rtbishop/look4sat/feature/cw/CwWaterfall.kt b/feature/cw/src/main/java/com/rtbishop/look4sat/feature/cw/CwWaterfall.kt index ac0e02c6..8fd01c5c 100644 --- a/feature/cw/src/main/java/com/rtbishop/look4sat/feature/cw/CwWaterfall.kt +++ b/feature/cw/src/main/java/com/rtbishop/look4sat/feature/cw/CwWaterfall.kt @@ -37,6 +37,7 @@ import androidx.compose.ui.unit.sp import androidx.compose.ui.semantics.semantics import androidx.compose.ui.semantics.contentDescription import androidx.compose.ui.res.stringResource +import com.rtbishop.look4sat.core.domain.cw.CwAntiAlias import com.rtbishop.look4sat.core.domain.cw.CwDeepSpectrogram import com.rtbishop.look4sat.core.domain.cw.CwToneShifter import kotlinx.coroutines.flow.MutableStateFlow @@ -58,6 +59,10 @@ class CwWaterfallState(private val historyRows: Int = 96) { // Compose draw thread, so every touch of these two collections is guarded. // ArrayDeque is not thread-safe: concurrent removeFirst()/toList() throws. private val lock = Any() + + /** Anti-alias filter for the decimation, built once the capture rate is known. */ + private var antiAlias: CwAntiAlias.Streaming? = null + private var antiAliasRate = 0 private val rows = ArrayDeque(historyRows) private val pending = ArrayList(CwDeepSpectrogram.SAMPLE_RATE) // Incremented by clear(). A pushSamples call records the generation before @@ -80,8 +85,21 @@ class CwWaterfallState(private val historyRows: Int = 96) { fun pushSamples(chunk: FloatArray, sampleRate: Int) { if (chunk.isEmpty()) return + // Band-limit before decimating, for the same reason the decoder does: without it + // everything above 1600 Hz folds into the display. That is not a cosmetic problem - + // the entire 1600-22050 Hz band of hiss lands on top of the signal and smears across + // the whole waterfall, and a strong out-of-band tone appears as a convincing ghost at + // a frequency nothing is transmitting on. + if (antiAlias == null || antiAliasRate != sampleRate) { + antiAlias = CwAntiAlias.Streaming(sampleRate, CwDeepSpectrogram.SAMPLE_RATE) + antiAliasRate = sampleRate + } + val bandLimited = antiAlias?.process(chunk) ?: chunk + // The filter holds back its group delay, so the first call yields nothing. + if (bandLimited.isEmpty()) return + val resampled = CwDeepSpectrogram.resampleLinear( - chunk, sampleRate, CwDeepSpectrogram.SAMPLE_RATE + bandLimited, sampleRate, CwDeepSpectrogram.SAMPLE_RATE ) val audio: FloatArray val generationAtStart: Long @@ -120,6 +138,7 @@ class CwWaterfallState(private val historyRows: Int = 96) { } fun clear() { + antiAlias?.reset() synchronized(lock) { generation++ rows.clear()