296 lines
10 KiB
Go
296 lines
10 KiB
Go
package audio
|
||
|
||
import (
|
||
"fmt"
|
||
"image"
|
||
"image/color"
|
||
"math"
|
||
"math/cmplx"
|
||
)
|
||
|
||
// Этот файл — единственное по-настоящему рискованное место во всём
|
||
// проекте, риск другого рода, чем обычно: код скомпилируется и
|
||
// отработает без паники, но КАЛИБРОВКА (тайминги строк Robot36 ниже) —
|
||
// по опубликованным спецификациям, без эталонной записи под рукой для
|
||
// сверки. Спектрограмма (см. spectrogram.go) не нуждается в точной
|
||
// калибровке — громче звук, ярче пиксель, и на любом входе получится
|
||
// осмысленная картинка. SSTV — противоположность: либо тайминги
|
||
// совпадают почти до сэмпла, либо на выходе шум/диагональные полосы
|
||
// вместо изображения. Если результат выглядит именно так — это первое
|
||
// место для подстройки, см. константы robot36*Ms ниже.
|
||
//
|
||
// Сюда же: авто-детект VIS-заголовка (11-битная последовательность
|
||
// тонов в начале передачи, определяющая режим SSTV) не реализован —
|
||
// это ещё один пласт сложности поверх итак рискованной части. Вместо
|
||
// этого просто ищем первый устойчивый синхроимпульс ~1200Гц и
|
||
// декодируем от него, предполагая заранее известный режим (Robot36).
|
||
|
||
const (
|
||
robot36SyncMs = 9.0
|
||
robot36PorchMs = 3.0
|
||
robot36YMs = 88.0
|
||
robot36SepMs = 4.5
|
||
robot36ChromaMs = 44.0
|
||
robot36LineMs = robot36SyncMs + robot36PorchMs + robot36YMs + robot36SepMs + robot36ChromaMs
|
||
|
||
robot36Width = 320
|
||
robot36Height = 240
|
||
)
|
||
|
||
// DecodeRobot36 декодирует Robot36 SSTV-сигнал в изображение 320x240.
|
||
func DecodeRobot36(wav *WAV) (*image.RGBA, error) {
|
||
if wav.SampleRate <= 0 {
|
||
return nil, fmt.Errorf("некорректная частота дискретизации: %d", wav.SampleRate)
|
||
}
|
||
|
||
padded := padToPow2(wav.Samples)
|
||
analytic := analyticSignal(padded)
|
||
freqRaw := instantaneousFrequency(analytic, wav.SampleRate)
|
||
freq := smooth(freqRaw[:len(wav.Samples)], 8) // обрезаем хвост от zero-padding обратно до исходной длины
|
||
|
||
samplesPerMs := float64(wav.SampleRate) / 1000
|
||
minSyncSamples := int(robot36SyncMs * 0.7 * samplesPerMs) // с запасом на неточность/шум — не требуем идеальных 9мс
|
||
|
||
startIdx := findSyncStart(freq, minSyncSamples)
|
||
if startIdx < 0 {
|
||
return nil, fmt.Errorf("не удалось найти начальный синхроимпульс SSTV (~1200Гц) в записи — возможно, это не Robot36, или запись не начинается прямо с сигнала")
|
||
}
|
||
|
||
img := image.NewRGBA(image.Rect(0, 0, robot36Width, robot36Height))
|
||
lineSamples := int(robot36LineMs * samplesPerMs)
|
||
|
||
pos := startIdx
|
||
var lastCr, lastCb []uint8
|
||
|
||
for line := 0; line < robot36Height; line++ {
|
||
if pos+lineSamples > len(freq) {
|
||
break // запись оборвалась раньше ожидаемого — отдаём то, что успели декодировать, не ошибку
|
||
}
|
||
|
||
yStart := pos + int((robot36SyncMs+robot36PorchMs)*samplesPerMs)
|
||
ySamples := decodeScan(freq, yStart, robot36YMs*samplesPerMs, robot36Width)
|
||
|
||
chromaStart := yStart + int(robot36YMs*samplesPerMs) + int(robot36SepMs*samplesPerMs)
|
||
chromaSamples := decodeScan(freq, chromaStart, robot36ChromaMs*samplesPerMs, robot36Width)
|
||
|
||
// Robot36 экономит полосу: Cr и Cb передаются по очереди через
|
||
// строку (2:1 вертикальный сабсэмплинг цветности), а не оба
|
||
// каждую строку — недостающий канал берём с предыдущей строки.
|
||
var cr, cb []uint8
|
||
if line%2 == 0 {
|
||
cr = chromaSamples
|
||
cb = lastCb
|
||
} else {
|
||
cb = chromaSamples
|
||
cr = lastCr
|
||
}
|
||
if cr == nil {
|
||
cr = chromaSamples
|
||
}
|
||
if cb == nil {
|
||
cb = chromaSamples
|
||
}
|
||
lastCr, lastCb = cr, cb
|
||
|
||
for x := 0; x < robot36Width; x++ {
|
||
r, g, b := ycrcbToRGB(ySamples[x], cr[x], cb[x])
|
||
img.Set(x, line, color.RGBA{R: r, G: g, B: b, A: 255})
|
||
}
|
||
|
||
pos += lineSamples
|
||
}
|
||
|
||
return img, nil
|
||
}
|
||
|
||
// decodeScan вытаскивает n значений яркости/цвета из отрезка частотного
|
||
// ряда длиной durationSamples, начиная с start — по одному значению на
|
||
// пиксель, усредняя частоту в соответствующем под-интервале (не берём
|
||
// один-единственный сэмпл — так меньше шума на выходе).
|
||
func decodeScan(freq []float64, start int, durationSamples float64, n int) []uint8 {
|
||
out := make([]uint8, n)
|
||
step := durationSamples / float64(n)
|
||
for i := 0; i < n; i++ {
|
||
s := start + int(float64(i)*step)
|
||
e := start + int(float64(i+1)*step)
|
||
if e <= s {
|
||
e = s + 1
|
||
}
|
||
if s < 0 {
|
||
s = 0
|
||
}
|
||
if e > len(freq) {
|
||
e = len(freq)
|
||
}
|
||
if s >= e {
|
||
continue
|
||
}
|
||
var sum float64
|
||
for j := s; j < e; j++ {
|
||
sum += freq[j]
|
||
}
|
||
out[i] = freqToValue(sum / float64(e-s))
|
||
}
|
||
return out
|
||
}
|
||
|
||
// freqToValue — стандартная для SSTV шкала: 1500Гц = 0, 2300Гц = 255,
|
||
// линейно между ними.
|
||
func freqToValue(f float64) uint8 {
|
||
v := (f - 1500) / (2300 - 1500) * 255
|
||
if v < 0 {
|
||
return 0
|
||
}
|
||
if v > 255 {
|
||
return 255
|
||
}
|
||
return uint8(v)
|
||
}
|
||
|
||
// ycrcbToRGB — стандартные коэффициенты ITU-R BT.601 (те же веса
|
||
// восприятия яркости, что и в stego.Grayscale).
|
||
func ycrcbToRGB(y, cr, cb uint8) (r, g, b uint8) {
|
||
yf := float64(y)
|
||
crf := float64(cr) - 128
|
||
cbf := float64(cb) - 128
|
||
return clamp8(yf + 1.402*crf), clamp8(yf - 0.344136*cbf - 0.714136*crf), clamp8(yf + 1.772*cbf)
|
||
}
|
||
|
||
func clamp8(v float64) uint8 {
|
||
if v < 0 {
|
||
return 0
|
||
}
|
||
if v > 255 {
|
||
return 255
|
||
}
|
||
return uint8(v)
|
||
}
|
||
|
||
// findSyncStart ищет первый устойчивый участок частоты в районе 1200Гц
|
||
// (±100Гц) длиной не меньше minSamples — начало SSTV-передачи. Без
|
||
// декодирования VIS-заголовка (см. doc файла) это единственная точка
|
||
// привязки: от нас найденного момента дальше строки идут строго по
|
||
// заранее известным Robot36-таймингам.
|
||
func findSyncStart(freq []float64, minSamples int) int {
|
||
run := 0
|
||
for i, f := range freq {
|
||
if f > 1100 && f < 1300 {
|
||
run++
|
||
if run >= minSamples {
|
||
return i - run + 1
|
||
}
|
||
} else {
|
||
run = 0
|
||
}
|
||
}
|
||
return -1
|
||
}
|
||
|
||
// --- Демодуляция: преобразование Гильберта → мгновенная фаза → её
|
||
// производная = мгновенная частота. Классическая техника для FM-сигналов
|
||
// (в отличие от спектрограммы, где годится оконное FFT — здесь окно
|
||
// физически не может быть узким и точным одновременно на масштабе
|
||
// одного пикселя SSTV, порядка десятков сэмплов).
|
||
|
||
// analyticSignal строит аналитический сигнал через преобразование
|
||
// Гильберта: FFT → обнулить отрицательные частоты и удвоить
|
||
// положительные (кроме DC/Найквиста) → обратное FFT. Стандартный рецепт,
|
||
// переиспользует уже написанный fft() из fft.go — в оба конца, через
|
||
// тождество IFFT(x) = conj(FFT(conj(x)))/N (см. ifft ниже), лишний код
|
||
// отдельной обратной реализации не пишем.
|
||
func analyticSignal(x []float64) []complex128 {
|
||
n := len(x)
|
||
buf := make([]complex128, n)
|
||
for i, v := range x {
|
||
buf[i] = complex(v, 0)
|
||
}
|
||
fft(buf)
|
||
|
||
for k := 1; k < n/2; k++ {
|
||
buf[k] *= 2
|
||
}
|
||
for k := n/2 + 1; k < n; k++ {
|
||
buf[k] = 0
|
||
}
|
||
|
||
ifft(buf)
|
||
return buf
|
||
}
|
||
|
||
func ifft(x []complex128) {
|
||
n := len(x)
|
||
for i := range x {
|
||
x[i] = cmplx.Conj(x[i])
|
||
}
|
||
fft(x)
|
||
for i := range x {
|
||
x[i] = cmplx.Conj(x[i]) / complex(float64(n), 0)
|
||
}
|
||
}
|
||
|
||
// instantaneousFrequency считает мгновенную частоту по аналитическому
|
||
// сигналу: разворачиваем фазу (убираем скачки на 2π при переходе через
|
||
// границу arctan), берём её производную, умножаем на sampleRate/(2π).
|
||
func instantaneousFrequency(analytic []complex128, sampleRate int) []float64 {
|
||
n := len(analytic)
|
||
if n == 0 {
|
||
return nil
|
||
}
|
||
|
||
unwrapped := make([]float64, n)
|
||
unwrapped[0] = cmplx.Phase(analytic[0])
|
||
prevPhase := unwrapped[0]
|
||
for i := 1; i < n; i++ {
|
||
phase := cmplx.Phase(analytic[i])
|
||
delta := phase - prevPhase
|
||
for delta > math.Pi {
|
||
delta -= 2 * math.Pi
|
||
}
|
||
for delta < -math.Pi {
|
||
delta += 2 * math.Pi
|
||
}
|
||
unwrapped[i] = unwrapped[i-1] + delta
|
||
prevPhase = phase
|
||
}
|
||
|
||
freq := make([]float64, n)
|
||
for i := 0; i < n-1; i++ {
|
||
freq[i] = (unwrapped[i+1] - unwrapped[i]) * float64(sampleRate) / (2 * math.Pi)
|
||
}
|
||
if n > 1 {
|
||
freq[n-1] = freq[n-2]
|
||
}
|
||
return freq
|
||
}
|
||
|
||
// smooth — простое скользящее среднее, сглаживает шум в сырой оценке
|
||
// мгновенной частоты перед тем, как переводить её в яркость пикселя.
|
||
func smooth(x []float64, window int) []float64 {
|
||
if window < 2 || len(x) == 0 {
|
||
return x
|
||
}
|
||
out := make([]float64, len(x))
|
||
var sum float64
|
||
for i := range x {
|
||
sum += x[i]
|
||
if i >= window {
|
||
sum -= x[i-window]
|
||
}
|
||
n := window
|
||
if i < window {
|
||
n = i + 1
|
||
}
|
||
out[i] = sum / float64(n)
|
||
}
|
||
return out
|
||
}
|
||
|
||
func padToPow2(x []float64) []float64 {
|
||
n := nextPowerOfTwo(len(x))
|
||
if n == len(x) {
|
||
return x
|
||
}
|
||
out := make([]float64, n)
|
||
copy(out, x)
|
||
return out
|
||
}
|