"""Exact variance of contiguous Hann-window periodogram-power averages.

Python 3.14, NumPy and Matplotlib; see pyproject.toml. Writes the PNG to cwd.
The model uses proper complex Gaussian white sample noise and the original
N=64 frequency grid. Blocks contain at most 32 bins, avoiding cyclic overlap.
"""
from pathlib import Path
import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt


def normalized_mean_variance(size):
    size = np.asarray(size, dtype=float)
    return (size + 2 * np.maximum(size - 1, 0) * (4 / 9)
            + 2 * np.maximum(size - 2, 0) * (1 / 36)) / size**2


def main():
    count = np.arange(1, 33)
    variance = normalized_mean_variance(count)
    fig, axes = plt.subplots(1, 2, figsize=(9, 4.2), dpi=120, layout='constrained')
    fig.patch.set_facecolor('white')
    axes[0].plot(count, count, '--', color='#555555', label='Independent powers')
    axes[0].plot(count, 1 / variance, color='#176da5', linewidth=2.3,
                 label='Periodic Hann powers')
    axes[0].scatter([10], [1 / normalized_mean_variance(10)], color='#176da5', s=25)
    axes[0].annotate('10 bins: about 5.42 independent powers',
                     xy=(10, 1 / normalized_mean_variance(10)), xytext=(3, 20),
                     arrowprops={'arrowstyle': '->', 'color': '#333333'}, fontsize=8.5)
    axes[0].set(xlabel='Number M of neighboring bins averaged',
                ylabel='Effective independent-power count', xlim=(1, 32), ylim=(0, 33))
    axes[0].legend(loc='upper left', fontsize=8)
    axes[1].semilogy(count, 1 / count, '--', color='#555555', label='Independent powers')
    axes[1].semilogy(count, variance, color='#176da5', linewidth=2.3,
                     label='Periodic Hann powers')
    axes[1].set(xlabel='Number M of neighboring bins averaged',
                ylabel='Variance of mean power / single-power variance', xlim=(1, 32))
    axes[1].legend(fontsize=8)
    for axis in axes:
        axis.grid(alpha=0.2)
    fig.suptitle('Window-induced correlation reduces spectral averaging gain', fontsize=12.5)
    fig.savefig(Path('hann-periodogram-power-averaging.png'), facecolor='white', transparent=False)
    plt.close(fig)


if __name__ == '__main__':
    main()