How does the effect of windowing change with the phase of the input signal?

Question

I am calculating the SNR (signal power to noise power) for a sine wave. I don't have an integer number of periods in the waveform being analyzed, so I am using a flattop window to reduce spectral leakage.

If I start with a long sine wave and select a subset of a given length, I find that the calculated SNR vs the starting index of the subset (i.e. the phase of the segment being analyzed) is a periodic function with a period equal to 1/2 period of the sine being analyzed.

You can see that the best SNR corresponds to a segment that starts and ends at zero crossings, and the worst case corresponds to a segment that starts and ends at a peak. The spectra of the two cases primarily differs by the power of the first few DFT points.

Why is this happening? I understood the effect of windowing in the frequency domain to be the convolution of the spectrum of the window with the spectrum of the signal—additional tones should not be generated.

Code and figures attached below.

import numpy as np
from scipy import signal
from matplotlib import pyplot as plt

# Number of samples to analyze.
N = 2**14
# Not quite an integer number of cycles in N samples. 
f = 12.98/N
# Complete data set is longer by 5 periods.
t = np.arange(0,N + int(5/f))
x = np.cos(2*np.pi*f*t) + 0.0001*np.random.randn(len(t))
w = signal.get_window('flattop',N)
# Main lobe width of flattop window
mlw = 9
half_mlw = (mlw-1)//2
starting_idx = np.arange(0, int(2/f), int(1/f / 100))
snr_list = np.zeros(len(starting_idx))
for i, start_idx in enumerate(starting_idx):
    subset = x[start_idx : start_idx + N]
    freq, ps = signal.periodogram(subset, window=w, detrend='constant')
    signal_idx = np.argmax(ps)
    signal_points = list(range(signal_idx - half_mlw, signal_idx + half_mlw + 1))
    signal_power = sum(ps[signal_points])
    noise_spectrum = ps.copy()
    noise_spectrum[signal_points] = 0
    noise_power = sum(noise_spectrum[1:-1])
    snr = 10*np.log10(signal_power/noise_power)
    snr_list[i] = snr

# Plot SNR vs phase of input waveform.
fig1, ax1 = plt.subplots()
ax1.plot(360 * starting_idx / (1/f), snr_list, '.-')
ax1.set_xlabel("Phase (deg)")
ax1.set_ylabel("SNR (dB)")
_, xmax = ax1.get_xlim()
ax1.set_xticks(np.arange(0, xmax, 90))
ax1.set_title("SNR vs Phase")

best_start = np.argmax(snr_list)
worst_start = np.argmin(snr_list)    

for start_idx, title in [(best_start, "Best SNR"), (worst_start, "Worst SNR")]:
    start = starting_idx[start_idx]
    subset = x[start : start + N]
    freq, ps = signal.periodogram(subset, window=w, detrend='constant')
    signal_idx = np.argmax(ps)
    signal_points = list(range(signal_idx - half_mlw, signal_idx + half_mlw + 1))
    fig,ax = plt.subplots(2,1)
    plt.suptitle(title)
    # Show points ID'd as signal points.
    # ax[0].plot(freq[signal_points], 10*np.log10(ps[signal_points]), 'mo')
    ax[0].plot(freq, 10*np.log10(ps),'.-')
    ax[0].set_xscale('linear')
    ax[0].set_xlim(0,3*freq[signal_idx])
    ax[1].plot(np.arange(N), subset)

Answer 1

Summary

The problem is caused by detrend='constant' in the call to periodogram, which subtracts the mean of the input before windowing, which actually injects a mean that's controlled by the spectral leakage of a rectangular window. The "vs Phase" pattern is sinusoidal due to the time-shift property of the Fourier transform: shift in time $\Leftrightarrow$ sinusoidal modulation in Fourier, and a phase shift is equivalently a time shift for a pure sine. However, it's also about tonal interference: no sines for a complex sinusoid.

detrend=False fixes the problem. For a sine, though, a better option is available - retrieving all sine parameters (frequency, phase, amplitude), using inverse problem solved from the sine DFT solution - the Two Bin Solution, which will be very accurate for such high SNR. It can also improve SNR estimation (just subtract the sine from the noisy signal). For a competitive and extensive comparison (frequency only), see Sine Frequency Estimation under White Noise.

