from __future__ import annotations import csv import hashlib import json import math import sqlite3 import struct import wave import zlib from pathlib import Path import matplotlib import numpy as np matplotlib.use("Agg") import matplotlib.pyplot as plt ROOT = Path(__file__).resolve().parent OUTPUTS = ROOT / "outputs" SAMPLE_RATE = 48_000 BPM = 96 BEAT_SECONDS = 60 / BPM SONG_BEATS = 64 SONG_SECONDS = SONG_BEATS * BEAT_SECONDS DATA_START_SECONDS = 2.5 BAUD = 1_000 SAMPLES_PER_BIT = SAMPLE_RATE // BAUD FREQUENCY_ZERO = 18_000 FREQUENCY_ONE = 20_000 CARRIER_AMPLITUDE = 0.018 REPETITION = 3 MAGIC = b"DBSONG1" HEADER_FORMAT = ">7sII32sI" PREAMBLE_BYTES = bytes([0x55] * 32) SYNC_BYTES = bytes.fromhex("d391c5a7") RNG_SEED = 20260720 def sha256_bytes(data: bytes) -> str: return hashlib.sha256(data).hexdigest() def bytes_to_bits(data: bytes) -> np.ndarray: return np.unpackbits(np.frombuffer(data, dtype=np.uint8)).astype(np.uint8) def bits_to_bytes(bits: np.ndarray) -> bytes: usable = bits[: len(bits) - (len(bits) % 8)] return np.packbits(usable).tobytes() def midi_frequency(note: int) -> float: return 440.0 * (2.0 ** ((note - 69) / 12.0)) def stereo_pan(signal: np.ndarray, pan: float) -> np.ndarray: angle = (pan + 1.0) * math.pi / 4.0 return np.column_stack((signal * math.cos(angle), signal * math.sin(angle))) def adsr(length: int, attack: float, decay: float, sustain: float, release: float) -> np.ndarray: if length <= 0: return np.zeros(0) attack_n = max(1, min(length, int(length * attack))) decay_n = max(1, min(length - attack_n, int(length * decay))) release_n = max(1, min(length - attack_n - decay_n, int(length * release))) sustain_n = max(0, length - attack_n - decay_n - release_n) return np.concatenate( [ np.linspace(0.0, 1.0, attack_n, endpoint=False), np.linspace(1.0, sustain, decay_n, endpoint=False), np.full(sustain_n, sustain), np.linspace(sustain, 0.0, release_n, endpoint=True), ] )[:length] def synth_note( note: int, seconds: float, amplitude: float, timbre: str = "pluck", phase: float = 0.0, ) -> np.ndarray: length = max(1, int(seconds * SAMPLE_RATE)) t = np.arange(length) / SAMPLE_RATE frequency = midi_frequency(note) if timbre == "bass": oscillator = ( np.sin(2 * np.pi * frequency * t + phase) + 0.22 * np.sin(2 * np.pi * frequency * 2 * t + phase) ) envelope = adsr(length, 0.02, 0.12, 0.62, 0.25) elif timbre == "pad": detune = frequency * 0.004 oscillator = ( np.sin(2 * np.pi * (frequency - detune) * t + phase) + np.sin(2 * np.pi * (frequency + detune) * t + phase) + 0.15 * np.sin(2 * np.pi * frequency * 2 * t) ) / 2.15 envelope = adsr(length, 0.18, 0.16, 0.72, 0.28) elif timbre == "lead": oscillator = ( np.sin(2 * np.pi * frequency * t + phase) + 0.28 * np.sin(2 * np.pi * frequency * 2 * t) + 0.08 * np.sin(2 * np.pi * frequency * 3 * t) ) / 1.36 envelope = adsr(length, 0.035, 0.12, 0.7, 0.24) else: oscillator = ( np.sin(2 * np.pi * frequency * t + phase) + 0.34 * np.sin(2 * np.pi * frequency * 2 * t) + 0.12 * np.sin(2 * np.pi * frequency * 3 * t) ) / 1.46 envelope = np.exp(-4.7 * t / max(seconds, 0.05)) envelope *= np.minimum(1.0, t / 0.008) return amplitude * oscillator * envelope def add_note( mix: np.ndarray, note: int, start_beat: float, beats: float, amplitude: float, pan: float, timbre: str, ) -> None: start = int(start_beat * BEAT_SECONDS * SAMPLE_RATE) signal = synth_note(note, beats * BEAT_SECONDS, amplitude, timbre) end = min(len(mix), start + len(signal)) if end > start: mix[start:end] += stereo_pan(signal[: end - start], pan) def add_kick(mix: np.ndarray, beat: float, amplitude: float = 0.34) -> None: length = int(0.34 * SAMPLE_RATE) t = np.arange(length) / SAMPLE_RATE phase = 2 * np.pi * (82 * t - 22 * t * t) signal = np.sin(phase) * np.exp(-13 * t) * amplitude start = int(beat * BEAT_SECONDS * SAMPLE_RATE) end = min(len(mix), start + length) mix[start:end] += stereo_pan(signal[: end - start], 0.0) def add_snare(mix: np.ndarray, beat: float, rng: np.random.Generator) -> None: length = int(0.22 * SAMPLE_RATE) t = np.arange(length) / SAMPLE_RATE noise = rng.normal(0, 1, length) noise = noise - np.convolve(noise, np.ones(15) / 15, mode="same") tone = np.sin(2 * np.pi * 180 * t) signal = (0.13 * noise + 0.08 * tone) * np.exp(-18 * t) start = int(beat * BEAT_SECONDS * SAMPLE_RATE) end = min(len(mix), start + length) mix[start:end] += stereo_pan(signal[: end - start], 0.08) def add_hat(mix: np.ndarray, beat: float, rng: np.random.Generator, open_hat: bool) -> None: seconds = 0.15 if open_hat else 0.055 length = int(seconds * SAMPLE_RATE) t = np.arange(length) / SAMPLE_RATE noise = rng.normal(0, 1, length) noise = noise - np.convolve(noise, np.ones(9) / 9, mode="same") signal = 0.035 * noise * np.exp((-28 if open_hat else -65) * t) start = int(beat * BEAT_SECONDS * SAMPLE_RATE) end = min(len(mix), start + length) mix[start:end] += stereo_pan(signal[: end - start], -0.18) def fft_lowpass(stereo: np.ndarray, cutoff_hz: float) -> np.ndarray: spectrum = np.fft.rfft(stereo, axis=0) frequencies = np.fft.rfftfreq(len(stereo), 1 / SAMPLE_RATE) spectrum[frequencies > cutoff_hz] = 0 return np.fft.irfft(spectrum, n=len(stereo), axis=0) def compose_song() -> np.ndarray: rng = np.random.default_rng(RNG_SEED) length = int(SONG_SECONDS * SAMPLE_RATE) mix = np.zeros((length, 2), dtype=np.float64) progression = [ ([60, 64, 67, 71], 48), ([57, 60, 64, 67], 45), ([53, 57, 60, 64], 41), ([55, 59, 62, 65], 43), ] melody = [ (0.0, 76, 1.0), (1.0, 79, 0.5), (1.5, 81, 0.5), (2.0, 79, 1.0), (3.0, 76, 1.0), (4.0, 72, 1.0), (5.0, 76, 0.5), (5.5, 79, 0.5), (6.0, 76, 1.0), (7.0, 72, 1.0), (8.0, 69, 1.0), (9.0, 72, 0.5), (9.5, 76, 0.5), (10.0, 74, 1.0), (11.0, 72, 1.0), (12.0, 71, 0.5), (12.5, 74, 0.5), (13.0, 79, 1.0), (14.0, 77, 0.5), (14.5, 74, 0.5), (15.0, 71, 1.0), ] for cycle in range(4): cycle_start = cycle * 16 for bar, (chord, bass_note) in enumerate(progression): bar_start = cycle_start + bar * 4 for note_index, note in enumerate(chord): add_note(mix, note, bar_start, 4.0, 0.055, -0.28 + note_index * 0.18, "pad") for beat_offset in range(4): add_note(mix, bass_note, bar_start + beat_offset, 0.85, 0.12, -0.12, "bass") arpeggio = chord + chord[1:3] for step in range(8): note = arpeggio[step % len(arpeggio)] + 12 add_note(mix, note, bar_start + step * 0.5, 0.42, 0.07, 0.25, "pluck") if cycle > 0: transpose = 12 if cycle == 3 else 0 variation = 2 if cycle == 2 else 0 for start, note, beats in melody: add_note( mix, note + transpose + (variation if int(start * 2) % 7 == 0 else 0), cycle_start + start, beats * 0.88, 0.09 if cycle < 3 else 0.075, 0.15, "lead", ) for beat in np.arange(0, SONG_BEATS, 0.5): add_hat(mix, float(beat), rng, open_hat=int(beat * 2) % 8 == 7) for beat in range(SONG_BEATS): if beat % 4 in (0, 2): add_kick(mix, beat, 0.30 if beat < 16 else 0.36) if beat % 4 in (1, 3): add_snare(mix, beat, rng) delay_samples = int(0.31 * SAMPLE_RATE) delayed = np.zeros_like(mix) delayed[delay_samples:] = mix[:-delay_samples] * 0.16 delayed[:, [0, 1]] = delayed[:, [1, 0]] mix += delayed mix = fft_lowpass(mix, 11_500) peak = float(np.max(np.abs(mix))) mix *= 0.78 / peak fade = int(0.8 * SAMPLE_RATE) mix[:fade] *= np.linspace(0, 1, fade)[:, None] mix[-fade:] *= np.linspace(1, 0, fade)[:, None] return mix def create_database(path: Path) -> None: if path.exists(): path.unlink() connection = sqlite3.connect(path) connection.execute("PRAGMA page_size=512") connection.execute("PRAGMA journal_mode=DELETE") connection.execute("VACUUM") connection.executescript( """ CREATE TABLE track ( id INTEGER PRIMARY KEY, title TEXT NOT NULL, bpm INTEGER NOT NULL, key TEXT NOT NULL, duration_seconds REAL NOT NULL, data_carrier TEXT NOT NULL ); CREATE TABLE arrangement ( section TEXT PRIMARY KEY, start_beat INTEGER NOT NULL, chord_progression TEXT NOT NULL, role TEXT NOT NULL ); CREATE TABLE liner_note ( id INTEGER PRIMARY KEY, note TEXT NOT NULL ); """ ) connection.execute( "INSERT INTO track VALUES (?, ?, ?, ?, ?, ?)", (1, "Side Channel", BPM, "C major / A minor", SONG_SECONDS, "18 kHz / 20 kHz stereo-side BFSK"), ) arrangement = [ ("intro", 0, "Cmaj7 Am7 Fmaj7 G7", "establish the progression"), ("theme", 16, "Cmaj7 Am7 Fmaj7 G7", "introduce the lead melody"), ("variation", 32, "Cmaj7 Am7 Fmaj7 G7", "alter selected melody tones"), ("finale", 48, "Cmaj7 Am7 Fmaj7 G7", "lift the melody by one octave"), ] connection.executemany("INSERT INTO arrangement VALUES (?, ?, ?, ?)", arrangement) notes = [ "The audible arrangement is original and synthesized by the experiment.", "The database is compressed, repeated three times per bit, and carried in stereo difference.", "Converting the track to mono is expected to remove the payload.", "If this table is readable, the song has successfully returned its own paperwork.", ] connection.executemany("INSERT INTO liner_note(note) VALUES (?)", [(note,) for note in notes]) connection.commit() connection.execute("VACUUM") connection.close() def make_packet(database_bytes: bytes) -> tuple[bytes, bytes]: compressed = zlib.compress(database_bytes, level=9) checksum = zlib.crc32(compressed) & 0xFFFFFFFF header = struct.pack( HEADER_FORMAT, MAGIC, len(database_bytes), len(compressed), hashlib.sha256(database_bytes).digest(), checksum, ) return header + compressed, compressed def encode_frame(packet: bytes) -> np.ndarray: prefix = bytes_to_bits(PREAMBLE_BYTES + SYNC_BYTES) body = np.repeat(bytes_to_bits(packet), REPETITION) return np.concatenate((prefix, body)) def modulate_bits(bits: np.ndarray) -> np.ndarray: frequencies = np.where(bits == 0, FREQUENCY_ZERO, FREQUENCY_ONE) sample_frequencies = np.repeat(frequencies, SAMPLES_PER_BIT) phase = 2 * np.pi * np.cumsum(sample_frequencies) / SAMPLE_RATE return CARRIER_AMPLITUDE * np.sin(phase) def embed_carrier(music: np.ndarray, bits: np.ndarray) -> tuple[np.ndarray, np.ndarray]: carrier = modulate_bits(bits) start = int(DATA_START_SECONDS * SAMPLE_RATE) if start + len(carrier) > len(music): raise ValueError("The encoded database does not fit inside the song.") result = music.copy() result[start : start + len(carrier), 0] += carrier result[start : start + len(carrier), 1] -= carrier return result, carrier def float_to_pcm16(stereo: np.ndarray) -> np.ndarray: clipped = np.clip(stereo, -0.999, 0.999) return np.round(clipped * 32767).astype(" None: pcm = float_to_pcm16(stereo) with wave.open(str(path), "wb") as output: output.setnchannels(2) output.setsampwidth(2) output.setframerate(SAMPLE_RATE) output.writeframes(pcm.tobytes()) def read_wav(path: Path) -> np.ndarray: with wave.open(str(path), "rb") as input_file: if input_file.getnchannels() != 2 or input_file.getsampwidth() != 2: raise ValueError("Expected 16-bit stereo PCM.") frames = input_file.readframes(input_file.getnframes()) return np.frombuffer(frames, dtype=" np.ndarray: start = int(DATA_START_SECONDS * SAMPLE_RATE) required = bit_count * SAMPLES_PER_BIT side = (stereo[:, 0] - stereo[:, 1]) / 2 segment = side[start : start + required] if len(segment) != required: raise ValueError("Audio ended before the encoded frame.") windows = segment.reshape(bit_count, SAMPLES_PER_BIT) t = np.arange(SAMPLES_PER_BIT) / SAMPLE_RATE def energy(frequency: float) -> np.ndarray: sine = np.sin(2 * np.pi * frequency * t) cosine = np.cos(2 * np.pi * frequency * t) return (windows @ sine) ** 2 + (windows @ cosine) ** 2 return (energy(FREQUENCY_ONE) > energy(FREQUENCY_ZERO)).astype(np.uint8) def decode_frame(stereo: np.ndarray, frame_bit_count: int, output_path: Path | None = None) -> dict: decoded = demodulate_bits(stereo, frame_bit_count) prefix_bits = len(PREAMBLE_BYTES + SYNC_BYTES) * 8 expected_prefix = bytes_to_bits(PREAMBLE_BYTES + SYNC_BYTES) prefix_errors = int(np.count_nonzero(decoded[:prefix_bits] != expected_prefix)) repeated = decoded[prefix_bits:] usable = repeated[: len(repeated) - (len(repeated) % REPETITION)] groups = usable.reshape(-1, REPETITION) majority = (np.sum(groups, axis=1) >= (REPETITION // 2 + 1)).astype(np.uint8) packet_bytes = bits_to_bytes(majority) header_size = struct.calcsize(HEADER_FORMAT) result = { "prefix_bit_errors": prefix_errors, "decoded": False, "sha256_match": False, "integrity_check": None, } if len(packet_bytes) < header_size: result["error"] = "header truncated" return result try: magic, original_length, compressed_length, expected_hash, expected_crc = struct.unpack( HEADER_FORMAT, packet_bytes[:header_size] ) if magic != MAGIC: result["error"] = "magic mismatch" return result compressed = packet_bytes[header_size : header_size + compressed_length] if len(compressed) != compressed_length: result["error"] = "compressed payload truncated" return result if (zlib.crc32(compressed) & 0xFFFFFFFF) != expected_crc: result["error"] = "compressed payload CRC mismatch" return result database_bytes = zlib.decompress(compressed) if len(database_bytes) != original_length: result["error"] = "database length mismatch" return result digest_match = hashlib.sha256(database_bytes).digest() == expected_hash result.update( { "decoded": True, "database_bytes": len(database_bytes), "sha256": sha256_bytes(database_bytes), "sha256_match": digest_match, } ) if output_path is not None and digest_match: output_path.write_bytes(database_bytes) connection = sqlite3.connect(output_path) integrity = connection.execute("PRAGMA integrity_check").fetchone()[0] recovered_title = connection.execute("SELECT title FROM track WHERE id = 1").fetchone()[0] liner_note_count = connection.execute("SELECT COUNT(*) FROM liner_note").fetchone()[0] connection.close() result.update( { "integrity_check": integrity, "recovered_title": recovered_title, "liner_note_count": liner_note_count, } ) return result except (struct.error, zlib.error, sqlite3.DatabaseError) as error: result["error"] = f"{type(error).__name__}: {error}" return result def add_noise(stereo: np.ndarray, snr_db: float, rng: np.random.Generator) -> np.ndarray: rms = float(np.sqrt(np.mean(stereo**2))) noise_rms = rms / (10 ** (snr_db / 20)) return np.clip(stereo + rng.normal(0, noise_rms, stereo.shape), -0.999, 0.999) def quantize(stereo: np.ndarray, bits: int) -> np.ndarray: levels = (2 ** (bits - 1)) - 1 return np.round(np.clip(stereo, -1, 1) * levels) / levels def add_echo(stereo: np.ndarray, delay_ms: float, gain: float) -> np.ndarray: delay = int(delay_ms * SAMPLE_RATE / 1000) output = stereo.copy() output[delay:] += stereo[:-delay] * gain peak = np.max(np.abs(output)) if peak >= 1: output *= 0.98 / peak return output def resample_roundtrip(stereo: np.ndarray, intermediate_rate: int) -> np.ndarray: old_positions = np.arange(len(stereo)) down_length = int(round(len(stereo) * intermediate_rate / SAMPLE_RATE)) down_positions = np.linspace(0, len(stereo) - 1, down_length) down = np.column_stack( [np.interp(down_positions, old_positions, stereo[:, channel]) for channel in range(2)] ) up_positions = np.linspace(0, down_length - 1, len(stereo)) return np.column_stack( [np.interp(up_positions, np.arange(down_length), down[:, channel]) for channel in range(2)] ) def make_spectrum_figure(stereo: np.ndarray, output_path: Path) -> None: start = int(DATA_START_SECONDS * SAMPLE_RATE) length = int(8 * SAMPLE_RATE) excerpt = stereo[start : start + length] mid = (excerpt[:, 0] + excerpt[:, 1]) / 2 side = (excerpt[:, 0] - excerpt[:, 1]) / 2 window = np.hanning(len(excerpt)) frequencies = np.fft.rfftfreq(len(excerpt), 1 / SAMPLE_RATE) mid_spectrum = 20 * np.log10(np.maximum(np.abs(np.fft.rfft(mid * window)), 1e-9)) side_spectrum = 20 * np.log10(np.maximum(np.abs(np.fft.rfft(side * window)), 1e-9)) mid_spectrum -= np.max(mid_spectrum) side_spectrum -= np.max(side_spectrum) figure, axes = plt.subplots(2, 1, figsize=(16, 9), constrained_layout=True) preview = stereo[: int(8 * SAMPLE_RATE) : 120] preview_time = np.arange(len(preview)) * 120 / SAMPLE_RATE axes[0].plot(preview_time, preview[:, 0], linewidth=0.7, label="left") axes[0].plot(preview_time, preview[:, 1], linewidth=0.7, alpha=0.72, label="right") axes[0].set(title="First eight seconds of the finished stereo song", xlabel="seconds", ylabel="amplitude") axes[0].legend() axes[0].grid(alpha=0.2) axes[1].plot(frequencies / 1000, mid_spectrum, label="stereo mid: audible music") axes[1].plot(frequencies / 1000, side_spectrum, label="stereo side: music + database carrier") axes[1].axvline(FREQUENCY_ZERO / 1000, color="#d62728", linestyle="--", linewidth=1) axes[1].axvline(FREQUENCY_ONE / 1000, color="#9467bd", linestyle="--", linewidth=1) axes[1].set( title="The database appears as energy around 18 kHz and 20 kHz in the stereo side channel", xlabel="frequency (kHz)", ylabel="relative magnitude (dB)", xlim=(0, 24), ylim=(-100, 5), ) axes[1].legend() axes[1].grid(alpha=0.2) figure.savefig(output_path, dpi=120) plt.close(figure) def run() -> None: OUTPUTS.mkdir(parents=True, exist_ok=True) source_database = OUTPUTS / "source.sqlite" recovered_database = OUTPUTS / "recovered.sqlite" song_path = OUTPUTS / "side-channel.wav" create_database(source_database) database_bytes = source_database.read_bytes() packet, compressed = make_packet(database_bytes) frame_bits = encode_frame(packet) music = compose_song() song, carrier = embed_carrier(music, frame_bits) write_wav(song_path, song) encoded_song = read_wav(song_path) clean_result = decode_frame(encoded_song, len(frame_bits), recovered_database) if not ( clean_result.get("sha256_match") and clean_result.get("integrity_check") == "ok" and recovered_database.read_bytes() == database_bytes ): raise RuntimeError(f"Clean recovery failed: {clean_result}") rng = np.random.default_rng(RNG_SEED + 1) cases = { "clean": encoded_song, "noise_50_db": add_noise(encoded_song, 50, rng), "noise_40_db": add_noise(encoded_song, 40, rng), "noise_30_db": add_noise(encoded_song, 30, rng), "noise_25_db": add_noise(encoded_song, 25, rng), "quantized_12_bit": quantize(encoded_song, 12), "quantized_8_bit": quantize(encoded_song, 8), "clipped_0_35": np.clip(encoded_song, -0.35, 0.35), "echo_12_ms_0_35": add_echo(encoded_song, 12, 0.35), "resample_44_1khz": resample_roundtrip(encoded_song, 44_100), "lowpass_15khz": fft_lowpass(encoded_song, 15_000), "mono_collapse": np.repeat(np.mean(encoded_song, axis=1)[:, None], 2, axis=1), } benchmark = [] for name, audio in cases.items(): decoded_bits = demodulate_bits(audio, len(frame_bits)) bit_errors = int(np.count_nonzero(decoded_bits != frame_bits)) result = decode_frame(audio, len(frame_bits)) benchmark.append( { "case": name, "bit_errors": bit_errors, "bit_error_rate": bit_errors / len(frame_bits), "database_decoded": bool(result.get("decoded")), "sha256_match": bool(result.get("sha256_match")), "failure": result.get("error", ""), } ) with (OUTPUTS / "benchmark.csv").open("w", newline="", encoding="utf-8") as output: writer = csv.DictWriter(output, fieldnames=benchmark[0].keys()) writer.writeheader() writer.writerows(benchmark) music_rms = float(np.sqrt(np.mean(music**2))) carrier_rms = float(np.sqrt(np.mean(carrier**2))) manifest = { "format": "database-song-experiment/v1", "title": "Side Channel", "sample_rate_hz": SAMPLE_RATE, "channels": 2, "duration_seconds": SONG_SECONDS, "tempo_bpm": BPM, "database_bytes": len(database_bytes), "database_sha256": sha256_bytes(database_bytes), "compressed_bytes": len(compressed), "packet_bytes": len(packet), "frame_bits": len(frame_bits), "fec": f"repetition-{REPETITION} majority vote", "carrier": { "modulation": "binary frequency-shift keying", "channel": "(left - right) / 2", "zero_hz": FREQUENCY_ZERO, "one_hz": FREQUENCY_ONE, "baud": BAUD, "amplitude_peak": CARRIER_AMPLITUDE, "start_seconds": DATA_START_SECONDS, "duration_seconds": len(carrier) / SAMPLE_RATE, "rms_db_below_music": 20 * math.log10(carrier_rms / music_rms), }, "audio": { "peak": float(np.max(np.abs(song))), "rms": float(np.sqrt(np.mean(song**2))), "wav_bytes": song_path.stat().st_size, }, "clean_recovery": clean_result, "benchmark_cases": len(benchmark), "benchmark_exact_recoveries": sum(row["sha256_match"] for row in benchmark), "limitations": [ "The decoder uses a known carrier start time and symbol rate.", "The benchmark applies deterministic digital transformations, not a speaker-to-microphone loopback.", "Mono conversion and low-pass filtering below the carrier intentionally destroy the payload.", "Listenability is a human judgment; the manifest records signal levels rather than claiming a perceptual score.", ], } (OUTPUTS / "manifest.json").write_text(json.dumps(manifest, indent=2), encoding="utf-8") make_spectrum_figure(encoded_song, OUTPUTS / "waveform-spectrum.png") print(json.dumps(manifest, indent=2)) if __name__ == "__main__": run()