DynamicalPowerspectrum and DynamicalCrossspectrum cannot select the per-segment mean (use_common_mean) #986

Description of the Issue

AveragedPowerspectrum and AveragedCrossspectrum accept use_common_mean, which
decides whether the normalisation uses the mean count rate of the whole observation or
of each segment. DynamicalPowerspectrum and DynamicalCrossspectrum build their
matrix from exactly those classes, but neither accepts the argument nor passes it on, so
every column is normalised by the common mean and the other setting is unreachable:

DynamicalPowerspectrum(lc, segment_size=8.0, norm="frac", use_common_mean=False)
# TypeError: DynamicalPowerspectrum.__init__() got an unexpected keyword argument 'use_common_mean'

In DynamicalPowerspectrum._make_matrix (and the equivalent in
DynamicalCrossspectrum):

avg = AveragedPowerspectrum(
    lc,
    dt=self.sample_time,
    segment_size=self.segment_size,
    norm=self.norm,
    gti=self.gti,
    save_all=True,          # use_common_mean is never passed, so it defaults to True
)

The class docstrings document norm but not which mean the normalisation uses, so there
is nothing to indicate that a choice is being made.

This is seems like a missing option that really should be there.

Justification for change

With norm="frac" and a common mean, each segment's power is divided by the mean rate
of the whole observation, so a changing count rate rescales the dynamical spectrum even
when the source's fractional variability is constant. X-ray binaries and AGN follow a
linear rms-flux relation (Uttley & McHardy 2001; Uttley, McHardy & Vaughan 2005), so
fractional rms is the measurement of variability that stays constant while the rate moves.
Here are three cases where the distinction changes a result:

• Dividing each segment by its own mean removes the rms-flux trend, which is what you
  want when asking whether the fractional variability changes with time.
• Tests of spectral stationarity compare the power in one segment against another. With
  a common mean, a rate change alone shifts every column, which such a test reads as
  non-stationarity. (I was doing work related to this upon discovering this issue).
• The Poisson noise level of a segment, 2/rate in frac normalisation, follows that
  segment's own rate. Subtracting a single level from a dynamical spectrum normalised by
  a common mean is inconsistent, because the noise level and the normalisation refer to
  different means.

The current workaround is to call AveragedPowerspectrum(..., use_common_mean=False, save_all=True) and transpose cs_all, which duplicates what the dynamical class already
does, including the GTI handling and the time axis.

Reproduction

A 512 s light curve with a 2 Hz sinusoid at a fixed fractional rms of 10 %, Poisson
sampled, whose mean count rate doubles at the halfway point. The fractional variability is
constant, so with norm="frac" the columns should not change.

import numpy as np
from stingray import Lightcurve
from stingray.powerspectrum import AveragedPowerspectrum, DynamicalPowerspectrum

rng = np.random.default_rng(0)
dt, T, f0, frac_rms = 1 / 64, 512.0, 2.0, 0.10
n = int(T / dt)
t = (np.arange(n) + 0.5) * dt
mean_rate = np.where(t < T / 2, 200.0, 400.0)
rate = mean_rate * (1 + frac_rms * np.sqrt(2) * np.sin(2 * np.pi * f0 * t))
lc = Lightcurve(t, rng.poisson(rate * dt).astype(float), dt=dt,
                skip_checks=True, gti=np.array([[0.0, T]]))

dyn = DynamicalPowerspectrum(lc, segment_size=8.0, norm="frac")
avg = AveragedPowerspectrum(lc, segment_size=8.0, norm="frac",
                            use_common_mean=False, save_all=True, silent=True)
per_seg = np.array(avg.cs_all).T

j = np.argmin(np.abs(dyn.freq - 2.0))
first = dyn.time < T / 2
for name, mat in (("common mean", dyn.dyn_ps), ("per-segment mean", per_seg)):
    lo, hi = mat[j, first].mean(), mat[j, ~first].mean()
    print(f"{name:18s} 2 Hz: {lo:.4f} -> {hi:.4f}, ratio {hi / lo:.2f}")

Output on stingray 2.3.2 and on current main:

common mean        2 Hz: 0.0416 -> 0.1578, ratio 3.79
per-segment mean   2 Hz: 0.0928 -> 0.0890, ratio 0.96

The count rate ratio between the halves is 1.99 and the expected power ratio at constant
fractional rms is 1.00. The dynamical spectrum changes by the square of the rate ratio;
the AveragedPowerspectrum workaround gives 0.96, consistent with 1.

Your suggestions for the Fix

Add use_common_mean=True to both __init__ methods, pass it into the averaged
spectrum, and say in the docstrings which mean is used. The default keeps the current
behaviour.

One thing worth noting: _make_matrix also stores

conv = avg.cs_all / avg.unnorm_cs_all
self.unnorm_conversion = np.nanmean(conv)

which compute_rms uses to recover the unnormalised periodogram. With a per-segment
mean this factor differs from segment to segment (by a factor of 4.7 across the light
curve above), so a single number no longer describes it well and compute_rms would return
the wrong rms. It has to become a per-segment array in that case.

I have a branch with the change and tests and will open a PR.