Related considerations / other approaches (marginal improvements):

1. Increase the "signal" interval (mlw -> signal_points) and compute it automatically based on some measure of "effective support".
2. To minimize leakage, use a window with high leakage rejection (high frequency resolution + rapid side-lobe dropoff), but also that most accurately reflects the continuous function it's sampling, i.e. ensure it's fully decayed. This way, if a window guarantees some leakage reduction in continuous time, it's actually realized.
3. Use higher frequency sine, and ensure the duration of the sliding FFT has an integer ratio relationship with the sine's frequency. This does have a considerable effect on the extreme (200+) dB scales, but if other stuff's done right, then minimal. (In your code, 12.98 -> 13 for f. Worst SNR goes from 55 to 293 dB, even with detrend='constant'!)
4. For a pure sine, there's no need to look at a spectrogram (what you're doing), unless using it for e.g. Welch's method - windowing the entire signal will provide the widest possible window and the most leakage reduction. The effect of a longer window upon leakages at extreme dB scales is very marginal, however - but that's not discussing noise.

Updated explanation + proof

A sine DFT closed form solution enabled proving all claims in this answer. Below are exactly true:

(1): Real and imaginary parts are sinusoidally modulated with frequency $f$, as a sine shifts. All imaginary bins are modulated same, reals are bin-dependent.

(2): $|X|^2$ is sinusoidally modulated with frequency $2f$, as a sine shifts.

(3): Bin energy ratios,

$$ \texttt{ER}_I\{X_\tau\} = \frac{ \sum_{k\in I}|X_\tau[k]|^2 } { \sum_{k\notin I}|X_\tau[k]|^2 }, $$

i.e. how the question measures SNR (and how "leakage" can be measured), are modulated with period $2\pi/(2f)$, sometimes near-sinusoidally:

$$ \boxed{ \begin{align} & \texttt{ER}_I\{X_\tau\} = \frac{C_0\cos(2\pi2\tau/N + q_0) + D_0} {C_1\cos(2\pi2\tau/N + q_1) + D_1} \\ & C_0, C_1, D_0, D_1, q_0, q_1 = \texttt{f}\{f, N, \phi, I\} \\ & |C_0| \leq D_0,\ |C_1| \leq D_1 \end{align} } $$

where $C, D$'s are sine-determined but same for the same sine, and $\tau$ is shift in samples. $\texttt{ER}$ is the least sinusoidal when the denominator is most sinusoidal, i.e. $|C_1| \approx D_1$. This happens when the numerator's interval, $I$, includes a larger interval near peak of the sine - turns out, the bins farthest away from peak are the most sinusoidally-modulated (have lowest DC offset). Lowering mlw from 9 to 8 in OP's code produces a sinusoidal pattern.

All proofs, and further explorations, are found in exploring the sine DFT modulus - see under "Proof: shifting modulates energy ratios ...".

What of windowing? Subtracting DC and then windowing, simply inserts a copy of the window at DC, scaled by the subtracted value. This scaling (DC) is sinusoidally modulated, hence the sum over the window is also sinusoidally modulated, exactly same as without windowing.

Turning off detrend forces a sinusoidal pattern with same mlw, since the windowing, without DC injection, drives the denominator's sinusoidal modulation's DC offset (i.e. $D_1$ relative to $|C_1|$) through the roof, while the numerator retains sinusoidal modulation (and it doesn't matter what its $D_0$ relative to $|C_0|$ is, as the linked proof explains, but it's not nearly as great).

The rest of the answer still has insights that this section doesn't cover (and the closed form solution can't explain), but OP's exact question is fully answered with this section alone. The rest of the answer also provides a perspective from continuous time (sampling, FT, DTFT).

---

Problem explanation

First, for the general reader - the "SNR" plot is generated by dividing a selected segment of ps (signal) by the non-selected segment (noise), where ps = abs(fft(subset))**2 with detrend=False or abs(fft(subset - mean(subset)))**2 with detrend='constant', and subset = x[start_idx:start_idx + N] where start_idx iterates starting_idx and of course len(x) > N (sub-ideal variable naming).

