Correct Windowing Effect at Amplitude Scale

Question

I am trying verify the noise floor returned by Matlab sinad(). I am able to get the results matching by summing power spectrum density with $\frac{f_\text{s}}{N}$. But I am not able to match the results by following Parseval's theorem to root sum square the $X[k]$.

$$ S_1=\sum_{n=0}^{N-1} \big|w[n]\big| $$ $$ S_2=\sum_{n=0}^{N-1} \big|w[n]^2\big| $$

where $w[n]$ is a window function.

$$ \sum_{n=0}^{N-1} \big|x[n]\big|^2 = \frac{1}{S_1}\sum_{k=0}^{N-1}\big|X[k]\big|^2 $$

My thought was to apply windows DC gain correction at $X[k]$ amplitude scale (not power scale). This relation works if the windowing function is Uniform, so obviously I am missing something.

Thanks to @Dan Boschen and @OverLordGoldDragon, I realized there are two errors:

1. Because this is non-coherent signals, the correction factor is $S_2$ rather than $S_1$. And correction is applied to power not amplitude.
2. The dividing factor should still be $frac{1}{N}$. Correction made within the summing operation.

Therefore the Parseval's equation becomes: $$ \sum_{n=0}^{N-1} \big|x[n]\big|^2 = \frac{1}{N}\sum_{k=0}^{N-1}\frac{X[k]^2}{S_2} $$

Here is revised test code. All results match.

rng(1)
fs = 500;   % sampling frequency
T = 1000;
t = 0:1/fs:T-1/fs;
Fs = 100;    % natural frequency
data = cos(2*pi*Fs*t) + 0.011*randn(size(t));
N = length(data);
F=fs*(0:(N/2))/N;

% Kaiser window used by sinad(), beta=38
win = kaiser(N,38);
data_fft_full_abs_k=abs(fft(win.*data'));
data_fft_full_abs_k = data_fft_full_abs_k(1:N/2+1);
S1=sum(win);
S2=sum(win.^2);
ENBW_hz=fs*S2/(S1^2);

% find the range of signal spectrum
[peak_fft,peak_freq_idx]=max(data_fft_full_abs_k);
idxLeft = peak_freq_idx-1;
idxRight = peak_freq_idx+1;
while idxLeft > 0 && data_fft_full_abs_k(idxLeft) <= data_fft_full_abs_k(idxLeft+1)
    idxLeft = idxLeft - 1;
end
while idxRight < N && data_fft_full_abs_k(idxRight-1) >= data_fft_full_abs_k(idxRight)
    idxRight = idxRight + 1;
end
idxLeft = idxLeft + 1;
idxRight = idxRight - 1;
% remove the signals
data_fft_full_abs_k(idxLeft:idxRight)=0;

% calculate rms with amplitude correction
% revert to double sided before calculating root sum square
data_fft_full_abs_k=vertcat(data_fft_full_abs_k,data_fft_full_abs_k(2:length(data_fft_full_abs_k)-1));
noise_fft_rss=rssq(sqrt(data_fft_full_abs_k.^2/S2))/sqrt(N)

% calculate rms with correct at psd
data_fft_full_power=data_fft_full_abs_k.^2;
data_fft_power=(2/S1^2)*data_fft_full_power(1:N/2+1);
data_fft_psd=data_fft_power./ENBW_hz; % 2/(fs*S2)*data_fft_power
noise_psd=sqrt(sum(data_fft_psd*fs/N))

% calculate rms with sinad()
[snrad_dB,dist_noise_pw_dB]=sinad(data,fs);
noise_floor_rms_sinad=sqrt(10^(dist_noise_pw_dB/10))

Answer 1

Some observations:

• Your Parseval equation is incorrect, should be $\frac{1}{N^3}\sum_{k=0}^{N-1}|X[k] * W[k]|^2$, where $*$ is circular convolution and FFT is such that iFFT does 1/N; $w[n]$ can in no way move outside the summation on right hand side
• Below I ignore sqrt everywhere
• noise_fft_rms is doing sum((data_fft_full_abs_k/S1)^2); this does not implement your equation, move S1 out
• noise_psd is doing sum(data_fft_full_abs_k.^2) / (N*S2), which is very different and won't match any formulation in terms of S1 in the general case

Parseval-Plancherel's is just a way of computing signal energy in frequency domain; I don't know what sinad is but it clearly does more than that. The norm by 1/S1 is for matching peak value with signal amplitude, as shown here - that's different from Parseval's. I'd check your definitions.