feat(cw): v3 spectrogram-based multi-channel Bayesian decoder

Complete rewrite inspired by Morse Expert / CW Skimmer (VE3NEA):
- CwFFT: radix-2 FFT (256-point) for time-frequency analysis
- CwSpectrogram: sliding-window waterfall (40 cols x 33 bins, 8ms resolution)
- CwBayesianDecoder: Gaussian probability replaces hard dit/dash thresholds
- CwChannelTracker: multi-channel peak detection (up to 3 signals)
- CwDecoder: integrates all components, monitors 200-1200 Hz simultaneously

Key advantages over v2 (ggmorse):
- Frequency-agnostic: full spectrum monitored, not locked to one tone
- Multi-channel: tracks multiple signals in parallel
- Bayesian: probability-based decisions, not hard ratios
- Doppler tolerant: frequency drift just moves energy between bins
This commit is contained in:
atsunatsu committed 2026-08-01 19:41:32 +08:00
1 parent b9e5ff70f5
commit 8af80866bc
7 files changed
+1685 -408

No files matched your search

@@ -0,0 +1,154 @@
/*
* 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 <https://www.gnu.org/licenses/>.
*/
package com.rtbishop.look4sat.core.domain.cw
/**
* Bayesian Morse timing decoder.
* Replaces hard thresholds with probability-based decision making.
*
* Inspired by VE3NEA's CW Skimmer approach:
* "Instead of making a hard decision at every input sample whether the signal
* is present or not, compute the probability that the signal is present."
*
* Uses Gaussian probability density centered on expected durations:
* P(dit | duration) = exp(-(duration - dotMs)^2 / (2 * variance^2))
* P(dash | duration) = exp(-(duration - 3*dotMs)^2 / (2 * variance^2))
*/
internal class CwBayesianDecoder {
// Morse timing parameters
private var dotDurationMs = 60f // initial 20 WPM
private var speedWpm = 20f
// Current symbol being accumulated
private var currentSymbol = StringBuilder()
private var textBuffer = StringBuilder()
// Recent dit lengths for speed estimation
private val recentDits = mutableListOf<Float>()
// Output
private var _decodedText = ""
val decodedText: String get() = _decodedText
/** Gaussian probability. */
private fun gaussianProb(durationMs: Float, expectedMs: Float, varianceMs: Float): Float {
if (varianceMs <= 0f) return 0f
val diff = durationMs - expectedMs
return kotlin.math.exp(-(diff * diff) / (2 * varianceMs * varianceMs))
}
/** Process a tone duration. Returns the symbol type with highest probability. */
fun processTone(durationMs: Float): ToneResult {
val ditProb = gaussianProb(durationMs, dotDurationMs, dotDurationMs * 0.4f)
val dashProb = gaussianProb(durationMs, dotDurationMs * 3f, dotDurationMs * 0.6f)
return if (ditProb > dashProb && ditProb > 0.05f) {
currentSymbol.append('0')
recentDits.add(durationMs)
updateSpeed()
ToneResult('0', ditProb)
} else if (dashProb > 0.05f) {
currentSymbol.append('1')
ToneResult('1', dashProb)
} else {
ToneResult(null, 0f)
}
}
/** Process a gap duration. Returns decoded character or null. */
fun processGap(durationMs: Float): Char? {
if (currentSymbol.isEmpty()) {
val wordProb = gaussianProb(durationMs, dotDurationMs * 7f, dotDurationMs * 1.2f)
if (wordProb > 0.2f) {
textBuffer.append(' ')
_decodedText = textBuffer.toString()
return ' '
}
return null
}
val interCharProb = gaussianProb(durationMs, dotDurationMs * 3f, dotDurationMs * 0.6f)
val wordProb = gaussianProb(durationMs, dotDurationMs * 7f, dotDurationMs * 1.2f)
if (wordProb > interCharProb && wordProb > 0.2f) {
val char = flushSymbol()
textBuffer.append(' ')
_decodedText = textBuffer.toString()
return char
}
if (interCharProb > 0.15f) {
val char = flushSymbol()
_decodedText = textBuffer.toString()
return char
}
return null
}
private fun flushSymbol(): Char? {
if (currentSymbol.isEmpty()) return null
val morse = currentSymbol.toString()
currentSymbol.clear()
val char = morseToChar(morse)
if (char != null) textBuffer.append(char)
return char
}
private fun updateSpeed() {
if (recentDits.size < 3) return
val sorted = recentDits.sorted()
val median = sorted[sorted.size / 2]
if (median > 0f) {
dotDurationMs = dotDurationMs * 0.7f + median * 0.3f
val wpm = 60.0f / (50.0f * dotDurationMs / 1000.0f)
if (wpm in 5f..55f) speedWpm = wpm
}
}
fun getSpeed(): Float = speedWpm
fun reset() {
dotDurationMs = 60f
speedWpm = 20f
recentDits.clear()
currentSymbol.clear()
textBuffer.clear()
_decodedText = ""
}
companion object {
private val MORSE_TABLE = mapOf(
"01" to 'A', "1000" to 'B', "1010" to 'C', "100" to 'D', "0" to 'E',
"0010" to 'F', "110" to 'G', "0000" to 'H', "00" to 'I', "0111" to 'J',
"101" to 'K', "0100" to 'L', "11" to 'M', "10" to 'N', "111" to 'O',
"0110" to 'P', "1101" to 'Q', "010" to 'R', "000" to 'S', "1" to 'T',
"001" to 'U', "0001" to 'V', "011" to 'W', "1001" to 'X', "1011" to 'Y',
"1100" to 'Z', "01111" to '1', "00111" to '2', "00011" to '3',
"00001" to '4', "00000" to '5', "10000" to '6', "11000" to '7',
"11100" to '8', "11110" to '9', "11111" to '0',
"010101" to '.', "110011" to ',', "001100" to '?', "011110" to '\'',
"101011" to '!', "10010" to '/', "10110" to '(', "101101" to ')',
"01000" to '&', "111000" to ':', "101010" to ';', "10001" to '=',
"01010" to '+', "100001" to '-', "001101" to '_', "010010" to '"',
"0001001" to '$', "011010" to '@'
)
fun morseToChar(morse: String): Char? = MORSE_TABLE[morse]
}
}
data class ToneResult(val symbol: Char?, val probability: Float)
@@ -0,0 +1,114 @@
/*
* 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 <https://www.gnu.org/licenses/>.
*/
package com.rtbishop.look4sat.core.domain.cw
/**
* Multi-channel CW signal tracker.
* Monitors the spectrogram for active frequency bins and extracts
* energy envelopes for each detected signal.
*
* Inspired by CW Skimmer's multi-channel approach:
* tracks all active signals in the passband simultaneously,
* selects the best one for decoded output.
*/
internal class CwChannelTracker(
private val spectrogram: CwSpectrogram,
private val maxChannels: Int = 3
) {
data class Channel(
val bin: Int,
val frequency: Float,
var active: Boolean = false,
var energy: Float = 0f,
val history: MutableList<Float> = mutableListOf(),
var confidence: Float = 0f
)
private val channels = Array(maxChannels) { Channel(0, 0f) }
/** Scan the current spectrogram column and update channel tracking. */
fun update(): List<Channel> {
val col = spectrogram.getCurrentColumn()
val peaks = findPeaks(col, threshold = 0.3f, minDistance = 2)
// Update existing channels
for (ch in channels) {
if (ch.active) {
if (peaks.contains(ch.bin)) {
ch.energy = col[ch.bin]
ch.history.add(ch.energy)
if (ch.history.size > 40) ch.history.removeAt(0)
ch.confidence = computeConfidence(ch.history)
} else {
// Signal lost — decay confidence
ch.history.add(0f)
if (ch.history.size > 40) ch.history.removeAt(0)
ch.confidence *= 0.9f
if (ch.confidence < 0.1f) ch.active = false
}
}
}
// Assign new peaks to inactive channels
var peakIdx = 0
for (ch in channels) {
if (!ch.active && peakIdx < peaks.size) {
val bin = peaks[peakIdx]
val freq = spectrogram.binToFreq(bin)
// Re-initialize channel
channels[peakIdx] = Channel(bin, freq, true, col[bin], mutableListOf(), 0.5f)
peakIdx++
}
}
return channels.filter { it.active }
}
/** Find peak bins in the spectrum. */
private fun findPeaks(spectrum: FloatArray, threshold: Float, minDistance: Int): List<Int> {
val peaks = mutableListOf<Int>()
for (i in 1 until spectrum.size - 1) {
if (spectrum[i] > spectrum[i - 1] && spectrum[i] > spectrum[i + 1] && spectrum[i] > threshold) {
if (peaks.isEmpty() || i - peaks.last() >= minDistance) {
peaks.add(i)
}
}
}
return peaks.sortedByDescending { spectrum[it] }
}
/** Compute confidence from energy history. Lower variance = higher confidence. */
private fun computeConfidence(history: List<Float>): Float {
if (history.size < 10) return 0.3f
val recent = history.takeLast(10)
val mean = recent.average().toFloat()
val variance = recent.map { (it - mean) * (it - mean) }.average().toFloat()
return if (mean > 0f) (mean / (mean + variance + 0.1f)).coerceIn(0f, 1f) else 0f
}
/** Get the channel with highest confidence. */
fun getBestChannel(): Channel? {
return channels.filter { it.active }.maxByOrNull { it.confidence }
}
fun reset() {
for (i in channels.indices) {
channels[i] = Channel(0, 0f)
}
}
}
@@ -19,52 +19,50 @@ package com.rtbishop.look4sat.core.domain.cw
import kotlinx.coroutines.flow.MutableStateFlow
import kotlinx.coroutines.flow.StateFlow
import kotlin.math.abs
/**
* CW (Morse code) decoder ported from ggerganov/ggmorse.
* CW (Morse code) decoder v3 — Spectrogram-based multi-channel Bayesian decoder.
*
* Key improvements over v1:
* - Automatic pitch detection (200-1200 Hz) via DFT
* - Automatic speed detection (5-55 WPM) via interval clustering
* - Adaptive threshold with signal statistics
* - Resampling to 4 kHz base rate for efficiency
* - Running Goertzel filter for tone detection
* - First-order IIR bandpass filter (HP + LP)
* Architecture inspired by Morse Expert / CW Skimmer (VE3NEA):
* 1. FFT spectrogram creates a frequency×time matrix
* 2. Multi-channel peak detector finds all active signals
* 3. Per-channel energy envelope extraction
* 4. Bayesian probability for symbol timing (Gaussian likelihood)
* 5. Best channel selected for output
*
* Algorithm flow:
* Audio buffer → Resample to 4 kHz → High-pass filter (200 Hz) →
* Low-pass filter (1200 Hz) → Pitch detection (DFT, 200-1200 Hz) →
* Running Goertzel at detected pitch → Adaptive threshold →
* Signal interval timing → Speed estimation → Morse character lookup
* Key advantages over v2 (ggmorse):
* - Frequency-agnostic: monitors all 200-1200 Hz simultaneously
* - Multi-channel: tracks up to 3 signals in parallel
* - Bayesian: probability-based decisions, not hard thresholds
* - Frequency drift tolerant: energy just moves between bins
*/
class CwDecoder(
val sampleRate: Int = 8000,
cwToneFreq: Float = -1f, // -1 = auto-detect
minFreq: Float = 200f,
maxFreq: Float = 1200f
cwToneFreq: Float = -1f // ignored in v3 (auto-detect via spectrogram)
) {
companion object {
private val MORSE_TABLE = mapOf(
"01" to 'A', "1000" to 'B', "1010" to 'C', "100" to 'D', "0" to 'E',
"0010" to 'F', "110" to 'G', "0000" to 'H', "00" to 'I', "0111" to 'J',
"101" to 'K', "0100" to 'L', "11" to 'M', "10" to 'N', "111" to 'O',
"0110" to 'P', "1101" to 'Q', "010" to 'R', "000" to 'S', "1" to 'T',
"001" to 'U', "0001" to 'V', "011" to 'W', "1001" to 'X', "1011" to 'Y',
"1100" to 'Z', "01111" to '1', "00111" to '2', "00011" to '3',
"00001" to '4', "00000" to '5', "10000" to '6', "11000" to '7',
"11100" to '8', "11110" to '9', "11111" to '0',
"010101" to '.', "110011" to ',', "001100" to '?', "011110" to '\'',
"101011" to '!', "10010" to '/', "10110" to '(', "101101" to ')',
"01000" to '&', "111000" to ':', "101010" to ';', "10001" to '=',
"01010" to '+', "100001" to '-', "001101" to '_', "010010" to '"',
"0001001" to '$', "011010" to '@'
)
fun morseToChar(morse: String): Char? = MORSE_TABLE[morse]
private val spectrogram = CwSpectrogram(
fftSize = 256,
hopSize = 64,
sampleRate = sampleRate,
minBin = 6,
maxBin = 38,
historyCols = 40
)
private val channelTracker = CwChannelTracker(spectrogram, maxChannels = 3)
private val bayesianDecoder = CwBayesianDecoder()
private const val BASE_SAMPLE_RATE = 4000f
private const val PITCH_DETECT_INTERVAL = 100 // frames between pitch scans
}
// Per-channel state
private data class ChannelTiming(
var isSignal: Boolean = false,
var toneSamples: Int = 0,
var gapSamples: Int = 0
)
private val timingStates = Array(3) { ChannelTiming() }
// Sample period in milliseconds (at 4 kHz effective rate for timing)
// The spectrogram processes at native sample rate, but timing analysis
// uses the spectrogram column rate: hopSize/sampleRate seconds per column
private val samplePeriodMs = 1000f * hopSize / sampleRate
// Output flows
private val _decodedTextFlow = MutableStateFlow("")
@@ -79,223 +77,85 @@ class CwDecoder(
private val _estimatedSpeed = MutableStateFlow<Float?>(null)
val estimatedSpeed: StateFlow<Float?> = _estimatedSpeed
// DSP components
private val resampler = CwResampler(sampleRate.toFloat(), BASE_SAMPLE_RATE)
private val hpFilter = CwFilter()
private val lpFilter = CwFilter()
private val pitchDetector = CwPitchDetector(BASE_SAMPLE_RATE, minFreq, maxFreq)
private val goertzel = CwGoertzel()
// Decoder state
private var decodedText = StringBuilder()
private var currentLetter = StringBuilder()
private var isSignal = false
private var signalOnSamples = 0
private var signalOffSamples = 0
private var pitchEstimate = if (cwToneFreq > 0f) cwToneFreq else -1f
private var pitchConfidenceCounter = 0
private var noiseFloor = 0.0f
private var signalPeak = 0.0f
private var speedEstimate = 20f // initial guess: 20 WPM
private var isPitchLocked = cwToneFreq > 0f
// Frame counter for periodic updates
private var frameCount = 0
init {
if (isPitchLocked) {
if (cwToneFreq > 0f) {
_estimatedPitch.value = cwToneFreq
}
}
// Interval history for speed estimation
private val intervalHistory = mutableListOf<Int>() // lengths of dits (type 0 only)
// Keep track of last processed sample for the goertzel filter
private var goertzelSampleCount = 0
companion object {
private const val hopSize = 64
}
fun processBuffer(buffer: FloatArray) {
// 1. Resample to 4 kHz base rate
val resampled = resampler.process(buffer)
// 1. Feed samples to spectrogram (generates FFT waterfall)
spectrogram.addSamples(buffer)
for (sample in resampled) {
// 2. Bandpass filter chain: 200 Hz HP → 1200 Hz LP
val hp = hpFilter.highPass(sample, 200f, BASE_SAMPLE_RATE)
val filtered = lpFilter.lowPass(hp, 1200f, BASE_SAMPLE_RATE)
val absVal = abs(filtered)
// 2. Update channel tracker (find active frequency bins)
val activeChannels = channelTracker.update()
// 3. Update noise floor and signal peak (running statistics)
noiseFloor = 0.999f * noiseFloor + 0.001f * absVal
if (absVal > signalPeak) {
signalPeak = absVal
} else {
signalPeak = 0.999f * signalPeak
}
// 3. For each active channel, extract timing
for ((idx, channel) in activeChannels.withIndex()) {
if (idx >= timingStates.size) break
val state = timingStates[idx]
val col = spectrogram.getCurrentColumn()
val energy = if (channel.bin in col.indices) col[channel.bin] else 0f
// 4. Adaptive threshold
val threshold = (noiseFloor + (signalPeak - noiseFloor) * 0.3f)
_signalStrength.value = if (signalPeak > 0f && threshold > 0f) {
((signalPeak - threshold) / signalPeak).coerceIn(0f, 1f)
} else {
0f
}
// Adaptive threshold: 30% above noise floor
val threshold = 0.3f + (energy - 0.3f) * 0.3f
// 5. Run Goertzel filter if pitch is locked
if (isPitchLocked && pitchEstimate > 0f) {
goertzel.process(filtered)
goertzelSampleCount++
}
// 6. Signal detection with adaptive threshold
if (absVal > threshold) {
if (!isSignal) {
// Rising edge — process the silence gap that just ended
if (signalOffSamples > 0) {
processGap(signalOffSamples)
if (energy > threshold) {
if (!state.isSignal) {
// Rising edge — process gap
if (state.gapSamples > 0) {
val gapMs = state.gapSamples * samplePeriodMs
bayesianDecoder.processGap(gapMs)
}
signalOffSamples = 0
isSignal = true
state.gapSamples = 0
state.isSignal = true
}
signalOnSamples++
state.toneSamples++
} else {
if (isSignal) {
// Falling edge — process the tone that just ended
processTone(signalOnSamples)
signalOnSamples = 0
isSignal = false
if (state.isSignal) {
// Falling edge — process tone
val toneMs = state.toneSamples * samplePeriodMs
bayesianDecoder.processTone(toneMs)
state.toneSamples = 0
state.isSignal = false
}
signalOffSamples++
state.gapSamples++
}
}
// 7. Periodic pitch detection (every ~100 frames)
if (!isPitchLocked) {
pitchConfidenceCounter++
if (pitchConfidenceCounter >= PITCH_DETECT_INTERVAL) {
pitchConfidenceCounter = 0
val pitch = pitchDetector.findPitch(resampled)
if (pitch != null) {
pitchEstimate = pitch
_estimatedPitch.value = pitch
isPitchLocked = true
goertzel.init(BASE_SAMPLE_RATE, pitch)
}
}
} else {
// Continuous pitch tracking: re-check periodically to handle Doppler drift
pitchConfidenceCounter++
if (pitchConfidenceCounter >= PITCH_DETECT_INTERVAL * 5) {
pitchConfidenceCounter = 0
// Narrow scan: ±100 Hz around current pitch estimate
val narrowDetector = CwPitchDetector(
BASE_SAMPLE_RATE,
(pitchEstimate - 100f).coerceAtLeast(200f),
(pitchEstimate + 100f).coerceAtMost(1200f),
5f
)
val pitch = narrowDetector.findPitch(resampled)
if (pitch != null && kotlin.math.abs(pitch - pitchEstimate) > 20f) {
pitchEstimate = pitch
_estimatedPitch.value = pitch
goertzel.init(BASE_SAMPLE_RATE, pitch)
}
}
}
// Push latest decoded text
_decodedTextFlow.value = decodedText.toString()
}
private fun processTone(samples: Int) {
val dotDuration = samplesForDot()
if (dotDuration <= 0) return
val ratio = samples.toFloat() / dotDuration
if (ratio < 1.5f) {
currentLetter.append('0') // 0 = dot
// Track dit lengths for speed estimation
intervalHistory.add(samples)
if (intervalHistory.size > 20) intervalHistory.removeAt(0)
} else if (ratio < 5.0f) {
currentLetter.append('1') // 1 = dash
}
// else: ignore very long tones (likely noise/interference)
// Update speed estimate from recent dits
updateSpeedEstimate()
}
private fun processGap(samples: Int) {
val dotDuration = samplesForDot()
if (dotDuration <= 0) return
val gapRatio = samples.toFloat() / dotDuration
if (currentLetter.isNotEmpty()) {
// Inter-character gap (3+ dot durations)
if (gapRatio >= 2.5f) {
val char = morseToChar(currentLetter.toString())
if (char != null) {
decodedText.append(char)
}
currentLetter.clear()
// Word gap (7+ dot durations)
if (gapRatio >= 7f) {
decodedText.append(' ')
}
}
} else {
// Word gap (7+ dot durations, no letter in progress)
if (gapRatio >= 7f) {
decodedText.append(' ')
}
}
}
private fun samplesForDot(): Int {
// Convert WPM to samples at 4 kHz base rate
// Using standard formula: dot = 60/(50*WPM) seconds
return ((BASE_SAMPLE_RATE * 60.0 / (50.0 * speedEstimate)).toInt()).coerceAtLeast(1)
}
private fun updateSpeedEstimate() {
if (intervalHistory.size < 3) return
// Use median of recent dit lengths for speed estimation
val sorted = intervalHistory.sorted()
val median = sorted[sorted.size / 2].toFloat()
if (median > 0f) {
val newSpeed = 60.0f / (50.0f * median / BASE_SAMPLE_RATE)
if (newSpeed in 5f..55f) {
// Smooth speed update (70% old, 30% new)
speedEstimate = speedEstimate * 0.7f + newSpeed * 0.3f
_estimatedSpeed.value = speedEstimate
// 4. Update outputs from best channel
frameCount++
if (frameCount % 5 == 0) { // Every 5 frames
val bestChannel = channelTracker.getBestChannel()
if (bestChannel != null) {
_estimatedPitch.value = bestChannel.frequency
_signalStrength.value = bestChannel.confidence
_estimatedSpeed.value = bayesianDecoder.getSpeed()
}
_decodedTextFlow.value = bayesianDecoder.decodedText
}
}
fun resetDecoder() {
isSignal = false
signalOnSamples = 0
signalOffSamples = 0
decodedText.clear()
currentLetter.clear()
intervalHistory.clear()
noiseFloor = 0.0f
signalPeak = 0.0f
speedEstimate = 20f
if (!isPitchLocked) {
pitchEstimate = -1f
pitchConfidenceCounter = 0
spectrogram.reset()
channelTracker.reset()
bayesianDecoder.reset()
for (state in timingStates) {
state.isSignal = false
state.toneSamples = 0
state.gapSamples = 0
}
goertzelSampleCount = 0
hpFilter.reset()
lpFilter.reset()
resampler.reset()
goertzel.reset()
frameCount = 0
_decodedTextFlow.value = ""
_signalStrength.value = 0f
_estimatedPitch.value = if (isPitchLocked) pitchEstimate else null
_estimatedPitch.value = null
_estimatedSpeed.value = null
}
}
@@ -0,0 +1,87 @@
/*
* 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 <https://www.gnu.org/licenses/>.
*/
package com.rtbishop.look4sat.core.domain.cw
import kotlin.math.cos
import kotlin.math.sqrt
/**
* Radix-2 FFT for real-valued input.
* Produces magnitude spectrum for the first N/2+1 bins.
* Used by CwSpectrogram for time-frequency analysis.
*/
internal class CwFFT(private val n: Int) {
init {
require(n > 0 && n and (n - 1) == 0) { "FFT size must be power of 2, got $n" }
}
private val cosTable = FloatArray(n / 2)
private val sinTable = FloatArray(n / 2)
init {
for (i in 0 until n / 2) {
val angle = -2.0 * kotlin.math.PI * i / n
cosTable[i] = cos(angle).toFloat()
sinTable[i] = kotlin.math.sin(angle).toFloat()
}
}
/** Compute magnitude spectrum for real input. Returns array of size n/2+1. */
fun magnitudeSpectrum(input: FloatArray): FloatArray {
require(input.size == n) { "Input size must be $n, got ${input.size}" }
val real = input.copyOf()
val imag = FloatArray(n)
// Bit-reversal permutation
var j = 0
for (i in 1 until n) {
var bit = n shr 1
while (j and bit != 0) { j = j xor bit; bit = bit shr 1 }
j = j xor bit
if (i < j) {
var tmp = real[i]; real[i] = real[j]; real[j] = tmp
}
}
// Radix-2 Cooley-Tukey FFT
var len = 2
while (len <= n) {
val half = len / 2
val step = n / len
for (i in 0 until n step len) {
for (k in 0 until half) {
val tReal = real[i + k + half] * cosTable[k * step] - imag[i + k + half] * sinTable[k * step]
val tImag = real[i + k + half] * sinTable[k * step] + imag[i + k + half] * cosTable[k * step]
real[i + k + half] = real[i + k] - tReal
imag[i + k + half] = imag[i + k] - tImag
real[i + k] += tReal
imag[i + k] += tImag
}
}
len = len shl 1
}
// Magnitude spectrum (first N/2+1 bins)
val mag = FloatArray(n / 2 + 1)
for (i in 0..n / 2) {
mag[i] = sqrt(real[i] * real[i] + imag[i] * imag[i]) / n
}
return mag
}
}
@@ -0,0 +1,152 @@
/*
* 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 <https://www.gnu.org/licenses/>.
*/
package com.rtbishop.look4sat.core.domain.cw
/**
* Sliding-window spectrogram for CW decoding.
* Maintains a time-frequency matrix updated with each audio frame.
*
* FFT size: 256, hop size: 64, sample rate: 4000 (or native)
* Frequency bins: 6..38 (187-1187 Hz, covers typical CW range)
* History: 40 columns (320 ms window)
* Time resolution: 64/4000 = 16 ms, Frequency resolution: 4000/256 = 15.625 Hz
*/
internal class CwSpectrogram(
private val fftSize: Int = 256,
private val hopSize: Int = 64,
private val sampleRate: Int = 4000,
private val minBin: Int = 6,
private val maxBin: Int = 38,
private val historyCols: Int = 40
) {
private val fft = CwFFT(fftSize)
val numBins: Int get() = maxBin - minBin + 1
// Hanning window
private val hanning = FloatArray(fftSize) {
(0.5 - 0.5 * kotlin.math.cos(2.0 * kotlin.math.PI * it / (fftSize - 1))).toFloat()
}
// Spectrogram data: [timeCol][freqBin]
private val spectrogram = Array(historyCols) { FloatArray(numBins) }
private var currentCol = 0
private var samplesBuffered = 0
private val buffer = FloatArray(fftSize)
// Per-bin running energy for normalization
private val binEnergy = FloatArray(numBins) { 1f }
private val alpha = 0.95f
/** Add audio samples, compute FFTs for each complete hop. */
fun addSamples(samples: FloatArray) {
var offset = 0
while (offset < samples.size) {
val needed = fftSize - samplesBuffered
val copyLen = minOf(needed, samples.size - offset)
System.arraycopy(samples, offset, buffer, samplesBuffered, copyLen)
samplesBuffered += copyLen
offset += copyLen
if (samplesBuffered >= fftSize) {
processFrame()
// Shift buffer: keep last (fftSize - hopSize) samples
System.arraycopy(buffer, hopSize, buffer, 0, fftSize - hopSize)
samplesBuffered = fftSize - hopSize
}
}
}
private fun processFrame() {
// Apply Hanning window
val windowed = FloatArray(fftSize) { buffer[it] * hanning[it] }
// Compute FFT magnitude spectrum
val mag = fft.magnitudeSpectrum(windowed)
// Update spectrogram column
val col = spectrogram[currentCol]
for (b in 0 until numBins) {
val binIdx = minBin + b
val rawMag = mag[binIdx]
// Running energy normalization
binEnergy[b] = alpha * binEnergy[b] + (1 - alpha) * rawMag
col[b] = if (binEnergy[b] > 1e-6f) rawMag / binEnergy[b] else 0f
}
currentCol = (currentCol + 1) % historyCols
}
/** Get the current spectrogram as a 2D array in chronological order. */
fun getSpectrogram(): Array<FloatArray> {
val result = Array(historyCols) { i ->
val srcIdx = (currentCol + i) % historyCols
spectrogram[srcIdx].copyOf()
}
return result
}
/** Get the most recent column (current energy across all frequencies). */
fun getCurrentColumn(): FloatArray {
val prevCol = (currentCol - 1 + historyCols) % historyCols
return spectrogram[prevCol].copyOf()
}
/** Find the frequency bin with peak energy. Returns -1 if no significant signal. */
fun findPeakBin(): Int {
val col = getCurrentColumn()
var maxBin = -1
var maxVal = 0f
for (i in col.indices) {
if (col[i] > maxVal) {
maxVal = col[i]
maxBin = i
}
}
return if (maxVal > 0.3f) maxBin else -1
}
/** Get energy at a specific bin over the last N columns in chronological order. */
fun getBinEnergy(bin: Int, numCols: Int): FloatArray {
val clamped = minOf(numCols, historyCols)
val result = FloatArray(clamped)
for (i in 0 until clamped) {
val colIdx = (currentCol - clamped + i + historyCols) % historyCols
result[i] = spectrogram[colIdx][bin]
}
return result
}
/** Get the bin index for a frequency in Hz. */
fun freqToBin(freqHz: Float): Int {
val bin = (freqHz * fftSize / sampleRate).toInt()
return (bin - minBin).coerceIn(0, numBins - 1)
}
/** Get the center frequency for a bin. */
fun binToFreq(bin: Int): Float {
return (minBin + bin).toFloat() * sampleRate / fftSize
}
fun reset() {
for (col in spectrogram) col.fill(0f)
currentCol = 0
samplesBuffered = 0
buffer.fill(0f)
binEnergy.fill(1f)
}
}
@@ -28,240 +28,274 @@ class CwDecoderTest {
@Test
fun morseToChar_basicLetters() {
assertEquals('A', CwDecoder.morseToChar("01"))
assertEquals('B', CwDecoder.morseToChar("1000"))
assertEquals('S', CwDecoder.morseToChar("000"))
assertEquals('O', CwDecoder.morseToChar("111"))
assertEquals('A', CwBayesianDecoder.morseToChar("01"))
assertEquals('S', CwBayesianDecoder.morseToChar("000"))
assertEquals('O', CwBayesianDecoder.morseToChar("111"))
}
@Test
fun morseToChar_numbers() {
assertEquals('1', CwDecoder.morseToChar("01111"))
assertEquals('5', CwDecoder.morseToChar("00000"))
assertEquals('0', CwDecoder.morseToChar("11111"))
assertEquals('1', CwBayesianDecoder.morseToChar("01111"))
assertEquals('0', CwBayesianDecoder.morseToChar("11111"))
}
@Test
fun morseToChar_unknown_returnsNull() {
assertNull(CwDecoder.morseToChar("......."))
assertNull(CwDecoder.morseToChar(""))
assertNull(CwDecoder.morseToChar("01-01"))
assertNull(CwBayesianDecoder.morseToChar("......."))
assertNull(CwBayesianDecoder.morseToChar(""))
}
// --- Resampler ---
// --- FFT ---
@Test
fun resampler_downsampleReducesSize() {
val resampler = CwResampler(8000f, 4000f)
val input = FloatArray(8000) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() }
val output = resampler.process(input)
assertTrue("Output size ${output.size} should be ~4000", output.size in 3800..4200)
}
@Test
fun resampler_emptyInput_returnsEmpty() {
val resampler = CwResampler(8000f, 4000f)
val output = resampler.process(FloatArray(0))
assertTrue(output.isEmpty())
}
@Test
fun resampler_sameRate_returnsSameSize() {
val resampler = CwResampler(4000f, 4000f)
val input = FloatArray(100) { it.toFloat() }
val output = resampler.process(input)
assertTrue("Output size should be ~100", output.size in 95..105)
}
@Test
fun resampler_resetClearsState() {
val resampler = CwResampler(8000f, 4000f)
val input = FloatArray(100) { 1f }
resampler.process(input)
resampler.reset()
// Should not crash
resampler.process(FloatArray(100) { 0f })
}
// --- Filter ---
@Test
fun filter_highPass_doesNotCrash() {
val filter = CwFilter()
val sampleRate = 4000f
// Test that the filter runs without crashing and produces finite values
val output = FloatArray(100) { filter.highPass(1.0f, 200f, sampleRate) }
output.forEach { assertFalse("Output should be finite: $it", it.isNaN() || it.isInfinite()) }
}
@Test
fun filter_lowPass_smoothsSignal() {
val filter = CwFilter()
val sampleRate = 4000f
// High frequency noise
val output = FloatArray(100) { filter.lowPass((sin(2.0 * PI * 1000.0 * it / sampleRate)).toFloat(), 500f, sampleRate) }
val maxVal = output.maxOrNull() ?: 1f
assertTrue("High freq should be attenuated, max=$maxVal", maxVal < 0.8f)
}
@Test
fun filter_reset() {
val filter = CwFilter()
filter.highPass(1f, 200f, 4000f)
filter.reset()
// Should not crash
assertEquals(0f, filter.highPass(0f, 200f, 4000f), 0.001f)
}
// --- Goertzel ---
@Test
fun goertzel_detectsPresentTone() {
val sampleRate = 4000f
val targetFreq = 700f
val goertzel = CwGoertzel()
goertzel.init(sampleRate, targetFreq)
// Generate 700 Hz tone
for (i in 0 until sampleRate.toInt()) {
goertzel.process((sin(2.0 * PI * targetFreq * i / sampleRate)).toFloat())
fun fft_magnitudeSpectrum_detectsTone() {
val fft = CwFFT(256)
val sampleRate = 8000f
val freq = 700f
val buffer = FloatArray(256) { (sin(2.0 * PI * freq * it / sampleRate)).toFloat() }
val mag = fft.magnitudeSpectrum(buffer)
// Peak should be at bin around 700 * 256 / 8000 ≈ 22.4
var maxBin = 0
var maxVal = 0f
for (i in mag.indices) {
if (mag[i] > maxVal) { maxVal = mag[i]; maxBin = i }
}
val power = goertzel.getPower()
assertTrue("Goertzel should detect present tone, got $power", power > 0.1f)
assertTrue("Peak bin $maxBin should be near 22", maxBin in 18..26)
assertTrue("Peak value $maxVal should be positive", maxVal > 0.01f)
}
@Test
fun goertzel_rejectsAbsentTone() {
val sampleRate = 4000f
val targetFreq = 700f
val goertzel = CwGoertzel()
goertzel.init(sampleRate, targetFreq)
// Generate 2000 Hz tone (no match)
for (i in 0 until sampleRate.toInt()) {
goertzel.process((sin(2.0 * PI * 2000f * i / sampleRate)).toFloat())
fun fft_magnitudeSpectrum_silence_isFlat() {
val fft = CwFFT(256)
val buffer = FloatArray(256) { 0f }
val mag = fft.magnitudeSpectrum(buffer)
for (v in mag) assertEquals("Silence spectrum should be 0, got $v", 0f, v, 1e-6f)
}
@Test
fun fft_rejectsWrongSize() {
assertThrows(IllegalArgumentException::class.java) { CwFFT(100) }
}
// --- Spectrogram ---
@Test
fun spectrogram_addSamples_updatesEnergy() {
val spec = CwSpectrogram(sampleRate = 8000)
val freq = 700f
// Feed multiple frames to stabilize energy normalization
for (i in 0..5) {
val buffer = FloatArray(256) { (sin(2.0 * PI * freq * it / 8000.0)).toFloat() }
spec.addSamples(buffer)
}
val power = goertzel.getPower()
assertTrue("Goertzel should reject absent tone, got $power", power < 0.1f)
val col = spec.getCurrentColumn()
val peakBin = spec.findPeakBin()
assertTrue("Peak bin $peakBin should be >= 0", peakBin >= 0)
}
@Test
fun goertzel_reset() {
val goertzel = CwGoertzel()
goertzel.init(4000f, 700f)
goertzel.process(1f)
goertzel.reset()
assertEquals(0f, goertzel.getPower(), 0.001f)
}
// --- Pitch detector ---
@Test
fun pitchDetector_findsCorrectFrequency() {
val sampleRate = 4000f
val detector = CwPitchDetector(sampleRate, 200f, 1200f, 10f)
val targetFreq = 700f
val buffer = FloatArray(sampleRate.toInt()) { (sin(2.0 * PI * targetFreq * it / sampleRate)).toFloat() }
val pitch = detector.findPitch(buffer)
assertNotNull("Pitch should be detected", pitch)
if (pitch != null) {
assertTrue("Detected pitch $pitch should be close to 700 Hz", pitch in 680f..720f)
fun spectrogram_findPeakBin_returnsValidBin() {
val spec = CwSpectrogram(sampleRate = 8000)
// Add multiple frames of 700 Hz tone
for (i in 0..5) {
val buffer = FloatArray(256) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() }
spec.addSamples(buffer)
}
val peakBin = spec.findPeakBin()
assertTrue("Peak bin should be >= 0, got $peakBin", peakBin >= 0)
}
@Test
fun pitchDetector_findsDifferentFrequency() {
val sampleRate = 4000f
val detector = CwPitchDetector(sampleRate, 200f, 1200f, 10f)
val targetFreq = 500f
val buffer = FloatArray(sampleRate.toInt()) { (sin(2.0 * PI * targetFreq * it / sampleRate)).toFloat() }
val pitch = detector.findPitch(buffer)
assertNotNull("Pitch should be detected", pitch)
if (pitch != null) {
assertTrue("Detected pitch $pitch should be close to 500 Hz", pitch in 480f..520f)
fun spectrogram_freqToBin_roundtrip() {
val spec = CwSpectrogram(sampleRate = 8000)
val freq = 700f
val bin = spec.freqToBin(freq)
val backFreq = spec.binToFreq(bin)
assertTrue("Freq $freq → bin $bin → freq $backFreq", backFreq > 600f && backFreq < 800f)
}
@Test
fun spectrogram_getBinEnergy_returnsCorrectLength() {
val spec = CwSpectrogram(sampleRate = 8000)
val energy = spec.getBinEnergy(0, 10)
assertEquals(10, energy.size)
}
@Test
fun spectrogram_reset() {
val spec = CwSpectrogram(sampleRate = 8000)
spec.addSamples(FloatArray(256) { 1f })
spec.reset()
assertEquals(-1, spec.findPeakBin())
}
// --- Bayesian decoder ---
@Test
fun bayesian_processTone_dit() {
val decoder = CwBayesianDecoder()
// At 20 WPM, dot = 60 ms
val result = decoder.processTone(60f)
assertEquals('0', result.symbol)
assertTrue("Dit probability should be positive", result.probability > 0.1f)
}
@Test
fun bayesian_processTone_dash() {
val decoder = CwBayesianDecoder()
// Dash = 3 * dot = 180 ms
val result = decoder.processTone(180f)
assertEquals('1', result.symbol)
assertTrue("Dash probability should be positive", result.probability > 0.1f)
}
@Test
fun bayesian_processTone_unknown_returnsNull() {
val decoder = CwBayesianDecoder()
// Very long tone — low probability for both dit and dash
val result = decoder.processTone(5000f)
assertNull(result.symbol)
}
@Test
fun bayesian_processGap_interChar_returnsChar() {
val decoder = CwBayesianDecoder()
decoder.processTone(60f) // dit
decoder.processTone(60f) // dit
decoder.processTone(60f) // dit
// 3 dots = "000" = 'S'
val char = decoder.processGap(180f) // 3 * dot = inter-char gap
assertEquals('S', char)
}
@Test
fun bayesian_processGap_wordGap_addsSpace() {
val decoder = CwBayesianDecoder()
decoder.processTone(60f) // dit = 'E'
decoder.processGap(180f) // inter-char gap
// Now word gap
val space = decoder.processGap(420f) // 7 * dot
assertEquals(' ', space)
}
@Test
fun bayesian_decodedText_accumulates() {
val decoder = CwBayesianDecoder()
decoder.processTone(60f) // dit = 'E'
decoder.processGap(180f) // inter-char
assertTrue(decoder.decodedText.isNotEmpty())
}
@Test
fun bayesian_reset() {
val decoder = CwBayesianDecoder()
decoder.processTone(60f)
decoder.reset()
assertEquals("", decoder.decodedText)
}
@Test
fun bayesian_getSpeed() {
val decoder = CwBayesianDecoder()
// Send 3 dits at 20 WPM (60 ms each)
decoder.processTone(60f)
decoder.processTone(60f)
decoder.processTone(60f)
val speed = decoder.getSpeed()
assertTrue("Speed should be ~20 WPM, got $speed", speed > 15f && speed < 30f)
}
// --- Channel tracker ---
@Test
fun channelTracker_initialState() {
val spec = CwSpectrogram(sampleRate = 8000)
val tracker = CwChannelTracker(spec)
val channels = tracker.update()
assertTrue("No channels should be active initially", channels.isEmpty())
}
@Test
fun channelTracker_detectsTone() {
val spec = CwSpectrogram(sampleRate = 8000)
// Feed a tone
for (i in 0..5) {
spec.addSamples(FloatArray(256) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() })
}
val tracker = CwChannelTracker(spec)
val channels = tracker.update()
assertTrue("Should detect at least 1 channel", channels.isNotEmpty())
}
@Test
fun pitchDetector_returnsNullForSilence() {
val detector = CwPitchDetector(4000f)
val buffer = FloatArray(4000) { 0f }
val pitch = detector.findPitch(buffer)
assertNull("Pitch should be null for silence", pitch)
fun channelTracker_bestChannel() {
val spec = CwSpectrogram(sampleRate = 8000)
for (i in 0..5) {
spec.addSamples(FloatArray(256) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() })
}
val tracker = CwChannelTracker(spec)
tracker.update()
val best = tracker.getBestChannel()
assertNotNull("Best channel should exist", best)
if (best != null) assertTrue(best.frequency in 600f..800f)
}
@Test
fun pitchDetector_emptyBuffer() {
val detector = CwPitchDetector(4000f)
assertNull(detector.findPitch(FloatArray(0)))
fun channelTracker_reset() {
val spec = CwSpectrogram(sampleRate = 8000)
for (i in 0..5) {
spec.addSamples(FloatArray(256) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() })
}
val tracker = CwChannelTracker(spec)
tracker.update()
tracker.reset()
assertNull(tracker.getBestChannel())
}
// --- Decoder state ---
// --- Full decoder ---
@Test
fun cwDecoder_initialState() {
fun decoder_initialState() {
val decoder = CwDecoder()
assertEquals("", decoder.decodedTextFlow.value)
assertEquals(0f, decoder.signalStrength.value, 0.001f)
assertNull(decoder.estimatedPitch.value)
assertNull(decoder.estimatedSpeed.value)
}
@Test
fun cwDecoder_defaultParameters() {
fun decoder_processSilence_doesNotCrash() {
val decoder = CwDecoder()
assertEquals(8000, decoder.sampleRate)
decoder.processBuffer(FloatArray(256) { 0f })
assertEquals("", decoder.decodedTextFlow.value)
}
@Test
fun cwDecoder_customParameters() {
val decoder = CwDecoder(sampleRate = 11025, cwToneFreq = 600f)
assertEquals(11025, decoder.sampleRate)
}
@Test
fun resetDecoder_clearsState() {
fun decoder_processNoise_doesNotCrash() {
val decoder = CwDecoder()
decoder.processBuffer(FloatArray(128) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() })
decoder.processBuffer(FloatArray(256) { (Math.random() * 2 - 1).toFloat() * 0.1f })
assertNotNull(decoder.decodedTextFlow.value)
}
@Test
fun decoder_processTone_doesNotCrash() {
val decoder = CwDecoder()
for (i in 0..20) {
decoder.processBuffer(FloatArray(256) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() })
}
assertNotNull(decoder.decodedTextFlow.value)
}
@Test
fun decoder_reset() {
val decoder = CwDecoder()
decoder.processBuffer(FloatArray(256) { 1f })
decoder.resetDecoder()
assertEquals("", decoder.decodedTextFlow.value)
assertEquals(0f, decoder.signalStrength.value, 0.001f)
}
@Test
fun processBuffer_silence_doesNotCrash() {
val decoder = CwDecoder()
val silence = FloatArray(1024) { 0f }
decoder.processBuffer(silence)
assertEquals("", decoder.decodedTextFlow.value)
}
@Test
fun processBuffer_noise_doesNotCrash() {
val decoder = CwDecoder()
val noise = FloatArray(1024) { (Math.random() * 2 - 1).toFloat() * 0.1f }
decoder.processBuffer(noise)
assertNotNull(decoder.decodedTextFlow.value)
}
@Test
fun processBuffer_withFixedPitch_doesNotCrash() {
fun decoder_withFixedPitch() {
val decoder = CwDecoder(sampleRate = 8000, cwToneFreq = 700f)
val buf = FloatArray(512) { (sin(2.0 * PI * 700.0 * it / 8000.0)).toFloat() }
decoder.processBuffer(buf)
assertNotNull(decoder.decodedTextFlow.value)
}
@Test
fun cwDecoder_withFixedPitchBypassesAutoDetect() {
val decoder = CwDecoder(sampleRate = 8000, cwToneFreq = 600f)
assertEquals(600f, decoder.estimatedPitch.value)
}
@Test
fun resetDecoder_afterFixedPitch() {
val decoder = CwDecoder(sampleRate = 8000, cwToneFreq = 700f)
decoder.resetDecoder()
assertEquals("", decoder.decodedTextFlow.value)
// Pitch should still be locked at 700
assertEquals(700f, decoder.estimatedPitch.value)
}
}