""" Vertical Tracking Angle (VTA) & Tracing Distortion Analyzer ============================================================ Software implementation of US Patent 4,359,768 (CBS Inc., 1982) "Vertical Tracking Angle Meter" — Abbagnaro & Gust Version 2.0 — correct sign handling, dual-channel spectrum, calibrated V400 Signal chain (Fig. 4 of patent): INPUT (RAW velocity pickup output — NO RIAA equalization applied) │ ├─── SIGNAL PATH A ──────────────────────────────────────────────────── │ [44] 2kHz HPF + 400Hz NOTCH → isolates 4kHz carrier + sidebands │ [46] Axis Crossing Detector → 4kHz carrier → 4kHz square wave │ [48] FM Discriminator → demodulates 4kHz → 400Hz dev signal │ [50] LPF → isolates 400Hz component of FM dev │ [52] Phase Compensator → corrects HPF residual phase shift │ → 400Hz DEVIATION SIGNAL (amplitude ∝ FM deviation F) │ └─── SIGNAL PATH B ──────────────────────────────────────────────────── [54] LPF → removes 4kHz, passes 400Hz only [55] Switch (0° or 90°) → selects tracking or tracing mode [57] 90° Phase Shift (opt.) → inserted for tracing distortion mode [56] Axis Crossing Detector → 400Hz sine → 400Hz square wave → 400Hz REFERENCE SQUARE WAVE CHOPPER [58]: deviation_signal × reference_square_wave → C₀ waveform (full-wave rectified) when signals in phase → TRACKING → C₉₀ waveform (zero average) when in quadrature → TRACING [60] LPF + OFFSET → DC proportional to VTA (offset sets zero at θR=16.0°) [62] METER → reads VTA directly in degrees VTA Formula (patent Section 3 / Col. 3-4): θP = arctan[ tan(θR) ± (9.1×10⁻³ · R · F · cos(φ)) / V₄₀₀ ] θP = vertical playback angle (degrees) θR = vertical recording angle (16.0° for CBS STR-112 Group 2B) R = groove radius (metres) F = peak 400Hz FM deviation of 4kHz tone (Hz) φ = phase lead of E2 w.r.t. E1 (degrees) [0° = in-phase] V₄₀₀= peak velocity of recorded 400Hz tone (m/s) − = left channel, + = right channel CBS STR-112 Test Record — Group 2B (vertical modulation, 400+4000 Hz): Band B1: +6 dB, Band B2: +9 dB, Band B3: +12 dB Vertical recording angle θR = 16.0° Usage: python vta_analyzer.py recording.wav [options] python vta_analyzer.py --demo # synthetic signal, no WAV needed Requirements: pip install numpy scipy matplotlib pip install soundfile # for WAV input (optional) """ import sys import argparse import warnings import numpy as np import matplotlib.pyplot as plt import matplotlib.gridspec as gridspec from scipy import signal as sig # ── Try soundfile, fall back gracefully ────────────────────────────────────── try: import soundfile as sf _HAS_SF = True except ImportError: _HAS_SF = False # ═════════════════════════════════════════════════════════════════════════════ # ███████╗███████╗████████╗████████╗██╗███╗ ██╗ ██████╗ ███████╗ # ██╔════╝██╔════╝╚══██╔══╝╚══██╔══╝██║████╗ ██║██╔════╝ ██╔════╝ # ███████╗█████╗ ██║ ██║ ██║██╔██╗ ██║██║ ███╗███████╗ # ╚════██║██╔══╝ ██║ ██║ ██║██║╚██╗██║██║ ██║╚════██║ # ███████║███████╗ ██║ ██║ ██║██║ ╚████║╚██████╔╝███████║ # ╚══════╝╚══════╝ ╚═╝ ╚═╝ ╚═╝╚═╝ ╚═══╝ ╚═════╝ ╚══════╝ # # ┌─────────────────────────────────────────────────────────────────────┐ # │ USER SETTINGS — EDIT HERE BEFORE RUNNING │ # │ Press F5 in Spyder to run │ # └─────────────────────────────────────────────────────────────────────┘ # ═════════════════════════════════════════════════════════════════════════════ # Set DEMO_MODE = True to run without a WAV file (generates a synthetic signal) # Set DEMO_MODE = False to analyse your own recording DEMO_MODE = True # ── Your WAV file (only used when DEMO_MODE = False) ───────────────────────── # Windows: r"C:\Users\YourName\recordings\track2.wav" # Mac/Linux: "/home/yourname/recordings/track2.wav" AUDIO_FILE = r"C:\path\to\your\recording.wav" # ── Channel polarity ───────────────────────────────────────────────────────── # The patent (Col 7) notes that L/R channel wiring affects the sign convention. # If your VTA results look wrong (too low, near 0° or negative) swap this. # How to check: run the script, if L channel reads ~6° and R reads ~24° # (or vice versa) set SWAP_CHANNELS = True to correct it. # This is equivalent to the front-panel channel switch on the hardware meter. SWAP_CHANNELS = False # ── Which stereo channel(s) to analyse ─────────────────────────────────────── # 'L' = left channel only # 'R' = right channel only # 'both' = analyse both and average (recommended — matches patent Col. 7) CHANNEL = 'both' # ── Tone frequencies of the track you recorded ────────────────────────────── # CBS STR-112 standard bands: F_MOD=400.0, F_CARRIER=4000.0 # CBS STR-112 200Hz bands: F_MOD=200.0, F_CARRIER=4000.0 # Other test records: set both to match your record F_MOD = 400.0 # Low modulating / reference frequency (Hz) F_CARRIER = 4000.0 # High carrier frequency (Hz) — typically 10x F_MOD # ── Groove radius in METRES ─────────────────────────────────────────────────── # Measure from the centre spindle hole to where the stylus sits on the track. # Use a ruler on the record while it is on the turntable. # Outer edge of record ~145 mm → 0.145 # Typical Band B1 → 0.100 # Mid record → 0.080 (default) # Inner bands ~60 mm → 0.060 GROOVE_RADIUS = 0.080 # ── Peak recorded velocity of the 400Hz tone (m/s) ─────────────────────────── # This is V400 in the patent formula. It must reflect what the cartridge # actually outputs, NOT the displacement stated on the record label. # # How it was derived for CBS STR-112 Band B2 (+6dB vertical): # 1. Label states 4kHz tone at -18dB re 11.2µm displacement # → 4kHz peak displacement = 1.41µm # → 4kHz peak VELOCITY = 2π × 4000 × 1.41e-6 = 3.544 cm/s # 2. Measured amplitude ratio from WAV: 400Hz is 9.31× louder than 4kHz # (cartridge outputs voltage proportional to velocity, and at 4kHz # produces 10× more voltage per µm than at 400Hz — so the WAV ratio # correctly reflects the true velocity ratio seen by the discriminator) # 3. V400 = 9.31 × 3.544 cm/s = 32.98 cm/s = 0.3298 m/s # # Validated: script gives ~24° VTA, matching cartridge spec of 23° and # independent DIN record measurement of 22-23°. ✓ # # ───────────────────────────────────────────────────────────────────── # QUICK REFERENCE — set PEAK_VEL to match your record and band: # # CBS STR-112 400Hz + 4kHz (flat recording — NO RIAA) # Vertical Band +6dB → PEAK_VEL = 0.3299 ← VALIDATED ✓ (matches 23° cart) # Vertical Band +9dB → PEAK_VEL = 0.4660 # Vertical Band +12dB → PEAK_VEL = 0.6583 # # CBS STR-112 200Hz + 4kHz (flat recording — NO RIAA) # Same displacement levels as 400Hz bands but V = 2pi*200*x (half velocity) # Band +6dB → PEAK_VEL = 0.02808 # Band +9dB → PEAK_VEL = 0.03967 # Band +12dB → PEAK_VEL = 0.05603 # Also set: F_MOD=200.0, F_CARRIER=4000.0 # # ANALOG MAGIC 60Hz + 7kHz (RIAA recorded — special case) # RIAA recording boosts 60Hz groove velocity ~+16dB vs flat recording. # V_low cannot be calculated from label data alone. # CALIBRATION PROCEDURE (use your known 23deg cartridge): # 1. Record the Analog Magic track as WAV, NO RIAA playback # 2. Set F_MOD=60, F_CARRIER=7000 # 3. Start with PEAK_VEL = 1.0 (arbitrary), note VTA reading X # 4. Correct value: PEAK_VEL = 1.0 * tan(X deg) / tan(23 deg) # (iterate once more to converge) # Also set: F_MOD=60.0, F_CARRIER=7000.0 # IMPORTANT: play Analog Magic WITHOUT RIAA — the cartridge output # already has the RIAA shape baked in from the groove, the script # processes it correctly as long as PEAK_VEL is calibrated to match. # ───────────────────────────────────────────────────────────────────── PEAK_VEL = 0.3299 # m/s — CHANGE THIS to match your record/band (see above) # ── Known vertical recording angle of your test record (degrees) ───────────── RECORDING_ANGLE = 16.0 # ── Plots to show ───────────────────────────────────────────────────────────── SHOW_SPECTRUM = True # FFT of input — verify tones are present and levels OK SHOW_SIGNAL_CHAIN = True # Each stage of the patent signal chain (diagnostic) SHOW_METER = True # Meter face + tracing distortion harmonics # ── Demo mode only: VTA error to inject (degrees) ──────────────────────────── DEMO_VTA_ERROR = 5.0 # ═════════════════════════════════════════════════════════════════════════════ # END OF USER SETTINGS — do not edit below this line # ═════════════════════════════════════════════════════════════════════════════ # ═════════════════════════════════════════════════════════════════════════════ # RECORD / INSTRUMENT CONSTANTS (adjust to match your setup) # ═════════════════════════════════════════════════════════════════════════════ THETA_R_DEG = 16.0 # Verified recording angle of this test record (degrees) # F_LOW and F_HIGH are set in the USER SETTINGS block at the top of this file # as F_MOD and F_CARRIER — do not add them here GROOVE_RADIUS_M = 0.080 # Groove radius, metres (~80 mm = mid-side typical) # Peak velocity is set in USER SETTINGS above as PEAK_VEL — do not redefine here # Calibrated value: 0.3298 m/s for CBS STR-112 Band B2 +6dB vertical # Calibrate from 4kHz-only reference band if available # ═════════════════════════════════════════════════════════════════════════════ # UTILITY: FILTERS # ═════════════════════════════════════════════════════════════════════════════ def _sos_bp(fs, fc, bw, order=4): lo = max(fc - bw/2, 1.0) hi = min(fc + bw/2, fs/2 - 1) return sig.butter(order, [lo, hi], btype='bandpass', fs=fs, output='sos') def _sos_hp(fs, fc, order=4): return sig.butter(order, min(fc, fs/2-1), btype='high', fs=fs, output='sos') def _sos_lp(fs, fc, order=4): return sig.butter(order, min(fc, fs/2-1), btype='low', fs=fs, output='sos') def _sos_notch(fs, fc, Q=30): b, a = sig.iirnotch(fc, Q, fs) return sig.tf2sos(b, a) def apply(sos, x): return sig.sosfiltfilt(sos, x) # ═════════════════════════════════════════════════════════════════════════════ # SIGNAL PATH A — FM discriminator chain # ═════════════════════════════════════════════════════════════════════════════ def path_a(audio: np.ndarray, fs: float, f_low=F_MOD, f_high=F_CARRIER) -> np.ndarray: """ Implements blocks [44]→[46]→[48]→[50]→[52] of Fig. 4. [44] 2kHz HPF + 400Hz notch → carrier band only [46] Axis crossing detector → 4kHz square wave (implicit in Hilbert FM) [48] FM discriminator → instantaneous frequency deviation [50] LPF → isolate 400Hz deviation component [52] Phase compensator → correct HPF residual phase (auto-calibrated) Returns: 400Hz deviation signal (amplitude ∝ F, the peak FM deviation Hz) """ # [44] 2kHz HPF + 400Hz notch (let sidebands around 4kHz pass) hpf = apply(_sos_hp(fs, 2000.0, order=4), audio) notch = apply(_sos_notch(fs, f_low, Q=20), hpf) # [46]+[48] Zero-crossing FM discriminator via instantaneous frequency # The axis-crossing detector converts the carrier to a square wave; # this is equivalent to taking the sign() then differentiating — in # DSP the Hilbert-transform approach gives the same instantaneous freq. analytic = sig.hilbert(notch) inst_ph = np.unwrap(np.angle(analytic)) inst_freq = np.concatenate([[0], np.diff(inst_ph)]) * fs / (2*np.pi) deviation = inst_freq - f_high # FM deviation in Hz (baseband) # [50] LPF — isolate f_low component of FM deviation dev_filt = apply(_sos_lp(fs, f_low * 2.5, order=4), deviation) # Note: block [52] phase compensator (hardware: RC network to correct HPF # residual phase at f_high) is NOT needed in DSP — the Hilbert-transform # FM discriminator has no residual phase error at the modulating frequency. # Applying software phase rotation here would destroy the phase information # that the chopper needs to separate tracking from tracing components. return dev_filt # units: Hz of FM deviation def _phase_compensate(dev: np.ndarray, fs: float, f_low: float) -> np.ndarray: """ Block [52]: correct residual phase introduced by the HPF section of [44]. Implemented as a fractional-sample circular shift via FFT phase rotation. """ t = np.arange(len(dev)) / fs ref = np.sin(2*np.pi*f_low*t) # find phase of f_low component cc = np.dot(dev, ref) / len(dev) ss = np.dot(dev, np.cos(2*np.pi*f_low*t)) / len(dev) phi_rad = np.arctan2(cc, ss) # phase error # rotate in frequency domain D = np.fft.rfft(dev) freqs = np.fft.rfftfreq(len(dev), 1/fs) # apply correction only near f_low; elsewhere leave unchanged # (hardware phase compensator is wideband — approximate that here) correction = np.exp(-1j * phi_rad * np.exp(-((freqs - f_low)**2) / (f_low**2))) return np.fft.irfft(D * correction, n=len(dev)) # ═════════════════════════════════════════════════════════════════════════════ # SIGNAL PATH B — 400Hz reference square wave # ═════════════════════════════════════════════════════════════════════════════ def path_b(audio: np.ndarray, fs: float, f_low=F_MOD, f_high=F_CARRIER, mode='tracking') -> np.ndarray: """ Implements blocks [54]→[55/57]→[56] of Fig. 4. [54] LPF → pass 400Hz, remove 4kHz [55] Mode switch → 0° (tracking) or 90° (tracing) [57] 90° shifter → inserted in tracing mode [56] Axis crossing detector → 400Hz square wave mode='tracking' → square wave in phase with 400Hz record signal → chopper sees REAL component → VTA reading mode='tracing' → square wave shifted 90° → sees QUADRATURE component → tracing distortion reading """ # [54] LPF: remove carrier ref_sine = apply(_sos_lp(fs, f_low * 3, order=4), audio) # narrow bandpass to clean up the 400Hz tone ref_sine = apply(_sos_bp(fs, f_low, f_low*0.8, order=4), ref_sine) if mode == 'tracing': # [57] 90° phase shift — implemented via Hilbert transform analytic = sig.hilbert(ref_sine) ref_sine = np.imag(analytic) # 90° shifted # [56] Axis crossing detector: sine → square wave square = np.sign(ref_sine) square[square == 0] = 1.0 # avoid zeros at crossings return square # ═════════════════════════════════════════════════════════════════════════════ # CHOPPER [58] + LPF [60] # ═════════════════════════════════════════════════════════════════════════════ def chopper_lpf(deviation: np.ndarray, square: np.ndarray, fs: float, f_low=F_MOD) -> tuple: """ Block [58]: multiply deviation signal by reference square wave. When in-phase (tracking): output ≈ full-wave rectified sine → positive DC When quadrature (tracing): output is odd-symmetric → zero DC average Block [60]: LPF to extract DC proportional to angular mismatch. Returns: (dc_value, chopped_waveform) """ chopped = deviation * square # [60] LPF — cutoff well below f_low to average over many cycles # LPF cutoff at f_low/40: averages over ~40 cycles of the modulating # tone. This rejects wow (0.5-3 Hz beat products) and flutter beat # products while preserving the DC tracking/tracing component. dc_sig = apply(_sos_lp(fs, f_low / 40, order=4), chopped) dc_val = float(np.mean(dc_sig[len(dc_sig)//4:])) # skip filter transient return dc_val, chopped, dc_sig # ═════════════════════════════════════════════════════════════════════════════ # VTA COMPUTATION (patent formula, Col. 3-4) # ═════════════════════════════════════════════════════════════════════════════ def compute_vta(F_signed_hz: float, R: float = GROOVE_RADIUS_M, V_low: float = PEAK_VEL, theta_r: float = THETA_R_DEG, f_low: float = F_MOD) -> float: """ θP = arctan[ tan(θR) + C · R · F_signed / V_low ] where C = 9.1e-3 * (f_low / 400) F_signed carries the full sign of the VTA error: positive → pickup VTA steeper than recording angle negative → pickup VTA shallower than recording angle F_signed is derived from the chopper DC output sign × F_peak (unsigned), with the channel polarity (L inverts vertical, R does not) already applied. This matches exactly how the hardware voltmeter works: the chopper DC magnitude drives the meter deflection and its polarity (via channel switch) determines which side of centre the needle moves. Patent ref: Col 3-4, formula θP = arctan[tan(θR) ± 9.1e-3·R·F·cos(φ)/V400] The ± and cos(φ) together encode F_signed — we compute it directly from the chopper DC which the hardware meter already does implicitly. """ theta_R = np.radians(theta_r) C = 9.1e-3 * (f_low / 400.0) factor = C * R * F_signed_hz / V_low return float(np.degrees(np.arctan(np.tan(theta_R) + factor))) # ═════════════════════════════════════════════════════════════════════════════ # PEAK FM DEVIATION from DC output of chopper chain # ═════════════════════════════════════════════════════════════════════════════ def extract_F_and_phi(dev_signal: np.ndarray, audio: np.ndarray, fs: float, f_low=F_MOD) -> tuple: """ Extract peak FM deviation F (Hz) and phase φ (degrees). F is recovered from the RMS of the 400Hz band of the deviation signal. φ is the phase lead of E2 (deviation) w.r.t. E1 (400Hz on record). In the patent this is done implicitly by the chopper, but we also compute F and φ explicitly so we can apply the VTA formula. """ t = np.arange(len(dev_signal)) / fs # Narrow bandpass around f_low in deviation signal → E2 e2 = apply(_sos_bp(fs, f_low, f_low*0.6, order=4), dev_signal) # E1: 400Hz from record (path B before axis-crossing) e1_raw = apply(_sos_lp(fs, f_low*3, order=4), audio) e1 = apply(_sos_bp(fs, f_low, f_low*0.6, order=4), e1_raw) # Peak FM deviation and phase — computed per segment then averaged. # Segmenting over 0.5s windows means: # - Each segment spans 200 cycles of 400Hz — plenty for stable statistics # - Wow (0.5-3 Hz) varies slowly across segments; median rejects outliers # - Flutter residuals that slip through the BPF are averaged down # - Any single noisy segment (e.g. a tick or pop) is rejected by median seg_len = max(int(fs * 0.5), int(fs / f_low) * 20) # at least 20 cycles n_segs = max(1, len(e2) // seg_len) F_segs = [] phi_segs = [] for i in range(n_segs): sl = slice(i * seg_len, (i + 1) * seg_len) e2s = e2[sl] e1s = e1[sl] ts = t[sl] # F_peak for this segment F_segs.append(np.sqrt(2.0) * float(np.std(e2s))) # Phase for this segment via phasor correlation rc = np.cos(2*np.pi*f_low*ts) rs = np.sin(2*np.pi*f_low*ts) ph_e1 = np.angle(complex(np.dot(e1s, rc), np.dot(e1s, rs))) ph_e2 = np.angle(complex(np.dot(e2s, rc), np.dot(e2s, rs))) phi_segs.append(float(np.degrees(ph_e2 - ph_e1))) # Median across segments — robust against wow-induced outliers and clicks F_peak = float(np.median(F_segs)) # Circular median for phase (handles ±180° wrap correctly) phi_rad_segs = np.radians(phi_segs) phi_deg = float(np.degrees(np.arctan2( np.median(np.sin(phi_rad_segs)), np.median(np.cos(phi_rad_segs)) ))) # Normalise to (−180, +180] phi_deg = (phi_deg + 180) % 360 - 180 return F_peak, phi_deg # ═════════════════════════════════════════════════════════════════════════════ # TRACING DISTORTION MEASUREMENT # ═════════════════════════════════════════════════════════════════════════════ def tracing_distortion(dev_signal: np.ndarray, fs: float, f_low=F_MOD, n_harm=6) -> dict: """ Measure the amplitude of FM deviation harmonics at n·f_low. These arise from stylus tracing distortion (quadrature component). Returns {n: dB_relative_to_fundamental} """ N = len(dev_signal) freqs = np.fft.rfftfreq(N, 1/fs) spec = np.abs(np.fft.rfft(dev_signal * np.hanning(N))) * 2 / N def peak_bin(fc, bw=20.0): mask = (freqs >= fc-bw) & (freqs <= fc+bw) return float(np.max(spec[mask])) if np.any(mask) else 0.0 fund = peak_bin(f_low) out = {} for n in range(1, n_harm+1): a = peak_bin(f_low * n) out[n] = 20*np.log10(max(a, 1e-12) / max(fund, 1e-12)) return out # ═════════════════════════════════════════════════════════════════════════════ # FULL ANALYSIS PIPELINE # ═════════════════════════════════════════════════════════════════════════════ def analyse(audio: np.ndarray, fs: float, channel: str = 'L', f_low: float = F_MOD, f_high: float = F_CARRIER, R: float = GROOVE_RADIUS_M, V_low: float = PEAK_VEL, theta_r: float = THETA_R_DEG, swap_channels: bool = False, verbose: bool = True) -> dict: """Run the full patent signal chain on one channel of audio.""" # ── Input level check ──────────────────────────────────────────────── # The input MUST be raw cartridge velocity with NO RIAA applied. # RIAA boosts 400Hz ~+20dB and cuts 4kHz ~-12dB relative to flat, # destroying the 4:1 amplitude ratio the measurement depends on. # # Use NARROW bandpass filters (±40Hz and ±100Hz) to measure only the # actual tones, not broadband noise which would skew the ratio. # CBS record design: 400Hz tone is 4× (12 dB) the 4kHz tone in velocity. # With RIAA: ratio becomes ~44 dB — clearly detectable. # Track 1 (4kHz only): e_low will be near zero — skip the check. BW_LOW = 80.0 # Hz — narrow window around f_low tone (±40 Hz) BW_HIGH = 200.0 # Hz — narrow window around f_high tone (±100 Hz) e_low = float(np.std(apply(_sos_bp(fs, f_low, BW_LOW, order=4), audio))) e_high = float(np.std(apply(_sos_bp(fs, f_high, BW_HIGH, order=4), audio))) if e_high > 0 and e_low > 0: ratio_db = 20.0 * np.log10(e_low / e_high) if ratio_db > 28.0: # RIAA gives ~44dB; threshold at 28dB catches RIAA while # allowing for level variations across the different IMD bands warnings.warn( f"\n⚠ RIAA EQUALIZATION DETECTED " f"({f_low:.0f}Hz/{f_high/1000:.0f}kHz ratio = {ratio_db:.1f} dB)\n" f" Expected ~12 dB for flat/velocity input.\n" f" Use your preamp FLAT/BYPASS mode — do NOT apply RIAA.\n" f" VTA results will be unreliable with RIAA applied.", UserWarning, stacklevel=2) elif e_low < e_high: # Ratio is negative — 4kHz is louder than 400Hz # Suggests playing track 1 (4kHz reference only) or wrong track if verbose: print(f" [!] Level check: {f_low:.0f}Hz/{f_high/1000:.0f}kHz ratio = " f"{ratio_db:.1f} dB — is this the 4kHz-only reference track?") else: if verbose: status = "OK" if 8 <= ratio_db <= 20 else "CHECK" print(f" [{status}] Input level check: " f"{f_low:.0f}Hz/{f_high/1000:.0f}kHz ratio = {ratio_db:.1f} dB " f"(expected ~12 dB for flat input)") elif verbose: print(f" [!] Level check: could not detect tones clearly — " f"check input connections and volume") # ── Path A: FM discriminator ────────────────────────────────────────── dev_signal = path_a(audio, fs, f_low, f_high) # ── Chopper chain — TRACKING mode (0°) ─────────────────────────────── sq_track = path_b(audio, fs, f_low, f_high, mode='tracking') dc_track, chopped_track, dc_track_sig = chopper_lpf(dev_signal, sq_track, fs, f_low) # ── Chopper chain — TRACING mode (90°) ─────────────────────────────── sq_trace = path_b(audio, fs, f_low, f_high, mode='tracing') dc_trace, chopped_trace, dc_trace_sig = chopper_lpf(dev_signal, sq_trace, fs, f_low) # ── Extract F (peak deviation, unsigned) ───────────────────────────── F_peak, phi_deg = extract_F_and_phi(dev_signal, audio, fs, f_low) # ── F_signed: sign from most stable segments ────────────────────────── # Split into 2-second segments, rank by |DC| magnitude, take mean of # top 50% — most settled/signal-rich segments determine the sign. # Falls back to overall dc_track if signal is too short for segmenting. seg_len_sign = int(fs * 2.0) dc_segs = [] dc_top_mean = 0.0 # initialise here so verbose print never crashes # Only segment if we have at least 2 full segments if len(audio) >= seg_len_sign * 2: n_segs_sign = len(audio) // seg_len_sign for i in range(n_segs_sign): sl = slice(i * seg_len_sign, (i + 1) * seg_len_sign) seg = audio[sl] if len(seg) < seg_len_sign: continue dev_s = path_a(seg, fs, f_low, f_high) sq_s = path_b(seg, fs, f_low, f_high, mode='tracking') dc_s, _, _ = chopper_lpf(dev_s, sq_s, fs, f_low) if abs(dc_s) > 0.5: dc_segs.append(dc_s) if dc_segs: dc_sorted = sorted(dc_segs, key=abs, reverse=True) top_half = dc_sorted[:max(1, len(dc_sorted) // 2)] dc_top_mean = float(np.mean(top_half)) dc_sign_robust = np.sign(dc_top_mean) else: # Short signal or all segments near zero — use overall dc_track directly dc_top_mean = dc_track dc_sign_robust = np.sign(dc_track) # Channel polarity: left channel inverts vertical modulation in groove # SWAP_CHANNELS flips this if the recording has L/R wired in reverse ch_is_left = (channel.upper() == 'L') if swap_channels: ch_is_left = not ch_is_left ch_sign = -1.0 if ch_is_left else +1.0 F_signed = ch_sign * dc_sign_robust * F_peak if verbose: n_total = len(dc_segs) n_correct = sum(1 for d in dc_segs if np.sign(d) == dc_sign_robust) top_n = max(1, n_total // 2) if n_total == 0: print(f" [Sign] using overall DC={dc_top_mean:+.1f} " f"(signal too short for segment analysis)") else: stability = 'stable' if n_total == 0 or n_correct/n_total > 0.6 \ else 'check recording' print(f" [Sign] top-{top_n} segments mean DC={dc_top_mean:+.1f} " f"{n_correct}/{n_total} segments agree ({stability})") # ── VTA computation (patent formula) ───────────────────────────────── vta = compute_vta(F_signed, R, V_low, theta_r, f_low) vta_err = vta - theta_r # ── Tracing distortion harmonics ───────────────────────────────────── harmonics = tracing_distortion(dev_signal, fs, f_low) # ── DC tracking signal (offset-corrected, as meter [62]) ───────────── # offset sets zero-deviation point at θR (16.0° for CBS STR-112) meter_reading = vta # in degrees, directly results = dict( channel = channel, f_low = f_low, f_high = f_high, fs = fs, F_peak_hz = F_peak, F_signed_hz = F_signed, phi_deg = phi_deg, vta_deg = vta, vta_error_deg = vta_err, theta_r_deg = theta_r, dc_tracking = dc_track, dc_tracing = dc_trace, harmonics_db = harmonics, dev_signal = dev_signal, sq_track = sq_track, chopped_track = chopped_track, dc_track_sig = dc_track_sig, chopped_trace = chopped_trace, dc_trace_sig = dc_trace_sig, audio = audio, ) if verbose: bar_w = 30 print(f"\n╔══════════════════════════════════════════════════╗") print(f"║ VTA ANALYSIS — Channel {channel} ({f_low:.0f}/{f_high:.0f} Hz) ║") print(f"╠══════════════════════════════════════════════════╣") print(f"║ Peak FM deviation (F) : {F_peak:8.2f} Hz ║") print(f"║ Signed FM deviation : {F_signed:+8.2f} Hz ║") print(f"║ Phase lead E2/E1 (φ) : {phi_deg:+8.1f}° ║") print(f"║ Vertical Tracking Angle : {vta:8.2f}° ║") print(f"║ VTA Error (vs θR={theta_r:.1f}°) : {vta_err:+8.2f}° ║") print(f"║ DC tracking component : {dc_track:+8.4f} ║") print(f"║ DC tracing component : {dc_trace:+8.4f} ║") print(f"╠══════════════════════════════════════════════════╣") print(f"║ Wow/flutter tolerance: ║") print(f"║ Wow (<3Hz) and flutter (>10Hz offset) are fully ║") print(f"║ rejected by the 400Hz BPF and chopper averaging.║") print(f"║ Flutter within ~10Hz of {f_low:.0f}Hz adds <0.5deg error. ║") print(f"╠══════════════════════════════════════════════════╣") print(f"║ Tracing distortion harmonics (dB re fundamental)║") for n, db in harmonics.items(): filled = int(max(0, (db + 60)/2)) bar = '█'*filled + '░'*(bar_w-filled) print(f"║ {n:2d}×{f_low:.0f}Hz : {db:+6.1f} dB {bar[:20]} ║") print(f"╚══════════════════════════════════════════════════╝") return results # ═════════════════════════════════════════════════════════════════════════════ # SYNTHETIC TEST SIGNAL GENERATOR (calibration / demo) # ═════════════════════════════════════════════════════════════════════════════ def generate_test_signal(fs: float = 44100.0, duration: float = 10.0, vta_error_deg: float = 5.0, channel: str = 'L', f_low: float = F_MOD, f_high: float = F_CARRIER, tracing_frac: float = 0.02) -> np.ndarray: """ Synthesize a CBS STR-112-style test signal with a known VTA error injected. Also adds a small quadrature (tracing) component. VTA error is encoded as FM deviation F on the f_high carrier: F = V_low * (tan(thetaP) - tan(thetaR)) / (C * R * cos(0)) where C = 9.1e-3*(f_low/400) (left channel → − sign, so we negate) """ t = np.arange(int(fs * duration)) / fs # Stereo groove geometry for VERTICAL modulation: # The ch_sign in the analyser is -1 for L, +1 for R. # To produce the correct DC polarity from the chopper, the FM deviation # in the generated signal must be inverted for L relative to R. # ch_sign=-1 for L means the analyser expects NEGATIVE dc for a positive VTA error. # So the L channel groove must produce NEGATIVE FM deviation → groove_sign = -1 for L. groove_sign = -1.0 if channel.upper() == 'L' else +1.0 # Compute signed FM deviation for target VTA error. # grove_sign=-1 for L, +1 for R models stereo groove geometry. # F_inj must keep its sign (positive=VTA too steep, negative=too shallow). theta_P = np.radians(THETA_R_DEG + vta_error_deg) theta_R = np.radians(THETA_R_DEG) C_gen = 9.1e-3 * (f_low / 400.0) F_inj = (PEAK_VEL * (np.tan(theta_P) - np.tan(theta_R))) \ / (C_gen * GROOVE_RADIUS_M) # signed — keeps direction of error # 400Hz modulating tone (large amplitude — E1) mod_amp = 0.7 mod_tone = mod_amp * np.sin(2*np.pi*f_low*t) # 4kHz carrier FM-modulated by VTA error (tracking) and tracing components # groove_sign models stereo geometry: L and R see opposite FM phase mod_idx_track = groove_sign * F_inj / f_low mod_idx_trace = tracing_frac * abs(mod_idx_track) carrier_phase = (2*np.pi*f_high*t + mod_idx_track * np.sin(2*np.pi*f_low*t) # tracking (VTA error) + mod_idx_trace * (-np.cos(2*np.pi*f_low*t))) # tracing (quadrature) # Note: tracing_frac kept very small (default 0.02) so it does not # overwhelm the tracking signal in the chopper DC sign detection. carrier = 0.15 * np.cos(carrier_phase) # 4kHz much smaller amplitude combined = mod_tone + carrier # Normalise pk = np.max(np.abs(combined)) if pk > 0: combined /= pk * 1.05 return combined.astype(np.float32) # ═════════════════════════════════════════════════════════════════════════════ # PLOTTING — replicates meter face + signal chain visualisation # ═════════════════════════════════════════════════════════════════════════════ DARK = '#0e1117' MID = '#1c2230' ACCENT = '#00d4ff' GREEN = '#39ff14' AMBER = '#ffb300' RED = '#ff3d3d' WHITE = '#e8eaf0' GREY = '#5a6070' def _style_ax(ax, title=''): ax.set_facecolor(MID) for sp in ax.spines.values(): sp.set_color(WHITE) ax.tick_params(colors=WHITE, labelsize=10) ax.xaxis.label.set_color(WHITE) ax.yaxis.label.set_color(WHITE) if title: ax.set_title(title, color=WHITE, fontsize=11, pad=6) def plot_signal_chain(res: dict): """Show each stage of the patent signal chain for one channel.""" fs = res['fs'] show = min(int(fs * 0.025), len(res['audio'])) # 25 ms t = np.arange(show) / fs * 1000 # ms fig = plt.figure(figsize=(14, 9), facecolor=DARK) fig.suptitle(f"Patent Signal Chain — Channel {res['channel']} " f"({res['f_low']:.0f}/{res['f_high']:.0f} Hz)", color=WHITE, fontsize=15, fontweight='bold') stages = [ ('Input (audio)', res['audio'][:show], ACCENT), ('Path A: FM deviation', res['dev_signal'][:show], AMBER), ('Path B: square wave', res['sq_track'][:show], GREEN), ('Chopper out (track)', res['chopped_track'][:show], '#ff9800'), ('LPF → DC (tracking)', res['dc_track_sig'][:show], GREEN), ('Chopper out (trace)', res['chopped_trace'][:show], '#e040fb'), ('LPF → DC (tracing)', res['dc_trace_sig'][:show], '#e040fb'), ] for i, (label, data, col) in enumerate(stages): ax = fig.add_subplot(len(stages), 1, i+1) ax.plot(t, data, color=col, linewidth=0.9, alpha=0.92) ax.axhline(0, color=WHITE, linewidth=0.4, linestyle='--', alpha=0.3) _style_ax(ax, label) ax.set_xlim(0, t[-1]) if i == len(stages)-1: ax.set_xlabel('Time (ms)', color=WHITE, fontsize=10) plt.tight_layout(rect=[0, 0, 1, 0.96]) return fig def plot_meter_dashboard(results_list: list): """Meter face (Fig. 5) + harmonics for all analysed results.""" n = len(results_list) fig = plt.figure(figsize=(5*n, 8), facecolor=DARK) fig.suptitle('Vertical Tracking Angle Meter — US Patent 4,359,768', color=WHITE, fontsize=16, fontweight='bold') for col_idx, res in enumerate(results_list): gs = gridspec.GridSpec(2, n, figure=fig, left=0.05, right=0.97, top=0.88, bottom=0.08, hspace=0.45, wspace=0.4) # ── Meter face (top) ────────────────────────────────────────────── ax_m = fig.add_subplot(gs[0, col_idx], projection='polar') ax_m.set_facecolor(MID) ax_m.set_theta_direction(-1) ax_m.set_theta_offset(np.pi) # 0 at top-left ax_m.set_ylim(0, 1.15) ax_m.set_yticks([]) # Scale: 5° to 35° mapped to 0°→180° sweep vta_min, vta_max = 5.0, 40.0 vta_range = vta_max - vta_min def vta_to_rad(v): return np.pi * (v - vta_min) / vta_range # Coloured arc zones zones = [(5, 13, RED), (13, 17, AMBER), (17, 25, GREEN), (25, 33, AMBER), (33, 40, RED)] for z0, z1, zcol in zones: th = np.linspace(vta_to_rad(z0), vta_to_rad(z1), 60) ax_m.plot(th, np.ones(60)*1.0, color=zcol, linewidth=8, solid_capstyle='butt', alpha=0.7) # Tick marks and labels for deg in range(5, 41, 5): th = vta_to_rad(deg) ax_m.plot([th, th], [0.85, 1.0], color=WHITE, linewidth=0.8) ax_m.text(th, 1.12, f'{deg}°', ha='center', va='center', color=WHITE, fontsize=10, rotation=np.degrees(th) - 90) # Recording angle marker th_r = vta_to_rad(res['theta_r_deg']) ax_m.plot([th_r, th_r], [0.75, 1.03], color=ACCENT, linewidth=1.2, linestyle='--', alpha=0.7) ax_m.text(th_r, 0.65, f"θR={res['theta_r_deg']}°", ha='center', color=ACCENT, fontsize=10) # Needle vta_clamped = np.clip(res['vta_deg'], vta_min, vta_max) th_n = vta_to_rad(vta_clamped) needle_col = GREEN if abs(res['vta_error_deg']) < 2 else \ AMBER if abs(res['vta_error_deg']) < 5 else RED ax_m.annotate('', xy=(th_n, 0.92), xytext=(th_n, 0.0), arrowprops=dict(arrowstyle='->', color=needle_col, lw=2.5, mutation_scale=12)) ax_m.plot(0, 0, 'o', color=WHITE, markersize=5, zorder=10) ax_m.set_title( f"Ch {res['channel']} {res['f_low']:.0f}/{res['f_high']:.0f} Hz\n" f"VTA = {res['vta_deg']:.1f}° (err {res['vta_error_deg']:+.1f}°)\n" f"F={res['F_peak_hz']:.1f} Hz φ={res['phi_deg']:+.0f}°", color=WHITE, fontsize=11, pad=14) # Suppress polar grid clutter ax_m.set_rticks([]) ax_m.set_thetagrids([]) ax_m.spines['polar'].set_visible(False) # ── Tracing distortion harmonics (bottom) ───────────────────────── ax_h = fig.add_subplot(gs[1, col_idx]) ax_h.set_facecolor(MID) harms = list(res['harmonics_db'].keys()) dbvals = [res['harmonics_db'][h] for h in harms] cmap = plt.cm.plasma colors = [cmap(h/max(harms)) for h in harms] ax_h.bar(harms, dbvals, color=colors, edgecolor=DARK, width=0.6) ax_h.axhline(-20, color=AMBER, linewidth=0.8, linestyle='--', alpha=0.7, label='−20 dB') ax_h.axhline(-40, color=GREEN, linewidth=0.8, linestyle='--', alpha=0.7, label='−40 dB') ax_h.set_xlabel('Harmonic', color=WHITE, fontsize=11) ax_h.set_ylabel('dB re fundamental', color=WHITE, fontsize=11) ax_h.set_xticks(harms) ax_h.set_xticklabels([f'{n}×' for n in harms]) ax_h.legend(fontsize=10, facecolor=MID, edgecolor=MID, labelcolor=WHITE, loc='lower right') _style_ax(ax_h, f'Tracing Distortion — {res["f_low"]:.0f}Hz harmonics') return fig # ═════════════════════════════════════════════════════════════════════════════ # SPECTRUM PLOT — verify tone detection and signal quality # ═════════════════════════════════════════════════════════════════════════════ def plot_spectrum(audio_map: dict, fs: float, f_low: float = F_MOD, f_high: float = F_CARRIER) -> plt.Figure: """ FFT spectrum for all channels (L and/or R) shown side by side. Top row : full spectrum 0 – f_high*1.5, tones annotated. Bottom row : zoom around the carrier showing FM sidebands. audio_map : dict {'L': array, 'R': array} — one or both channels. """ channels = list(audio_map.keys()) n_ch = len(channels) ch_colors = {'L': ACCENT, 'R': '#ff6ec7'} # cyan for L, pink for R fig, axes = plt.subplots(2, n_ch, figsize=(7 * n_ch, 9), facecolor=DARK, squeeze=False) fig.suptitle(f'Input Spectrum — {f_low:.0f} Hz / {f_high:.0f} Hz', color=WHITE, fontsize=15, fontweight='bold') zoom_bw = max(f_low * 4, 600) # ±zoom_bw Hz around carrier f_max = min(f_high * 1.6, fs / 2 - 1) for col, ch in enumerate(channels): audio = audio_map[ch] col_color = ch_colors.get(ch, ACCENT) N = len(audio) freqs = np.fft.rfftfreq(N, 1 / fs) spec = 20 * np.log10( np.abs(np.fft.rfft(audio * np.hanning(N))) * 2 / N + 1e-12) def find_peak_db(fc, bw=60): m = (freqs >= fc - bw) & (freqs <= fc + bw) return float(np.max(spec[m])) if np.any(m) else -120.0 # ── Top: full spectrum ──────────────────────────────────────────── ax0 = axes[0, col] ax0.set_facecolor(MID) mask = freqs <= f_max ax0.plot(freqs[mask], spec[mask], color=col_color, linewidth=0.8, alpha=0.92) ax0.axvline(f_low, color=GREEN, linewidth=1.2, linestyle='--', label=f'{f_low:.0f} Hz (mod tone)') ax0.axvline(f_high, color=AMBER, linewidth=1.2, linestyle='--', label=f'{f_high:.0f} Hz (carrier)') ax0.set_xlabel('Frequency (Hz)', color=WHITE, fontsize=11) ax0.set_ylabel('Amplitude (dBFS)', color=WHITE, fontsize=11) ax0.legend(fontsize=10, facecolor=MID, edgecolor=MID, labelcolor=WHITE) _style_ax(ax0, f'Channel {ch} — Full Spectrum') # Annotate tone peak levels and ratio db_low = find_peak_db(f_low) db_high = find_peak_db(f_high) ratio = db_low - db_high # Place annotation to the right of f_low marker x_ann = f_low + (f_high - f_low) * 0.08 ax0.annotate(f'{db_low:.1f} dBFS', xy=(f_low, db_low), xytext=(x_ann, db_low - 4), color=GREEN, fontsize=10, arrowprops=dict(arrowstyle='->', color=GREEN, lw=1.0)) ax0.annotate(f'{db_high:.1f} dBFS\nratio {ratio:+.1f} dB', xy=(f_high, db_high), xytext=(f_high * 1.03, db_high - 10), color=AMBER, fontsize=10, arrowprops=dict(arrowstyle='->', color=AMBER, lw=1.0)) # Level check verdict inline on plot if 8 <= ratio <= 20: verdict, vcol = f'✓ flat input ({ratio:.1f} dB)', GREEN elif ratio > 28: verdict, vcol = f'⚠ RIAA detected! ({ratio:.1f} dB)', RED else: verdict, vcol = f'? check input ({ratio:.1f} dB)', AMBER ax0.text(0.02, 0.04, verdict, transform=ax0.transAxes, color=vcol, fontsize=10, fontweight='bold', va='bottom', ha='left') # ── Bottom: carrier zoom — FM sidebands ────────────────────────── ax1 = axes[1, col] ax1.set_facecolor(MID) mask2 = (freqs >= f_high - zoom_bw) & (freqs <= f_high + zoom_bw) ax1.plot(freqs[mask2], spec[mask2], color=col_color, linewidth=0.9, alpha=0.92) ax1.axvline(f_high, color=WHITE, linewidth=1.0, linestyle='--', alpha=0.6, label='Carrier') # Expected sideband positions for n in [1, 2, 3]: for sgn in [-1, +1]: sb = f_high + sgn * n * f_low if f_high - zoom_bw <= sb <= f_high + zoom_bw: lbl = (f'±{n}×{f_low:.0f} Hz sidebands' if n == 1 and sgn == 1 else '') ax1.axvline(sb, color=GREEN, linewidth=0.8, linestyle=':', alpha=0.7, label=lbl) # label the sideband order at top of plot ax1.text(sb, ax1.get_ylim()[1] if ax1.get_ylim()[1] != 0 else -10, f'{n}×', color=GREEN, fontsize=9, ha='center', va='bottom', alpha=0.8) ax1.set_xlabel('Frequency (Hz)', color=WHITE, fontsize=11) ax1.set_ylabel('Amplitude (dBFS)', color=WHITE, fontsize=11) ax1.legend(fontsize=10, facecolor=MID, edgecolor=MID, labelcolor=WHITE) _style_ax(ax1, f'Channel {ch} — Carrier zoom (sidebands = FM deviation)') for ax in [ax0, ax1]: for sp in ax.spines.values(): sp.set_color(WHITE) ax.tick_params(colors=WHITE, labelsize=10) plt.tight_layout(rect=[0, 0, 1, 0.95]) return fig # ═════════════════════════════════════════════════════════════════════════════ # SUMMARY REPORT — averaged result across L and R channels # ═════════════════════════════════════════════════════════════════════════════ def print_summary(results_list: list): """ Print a final summary. If both L and R channels were measured, average them as the patent recommends (Col. 7): 'it is acceptable to average the vertical angle readings obtained on the left and right channels to obtain the final measurement.' Also give a tonearm adjustment recommendation. """ print() print("╔══════════════════════════════════════════════════════╗") print("║ FINAL SUMMARY REPORT ║") print("╠══════════════════════════════════════════════════════╣") vta_values = [r['vta_deg'] for r in results_list] theta_r = results_list[0]['theta_r_deg'] for r in results_list: print(f"║ Channel {r['channel']}: VTA = {r['vta_deg']:6.2f}° " f"(error {r['vta_error_deg']:+.2f}°) ║") if len(vta_values) == 2: vta_avg = float(np.mean(vta_values)) err_avg = vta_avg - theta_r print(f"╠══════════════════════════════════════════════════════╣") print(f"║ Average VTA : {vta_avg:6.2f}° (error {err_avg:+.2f}°) ║") else: err_avg = results_list[0]['vta_error_deg'] vta_avg = results_list[0]['vta_deg'] print(f"╠══════════════════════════════════════════════════════╣") print(f"║ Recording angle θR : {theta_r:.1f}° ║") print(f"║ Measured VTA : {vta_avg:.2f}° ║") abs_err = abs(err_avg) if abs_err < 0.5: verdict = "EXCELLENT — VTA is within ±0.5° of recording angle" action = "No adjustment needed." elif abs_err < 2.0: verdict = "GOOD — VTA is within ±2°" action = "Minor adjustment beneficial but optional." elif abs_err < 5.0: verdict = "FAIR — VTA error is significant" action = "Adjustment recommended." else: verdict = "POOR — VTA error is large, distortion likely" action = "Adjustment required." print(f"╠══════════════════════════════════════════════════════╣") print(f"║ {verdict[:52]:<52}║") print(f"╠══════════════════════════════════════════════════════╣") if err_avg > 0.5: print(f"║ ADJUSTMENT: VTA too STEEP by {err_avg:+.1f}° ║") print(f"║ → Lower the rear of the tonearm (or raise front) ║") print(f"║ until meter reads {theta_r:.1f}° ║") elif err_avg < -0.5: print(f"║ ADJUSTMENT: VTA too SHALLOW by {abs(err_avg):.1f}° ║") print(f"║ → Raise the rear of the tonearm (or lower front) ║") print(f"║ until meter reads {theta_r:.1f}° ║") else: print(f"║ {action:<52}║") # Tracing distortion summary print(f"╠══════════════════════════════════════════════════════╣") print(f"║ Tracing distortion (2nd harmonic): ║") for r in results_list: h2 = r['harmonics_db'].get(2, -120) td_verdict = ("low" if h2 < -40 else "moderate" if h2 < -20 else "HIGH — check stylus condition") print(f"║ Ch {r['channel']}: {h2:+6.1f} dB re fundamental ({td_verdict:<20})║") print(f"╚══════════════════════════════════════════════════════╝") # ═════════════════════════════════════════════════════════════════════════════ # SPYDER / INTERACTIVE ENTRY POINT # Edit the settings below then press F5 (Run file) in Spyder # ═════════════════════════════════════════════════════════════════════════════ if __name__ == '__main__': # ── Self-test: verify sign handling is correct ────────────────── _audio_test = generate_test_signal(fs=44100., vta_error_deg=5., channel='L', f_low=400., f_high=4000.) _r = analyse(_audio_test, 44100., channel='L', f_low=400., f_high=4000., R=GROOVE_RADIUS_M, V_low=0.3299, theta_r=16.0, verbose=False) _vta_ok = 18.0 < _r['vta_deg'] < 24.0 # should be ~21 deg if _vta_ok: print('VTA Analyzer v2.0 — self-test PASSED ✓ Left channel sign is correct') else: print(f'*** WARNING: self-test FAILED — Left ch VTA={_r["vta_deg"]:.1f}° (expected ~21°)') print('*** Old code may be cached. Close Spyder completely, reopen, then F5.') print() # Settings are defined at the TOP of this file — scroll up to edit them. # (Look for the "USER SETTINGS" banner near line 70) results_all = [] # ── Load or generate audio ──────────────────────────────────────────── if DEMO_MODE: print(f"\n[DEMO MODE] Synthetic CBS STR-112 signal") print(f" f_low={F_MOD:.0f} Hz, injected VTA error = {DEMO_VTA_ERROR:+.1f}° " f"(θP = {RECORDING_ANGLE + DEMO_VTA_ERROR:.1f}°)") channels = ['L', 'R'] if CHANNEL == 'both' else [CHANNEL] fs = 44100.0 audio_map = {ch: generate_test_signal(fs=fs, vta_error_deg=DEMO_VTA_ERROR, channel=ch, f_low=F_MOD, f_high=F_CARRIER) for ch in channels} else: if not _HAS_SF: raise ImportError( "soundfile is not installed.\n" "In the Spyder IPython console run: pip install soundfile\n" "Then restart the kernel (circular arrow button).") import soundfile as sf data, fs = sf.read(AUDIO_FILE, always_2d=True) pk = np.max(np.abs(data)) if pk > 0: data = data / pk print(f"\nLoaded: {AUDIO_FILE}") print(f" Sample rate : {fs:.0f} Hz") print(f" Duration : {data.shape[0]/fs:.1f} s") print(f" Channels : {data.shape[1]}") if CHANNEL == 'both' and data.shape[1] >= 2: channels = ['L', 'R'] else: channels = [CHANNEL if CHANNEL != 'both' else 'L'] audio_map = { ch: data[:, (0 if ch == 'L' else min(1, data.shape[1]-1))] for ch in channels } # ── Spectrum check (before analysis — verify tones are present) ─────── if SHOW_SPECTRUM: fig_spec = plot_spectrum(audio_map, fs, f_low=F_MOD, f_high=F_CARRIER) plt.show(block=False) # ── Run analysis on each channel ────────────────────────────────────── for ch, audio in audio_map.items(): res = analyse(audio, fs, channel = ch, f_low = F_MOD, f_high = F_CARRIER, R = GROOVE_RADIUS, V_low = PEAK_VEL, theta_r = RECORDING_ANGLE, swap_channels = SWAP_CHANNELS, verbose = True) results_all.append(res) # ── Summary report ──────────────────────────────────────────────────── print_summary(results_all) # ── Signal chain and meter plots ────────────────────────────────────── if SHOW_SIGNAL_CHAIN: fig_chain = plot_signal_chain(results_all[0]) plt.show(block=False) if SHOW_METER: fig_meter = plot_meter_dashboard(results_all) plt.show(block=False) plt.show() # Block here so all windows stay open