doppler_spread.m realizes ~0.71x the specified Doppler spread

[spinkham]

Summary

octave/doppler_spread.m hands the Gaussian Doppler power spectrum directly
to fir2 as the filter's amplitude target:

sigma = dopplerSpreadHz/2;
x = 0:lowFs/100:lowFs/2;
y = (1/(sigma*sqrt(2*pi)))*exp(-(x.^2)/(2*sigma*sigma));
b = fir2(Ntaps-1, x/(lowFs/2), y);

Since the output PSD of filtered white noise is |H(f)|^2, the realized
fading-process power spectrum is the square of that Gaussian — a Gaussian
with standard deviation sigma/sqrt(2). The realized 2-sigma Doppler spread
is therefore 1/sqrt(2) of the specified value. Measured against current
main with the script below: asking for a 1.0 Hz spread yields 0.71 Hz
(PSD second moment 0.715, autocorrelation fit 0.709 — both right on the
1/sqrt(2) = 0.707 theory). All copies I checked carry the same line:
codec2/octave/doppler_spread.m, the legacy codec2-dev history back to
its 2016 introduction, misc/freedv_low/doppler_spread.m (the version
issues #16/#18 reference as bug-fixed — those fixes address the resampler,
not the spectral shape), and radae/doppler_spread.m.

The existing unit test (doppler_spread_ut.m) compares the designed
filter's amplitude response against y, so it passes by construction
and cannot see this.

Why the spread is defined on the power spectrum

Every standard and implementation I could check specifies the Gaussian and
its 2-sigma "frequency spread" / "Doppler spread" on the power spectrum
of the tap-gain process, which means a shaping filter needs
|H| = sqrt(PSD):

• ITU-R F.1487 (Annex 1, eq. 2): "Each tap-gain function has a power
  spectrum, f_i(ν) … each of which is a Gaussian function of frequency",
  and (Annex 3) "the frequency spread is the 2σ value as used in
  equation (2)".
• MIL-STD-188-110C, Appendix E (§E.5.4) writes the filter requirement
  with the square explicit, |H_j(f)|² = e^(−2f²/d_j²)/sqrt(π d_j²/2),
  and gives the corresponding time-domain taps
  f_j(t) = k·sqrt(2)·e^(−π²t²d²) (truncated where the tap falls to 1% of
  peak).
  https://everyspec.com/MIL-STD/MIL-STD-0100-0299/download.php?spec=MIL-STD-188_110C.037889.PDF
  (pp. 218–226)
• NTIA/ITS "HF Simulator" manual (E. Johnson, 1991, ITS-sponsored,
  citing Watterson et al. 1970): random processes "are produced by passing
  Rayleigh-distributed noise through a filter whose |H(jΩ)|² has the
  required Gaussian shape."
  https://its.ntia.gov/umbraco/surface/download/publication?reportNumber=HF+Simulator.pdf
• PathSim — the reference the code comment credits ("Used PathSim
  technical guide as a reference - thanks Moe!") — is also on this side,
  though its prose is easy to misread: it builds a Gaussian impulse
  response with σ_taps = Fs·sqrt(2)/(2π·F2σ) (pathsimtech100 §4.1.5),
  which gives an amplitude response of width F2σ/sqrt(2), exactly
  sqrt(PSD) of a Gaussian PSD with 2-sigma width F2σ. The sqrt lives
  in the width constant. (I've sent the modernized PathSim port a doc PR
  spelling this out at the constant, so the guide's prose stops catching
  people: docs(GaussFIR): explain the sqrt(2) in the Doppler filter width constant bubnikv/pathsim#4.)

Reproduction (Octave, run from codec2/octave)

Two independent estimators, both needed: the FFT bins must be much finer
than the spectrum (coarse bins smear the second moment upward), and the
autocorrelation fit is window-free.

pkg load signal
d = 1.0; Fs = 8000; N = Fs*3600;                % ask for 1.0 Hz spread
s = doppler_spread(d, Fs, N);
% (a) PSD second moment, fine bins (2^19 -> 0.015 Hz/bin)
seg = 2^19; nseg = floor(length(s)/seg); P = zeros(seg,1);
for k=1:nseg
  xk = s((k-1)*seg+1:k*seg)(:) .* hanning(seg);
  P = P + abs(fft(xk)).^2;
end
f = ((0:seg-1)'/seg)*Fs; f(f>Fs/2) -= Fs;
m = abs(f) < 5*d;                               % exclude resampler images
sigma_psd = sqrt(sum(P(m).*f(m).^2)/sum(P(m)));
% (b) autocorrelation small-lag fit: R(tau) = exp(-2*pi^2*sig^2*tau^2)
sig_ac = 0; L = 0;
for lag_s = [0.05 0.1 0.15 0.2 0.25]
  k = round(lag_s*Fs);
  R = abs(sum(conj(s(1:end-k)).*s(1+k:end)) / sum(abs(s).^2));
  sig_ac += sqrt(-log(R)/(2*pi^2*(k/Fs)^2)); L += 1;
end
sig_ac /= L;
printf("specified 2-sigma: %.1f Hz | realized PSD: %.3f Hz | autocorr: %.3f Hz\n", ...
       d, 2*sigma_psd, 2*sig_ac);

Output on current main (Octave 8.x):

specified 2-sigma: 1.0 Hz | realized PSD: 0.715 Hz | autocorr: 0.709 Hz

i.e. right on the 1/sqrt(2) = 0.707 prediction.

Impact

Every fading simulation built on doppler_spread (the ch tool's
MPP/MPD-style channels, freedv/FreeDATA test campaigns) runs its fading
~29% slower than labeled — e.g. a "1 Hz poor channel" is effectively a
~0.7 Hz channel. Relative comparisons made with the same tool are
unaffected; absolute comparisons to other simulators (PathSim, IONOS,
MIL-STD-compliant instruments) or to published performance curves are
biased optimistic, since slower fading is generally easier on a modem.

Suggested fix

Minimal (one token — keeps the design approach):

b = fir2(Ntaps-1, x/(lowFs/2), sqrt(y));

Or the MIL-STD-188-110C App E closed-form time-domain taps, which also
remove the frequency-sampling design approximation entirely (note this
changes the filter length, so the Ntaps used for the memory-fill
bookkeeping needs to follow):

tau = sqrt(log(100))/(pi*d);                    % 1%-of-peak truncation
t = (-ceil(tau*lowFs):ceil(tau*lowFs))/lowFs;
b = sqrt(2)*exp(-pi^2 * t.^2 * d^2);
b = b/sqrt(sum(b.^2));                          % unit noise-power gain

Either way, a regression test on the realized spread (autocorrelation fit
or fine-bin second-moment PSD, per the repro above) would keep this pinned;
the current amplitude-response unit test cannot catch it.

Happy to send a PR for whichever form you prefer.

[drowe67]

@spinkham - thank you that's a good catch, and will impact several projects such as the RADE/BBFM work. Feel free to submit a PR for the simplest fix - it would be interesting to see if any tests break. Thanks also for the suggested test, that's a fine idea.

I'm busy with other work and personal commitments so it may be some time before this one is fully run to ground across several projects. It has come at a good time for the RADE V2 work, as I will be doing some retraining over the next few months.