I am using an implementation of the averaged cyclic periodogram (section 3.2.4 in Antoni, Jérôme. "Cyclic spectral analysis in practice." Mechanical Systems and Signal Processing 21.2 (2007): 597-630.). The cyclic spectrum is estimated based on two $L$ ln $$\hat{P}_{YX}^{(W)}\left(f,\alpha;L\right)=\frac{1}{K\Delta}\sum_{k=0}^{K-1}X_{N_w}^{k}\left(f+0.5\alpha \right)X_{N_w}^{k}\left(f-0.5\alpha \right)^*$$ where $$X_{N_w}^{k}\left(f\right)=\Delta\sum_{n=kR}^{n=kR+N_w-1}w_k[n]y[n]e^{-j2\pi fn\Delta}$$
This implementation is based on the attached scipy pull request.
Now, if I apply this implementation with the Hann window, on normal noise, I receive an artifact elevation at the modulation frequency representing the half of the window length ($\frac{2F_s}{N_w}$). Seems like a Boxcar window eliminates this problem. Is it the only and the best solution?
import numpy as np
from scipy.signal import csd, windows
import matplotlib.pyplot as plt
fs = 20000
s = np.random.normal(size=fs)
alpha = np.arange(5, 200)
# window = windows.boxcar(256)
window = windows.hann(256)
BiSpectrum = []
for alpha1 in alpha:
x = s * np.exp(-1j * np.pi * (alpha1 / fs) * np.arange(s.shape[-1]))
y = s * np.exp(1j * np.pi * (alpha1 / fs) * np.arange(s.shape[-1]))
f, Pxy = csd(x, y, fs=fs, window=window)
Pxy = Pxy[f > 0]
f = f[f > 0]
BiSpectrum.append(Pxy)
BiSpectrum = np.abs(BiSpectrum).T
plt.figure(figsize=[19, 9])
plt.pcolormesh(alpha, f, BiSpectrum, vmax=np.percentile(BiSpectrum, 99.99))
plt.xlabel('Modulation frequency $\\alpha$ [Hz]')
plt.ylabel('Carrier Frequency $f$ [Hz]')
plt.colorbar()
plt.show()
