@@ -320,3 +320,95 @@ def estimate_cfo(iq, sample_rate, nfft=None):
320320 spec = np .abs (np .fft .fftshift (np .fft .fft (iq [:nfft ], nfft )))
321321 freqs = np .fft .fftshift (np .fft .fftfreq (nfft , 1.0 / sample_rate ))
322322 return float (freqs [int (np .argmax (spec ))])
323+
324+
325+ def estimate_fm_cfo (iq , sample_rate , * , channel_bw = 100_000.0 , trim_hz = None ):
326+ """Estimate a broadcast-FM capture's carrier frequency offset, in Hz.
327+
328+ Method: the phase discriminator's DC term. Under a constant carrier
329+ offset df, ``angle(x[n] * conj(x[n-1]))`` picks up that same df at every
330+ sample, while the program audio (the wanted FM message) averages to zero
331+ over a long window. So the MEAN of the demodulated instantaneous frequency
332+ IS the carrier offset. This is the fine, accurate estimator for FM, and is
333+ distinct from ``estimate_cfo``: that one finds where energy sits in the
334+ passband (a coarse, spread-out measurement for a wideband FM signal),
335+ whereas this recovers the sub-kHz crystal offset a real receiver must
336+ correct.
337+
338+ The capture is channelized (lowpass + decimate to the FM channel) before
339+ discriminating, for exactly the reason ``fm_pilot_excess_db`` documents: a
340+ discriminator fed the full capture bandwidth is dominated by out-of-channel
341+ noise and adjacent stations, and the mean of that is meaningless. On the
342+ channelized signal the mean is the offset.
343+
344+ Why NOT the 19 kHz pilot: with a differential discriminator a carrier
345+ offset adds a DC term but does NOT move the pilot -- differentiation is
346+ invariant to a constant frequency shift -- so the pilot's position carries
347+ no CFO information (this is the trap the TODO's "or the pilot frequency"
348+ wording invites). The pilot's real role here is a presence gate: confirm a
349+ station is actually there with ``fm_pilot_excess_db`` before trusting this
350+ number, because the mean of noise is not a carrier offset.
351+
352+ Args:
353+ channel_bw: FM channel half-not-needed; used only to size the click
354+ trim. Broadcast FM is ~100 kHz half-bandwidth of program content.
355+ trim_hz: discriminator "clicks" (phase wraps in deep fades / low SNR)
356+ are impulsive outliers that bias a plain mean. Samples whose
357+ demodulated frequency exceeds trim_hz in magnitude are dropped
358+ before averaging. Default 1.5 * channel_bw keeps the full +/-75 kHz
359+ deviation swing while rejecting the near-+/-Nyquist click spikes.
360+ Pass None to disable trimming.
361+
362+ Returns the offset in Hz (positive = the signal sits above the tuned
363+ center), or None when the capture is too short or the rate too low to
364+ channelize to the FM channel.
365+
366+ OUR code.
367+ """
368+ iq = np .asarray (iq , dtype = np .complex64 )
369+ if sample_rate <= 4.0 * channel_bw or iq .size < 4096 :
370+ return None
371+ # Channelize to the FM channel, mirroring fm_pilot_excess_db: lowpass and
372+ # decimate to ~500 kSps before discriminating so out-of-channel energy
373+ # does not pollute the mean. Taps scale with the decimation.
374+ decim = max (1 , int (sample_rate // 480_000 ))
375+ if decim > 1 :
376+ from .filters import design_lowpass , fir_apply
377+ taps = design_lowpass (120_000.0 , sample_rate , num_taps = 8 * decim + 1 )
378+ iq = fir_apply (iq , taps )[::decim ]
379+ sample_rate = sample_rate / decim
380+ if iq .size < 2 :
381+ return None
382+ # phase discriminator -> instantaneous frequency in Hz
383+ inst_hz = (np .angle (iq [1 :] * np .conj (iq [:- 1 ])).astype (np .float64 )
384+ * sample_rate / (2.0 * np .pi ))
385+ trim = 1.5 * float (channel_bw ) if trim_hz is None else float (trim_hz )
386+ if trim > 0 :
387+ keep = np .abs (inst_hz ) <= trim
388+ if keep .any ():
389+ inst_hz = inst_hz [keep ]
390+ if inst_hz .size == 0 :
391+ return None
392+ return float (np .mean (inst_hz ))
393+
394+
395+ def correct_fm_cfo (iq , sample_rate , cfo_hz = None , ** kwargs ):
396+ """Remove a broadcast-FM carrier frequency offset. OUR code.
397+
398+ If ``cfo_hz`` is None it is measured with ``estimate_fm_cfo`` (any keyword
399+ arguments are forwarded to it). The correction is a frequency shift by
400+ ``-cfo_hz`` -- ``frequency_shift(iq, -cfo_hz, sample_rate)`` -- which brings
401+ the carrier back to center. Estimating and correcting are separated on
402+ purpose (``estimate_fm_cfo`` never mutates data), and this is the
403+ convenience that does both in one call for the common case.
404+
405+ Returns ``(corrected_iq, cfo_hz)`` so the applied correction is visible to
406+ the caller. When the offset cannot be estimated (capture too short or rate
407+ too low), returns the input unchanged with ``cfo_hz = 0.0``.
408+ """
409+ if cfo_hz is None :
410+ cfo_hz = estimate_fm_cfo (iq , sample_rate , ** kwargs )
411+ if cfo_hz is None :
412+ return np .asarray (iq , dtype = np .complex64 ), 0.0
413+ from .mixing import frequency_shift
414+ return frequency_shift (iq , - float (cfo_hz ), sample_rate ), float (cfo_hz )
0 commit comments