This image from Spectral leakage Wiki nicely summarizes most of it:

(an important distinction is that the image's blue vs red shows different durations, hence effective frequencies, as opposed to different shifts/phases, so this doesn't reflect "vs shift")

Since DC is subtracted before windowing, it's same as subtracting DC after a rectangular window. A rectangular window generates sincs centered at each of the pure sine's impulses in frequency domain. Those sincs overlap and add over all of frequency, and they have an awful, $\propto 1/f$ decay, making most workarounds have marginal effects at extreme dB scales.

Applying the actual window convolves the window's spectrum with the DC-injected input's, which gets counted as "noise" in the SNR computation.

DC, however, isn't special - subtracting bins next to it will have a similar effect, better or worse depending on the exact configuration. In your code, adding means.append(subset.mean()) and bin1s.append(fft(subset)[1]/len(subset)) tracks the mean and (normalized) bin 1, and yields

The severity of the "injection" depends on the configuration (sine frequency, window duration, hop locations, etc), and in best case (very unlikely) can be a non-issue.

As a summary, I like this sentence from @JasonC's answer:

In effect, the detrend gives us the spectral leakage characteristics at DC of a rectangular window even when we use a different window.

Why the sinusoidal pattern?

Sliding the window is same as keeping the window in place and sliding the sine. Sliding the sine is time-shifting the sine, and $x(t - c) \Leftrightarrow e^{-j\omega c} X(\omega)$. So, the sincs are sinusoidally modulated, with shift-dependence. That is, the shifts change the frequency of said modulation, and if the shifts are uniform (same hop size), it traces out a complex sinusoid of a fixed frequency for any given bin (except DC, whose imaginary part must be zero for a real-valued input).

This can be proven and made more explicit. I derived an analytic solution for DFT of a sine, which can be used to predict the behavior of the DC (or any other) bin for any sliding FFT configuration; plugging in for the DC bin ($k = 0$):

$$ \texttt{DFT}\{\cos(2\pi f_0 (t + \tau))\}[0] = \frac{T}{2s} \sum_{l=-\infty}^{\infty} \bigg( \\ \begin{align} & \left[ \exp\left\{ j2\pi \left( \frac{T - s}{2} \left(\frac{l}{s} + f_0\right) + \tau f_0 \right) \right\} \text{sinc}\left(-\pi T \left(\frac{l}{s} + f_0\right) \right) + \right. \\ & \left. \ \ \exp\left\{ j2\pi \left( \frac{T - s}{2} \left(\frac{l}{s} - f_0\right) - \tau f_0 \right) \right\} \text{sinc}\left(-\pi T \left(\frac{l}{s} - f_0\right) \right) \right] \bigg) \end{align} $$

where $T = N \cdot s$ is the duration of subset, s = t[1] - t[0], N = len(subset), and $\tau$ corresponds to your start_idx. Indeed I originally formulated the derivation around your exact setup, and your exact setup is still exactly described by the result.

The expression shows dependence on size of subset via $T$, sine's frequency via $f_0$, and choice of starting_idx via $\tau$. And of course the sampling rate, $1/s$. Simplifying a little and setting $l=0$,

$$ \begin{align} & \left[ \exp\left\{+j2\pi f_0 (T + \tau)\right\} \text{sinc}(-\pi T f_0) + \right. \\ & \left. \ \exp\left\{-j2\pi f_0 (T + \tau)\right\} \text{sinc}(+\pi T f_0) \right] \end{align} $$

The only part that depends on the shift, $\tau$, is the sinusoidal modulation. Indeed, as a function of $\tau$, we expect modulation with frequency $f_0$, phase shift $T$, and constant scaling by $\text{sinc}$ that depends on $T$ and $f_0$. The other $l$ are simply shifted copies of $l=0$. It's true that it's no longer clear whether the result remains an exact sine, but it's certainly close enough in the examples inspected (I inspected many other $f_0$, $T$, and $\tau$).

