PSD window correction factor?

[lukexor]

In the section on IQ sampling there is the following code for calculating the PSD and it follows with a note about optionally applying a window.

Shouldn't there be a window correction factor applied to account for the reduction in amplitude resulting from the window?

I've seen a number of discussions on window correction factors and I'm unsure which applies in this case (coherent gain or equivalent noise bandwidth). E.g. dividing by the sum of the square of the window, the square of the sum, of the mean, etc.

I couldn't find any reference to it in this section, unless I missed it.

Fs = 1e6 # lets say we sampled at 1 MHz
# assume x contains your array of IQ samples
N = 1024
x = x[0:N] # we will only take the FFT of the first 1024 samples, see text below
x = x * np.hamming(len(x)) # apply a Hamming window
PSD = np.abs(np.fft.fft(x))**2 / (N*Fs)
PSD_log = 10.0*np.log10(PSD)
PSD_shifted = np.fft.fftshift(PSD_log)

[777arc]

Hmm that's a good point, I think the correction factor for preserving signal power would be the one to do, so I believe it's dividing by the sum of the squared windowing function. Although that would mean for no/rect window you would have to divide by N, which we are already doing later, so perhaps the scaling factor is 1/(np.sum(window**2) / N)

If you're able to figure it out with some example code to prove it works, it would be a great contribution to the open source textbook 😄

[lukexor]

Yeah I'll post back if I come to something. I've been using a variant of this correction in a Rust/wasm spectrum analyzer I'm building and while it looks fine, I'm still fairly new to dsp and not sure how to validate it's correct but the amplitudes seem similar to what you get in analyzers like SDR±±

[lukexor]

Not sure if this is accurate - but it feels correct:

ENBW = np.sum(window**2) / np.sum(window)**2  # Equivalent Noise Bandwidth (normalized)
PSD = np.abs(np.fft.fft(x*window))**2 / (N*Fs*ENBW*N)

For the example at the end of the IQ Sampling section:

idx_50Hz = np.argmin(np.abs(f - 50))

print(f"PSD at 50 Hz peak (hanning, ENBW corrected): {PSD_shifted[idx_50Hz]:.2f} dB")
print(f"ENBW correction: {10*np.log10(ENBW*N):.2f} dB")

It prints:

PSD at 50 Hz peak (no window): 6.69 dB
PSD at 50 Hz peak (hanning, ENBW corrected): -0.07 dB
ENBW correction: 1.76 dB

I also compared to https://www.gaussianwaves.com/2020/09/equivalent-noise-bandwidth-enbw-of-window-functions/ and checked their results for a Hamming window and get similar numbers.

[777arc]

Interesting, ends up being normalized by N^2 because of the extra N from the ENBW

[lukexor]

So I think I was looking for the coherent power gain: https://www.gaussianwaves.com/2020/09/window-function-figure-of-merits/

I suspect you need to include both for full accuracy. Gonna test some more and compare

[777arc]

I would imagine a quick way to test is take the same signal and put it through different windows (incl rect) and see if the levels remain constant

[lukexor]

So far I've just been eyeballing signals from my SDR - but if I find the time later I'll try and make some known test cases.

This is currently what I'm doing in my Rust/WASM spectrum analyzer:

    fn coherent_power_gain_db(window: &[f32]) -> f32 {
        let cpg = window.iter().sum::<f32>() / window.len() as f32;
        20.0 * cpg.log10()
    }

    fn equivalent_noise_bandwidth_db(window: &[f32]) -> f32 {
        let len = window.len() as f32;
        let window_sq_sum = window.iter().map(|v| v.powi(2)).sum::<f32>();
        let window_sum_sq = window.iter().sum::<f32>().powi(2);
        let enbw = len * window_sq_sum / window_sum_sq;
        10.0 * enbw.log10()
    }

Then in my FFT processing function

        self.fft
            .process_with_scratch(&mut self.fft_buf, &mut self.fft_scratch);

        fft_out
            .iter_mut()
            .zip(&self.fft_buf)
            .for_each(|(psd_log, &sample)| {
                // https://kluedo.ub.rptu.de/frontdoor/deliver/index/docId/4293/file/exact_fft_measurements.pdf
                *psd_log =
                    10.0 * (sample * self.norm_factor).norm_sqr().log10() - self.cpg - self.enbw;
            });

norm_factor is just 1/fft_size