A bold prediction is that mean(subset) should be of same frequency as subset... and a correct one:

(I equalized amplitudes for plotting)

This covers the gist, but is actually incomplete; there's no sinusoidal pattern for SNR of a complex sinusoid! In fact we've only shown that the mean is sinusoidal, not SNR. See "Completed explanation (part 2)" below.

Maximum SNR?

• detrend=False
• Remove noise
• Gaussian window (not saying it's best for your purpose)
• large and integer f
• larger mlw (number of points counted as "signal")

(didn't bother labeling phase) That's $6 \times 10^{28}$ max SNR. We know we're abusing floats because it's no longer sinusoidal; if we're a little less extreme, the sine persists. Now with OP's original with noise:

which matches the ground truth 10*log10(x.var() / (0.0001*randn).var()) = 10*log10(0.5 / 0.0001**2) of 77.0.

If must detrend, before or after windowing?

After, but ideally not at all.

The need to detrend really indicates the need for filtering others frequencies also, as it's unlikely to be just DC. If DC and frequencies around it are a problem, a viable option is a CWT-based power spectrum; wavelets (if done right) are mean-rejecting by definition and have phenomenal decay leading up to the mean. It's not only unlike STFT, but STFT is straight up incapable of doing this. Make sure it's L1 normed and avoid PyWavelets & scipy.

Completed explanation (part 1)

Above's not actually done yet. So what if we "inject DC" in continuous time? Per DFT, we're literally removing it. In short, the effect in this context is the same.

To see what's up, we amend the DFT derivation by replacing the result for $\hat x(f)$ with

$$ \hat x(f) = \mathcal{F}\{w(t - (T - s)/2) (\cos(2\pi f_0 (t + \tau)) - \mu)\} $$

where $w(t)$ is zero outside of $I=[-s/2, T-s/2]$ (the DFT interval in continuous time), and $\mu$ is the mean of the same cosine over $I$. We're already familiar with the mean, so let's just say it's some small constant that depends on $\tau$ (the only thing we care about). This modifies the integral in frequency as: (1) replace sinc with $\hat w(f)$, (2) add $-\mu \hat w(f)$ to end-result (in the integral, $-\mu \delta(f)$ is added to the right side).

We already knew (1), the key is (2). This $\mu \hat w(f)$ won't cancel in the infinite sum and will pop up in the DFT, so we expect a mean-scaled copy of the window's frequency response planted at DC, which is what your plots show.

Still, purely in terms of DFT, since multiplying by a window is convolving with its spectrum in frequency, shouldn't removing a bin also remove its copy at the bin? Yes, but only if all other bins are zero within the support of fft(window), yet we established that DC isn't special and is very similar to adjacent bins. This is explained in part 3.

Completed explanation (part 2)

Next, we've just shown the mean is sinusoidal - what about the ratio ("SNR")? One can show, the signal bin is also sinusoidal with same frequency, and $\sin(\omega t) / (\sin(\omega t) + C)$ is approximately proportional to $\sin(\omega t)$ for large $C$ (i.e. ratio is sinusoidal with same frequency). The $+C$ of course comes from the fact that we're not dividing individual complex bins, but sums of square moduli of bins. That's a lot messier to show exactly mathematically, but not on high level:

If all bins are scaled by the same factor for all $\tau$, and this scaling is sinusoidal as a function of $\tau$, then so is the sum, and it has the same frequency as the DC or any other bin (by "same scaling" premise). To check the sameness premise, we check dependence on $k$ (bin index); that makes our previous expression (one term) $\exp\left\{+j2\pi f_0 (T + \tau - k)\right\} \text{sinc}(-\pi T (f_0 - k))$. We need the ratio for different $\tau$ to be the same for all $k$ - and it is: the sinc cancels, and the $k$ (along $T$) in the exp cancels. Yet... the modulus for all bins is independent of $\tau$ (the exp does nothing)!

This proves that the sinusoidal SNR is due to tonal interference; indeed, replacing np.cos( in code with np.exp(-1j* turns the SNR into a flat line:

What about the mean that we've shown earlier? Yes, it too. Frankly when I first answered I forgot that it's about modulus. The detrend='constant' will now simply lower SNR.

But wait, what of modulus of the interference? Since "same scaling" holds of each tone, their complex sum - and hence their modulus, and thereby sum of modulus - will retain sinusoidal dependence on $\tau$.

Completed explanation (part 3)

Left is spectrum of windowed segment, right is of the segment and the window:

One frame uses detrend='constant', other detrend=None. win * segment convolves fft(win) with fft(x) (x here is segment), yet inspecting the GIF, the result doesn't add up: the DC bin is nearly equal to adjacent bins. This is where phase kicks in; now real and imag parts:

Convolution of a and b is just sum of a * shift(b, i) for all i. So the DC bin of fft(win * x) is fft(x) times fft(win), where fft(win) is centered over the DC bin of fft(x):

fft(x)[0] is the plot summed. Note the primary signal here is obliterated by the window's excellent frequency resolution, and what remains is leakage. The image shows the 'constant' case, with DC zeroed; its removal matters because fft(win) effectively behaves like [1, -1, 1, -1], so every non-zero input makes a huge difference. Yet, removal of any bin adjacent to DC would have a nearly identical effect.

This explains the behavior in terms of DFT and resolves the apparent "paradox". I don't think it's possible, however, to do this for more than just one case at a time and prove dependencies on $\tau$ as we have. Sometimes the DFT has a DC, other times not, and adjacent bins also vary. It's what motivated my analytic solution.

Answer 2

What's going on?

I think the main issue is the use of detrend='constant' in the calls to signal.periodogram. If I change that to detrend=False (at both call sites), I get a much smaller variation in the SNR vs phase plot:

Why does detrend='constant' cause a problem?

Well, if you dig through the various helper functions called by signal.periodogram, you will eventually find yourself in a function named _fft_helper in a file named "_spectral_py.py" which contains this code:

    # Detrend each data segment individually
    result = detrend_func(result)

    # Apply window by multiplication
    result = win * result

    # Perform the fft. Acts on last axis by default. Zero-pads automatically
    if sides == 'twosided':
        func = sp_fft.fft
    else:
        result = result.real
        func = sp_fft.rfft
    result = func(result, n=nfft)

The detrend is applied before windowing! So the unwindowed data has its DC bias removed, but after windowing, all bets are off! The detrend can actually (depending on phase) create a much larger DC term than would exist without the detrend, accomplishing the opposite of the intended effect.

The energy injected at DC is being interpreted as noise by your SNR calculations, which is why your SNR is varying so much with phase.

Put another way, when the mean of the signal is calculated for the detrend operation, it is as if the signal has already been windowed with a rectangular window. The spectral leakage is quite bad, and the energy at all frequencies varies due to the phase-dependent constructive/destructive interference of the two (remember, there is negative frequency too!) input tones.

Alternatively, we can consider the time domain: The input tone will have some integer N number of "whole cycles" which will have a mean of zero, and then there will be a fraction of a cycle at the end, which will have a phase-dependent mean.

Normally, this energy would get mostly cancelled when we apply our actual window, but if we detrend first, the magnitude of the frequency domain convolution at DC changes by the amount of energy we subtract. In effect, the detrend gives us the spectral leakage characteristics at DC of a rectangular window even when we use a different window.

This behavior of SciPy seems pretty bad to me, but I am very ignorant about statistical spectrum analysis, Welch's method, etc., so perhaps there is a good reason for this.

Without detrending, why is there still some variation in SNR vs phase?

The SNR vs phase plot takes the form of a sinusoidal shape superimposed on a random-walk-like pattern, with the latter varying from run to run.

The sinusoidal variation is due to the FFT being complex-valued, and the negative frequency component of your real-valued test tone is contributing energy to the positive frequency FFT bins in a phase-dependent way.

The random walk is due to the signal having noise. When the noise is integrated by overlapped FFTs, the overlap leads to a correlation in the noise, creating the random walk effect.

To reduce the impact of these effects, you could increase the FFT length, switch to an FFT window with better sidelobe rolloff (like a Kaiser window), or both.