2
$\begingroup$

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))
$\endgroup$
4
  • 1
    $\begingroup$ Does this answer your question? dsp.stackexchange.com/a/52960/21048 $\endgroup$ Commented Mar 16, 2023 at 21:17
  • $\begingroup$ Thanks. It helps me to realize correction factor should be the non-coherent correction factor. Updated my post to reflect this. $\endgroup$
    – John L
    Commented Mar 19, 2023 at 17:01
  • $\begingroup$ So is there still a question? If the problem is solved and the answer helped, great - you should +1 and accept the answer. Otherwise edit the question, but note that StackExchange isn't for extended back-forths, and it's preferred to ask a new question for new/developed problems. And FYI, the right side of your equations may be correct for what you're trying to achieve, but none of the equations themselves are correct; it's $\sum |x|^2 = \frac{1}{N} \sum |X|^2$. $\endgroup$ Commented Mar 20, 2023 at 14:37
  • 1
    $\begingroup$ Thank you. I just up voted your answer and Dan Boschen's comment. I don't think your answer itself completely answered the question so I will not accept it as an answer to not confuse future viewer. $\endgroup$
    – John L
    Commented Mar 20, 2023 at 15:51

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ Both my original Parseval equation and your edited Parseval equation are correct. Mine is implied be with Uniform window, yours is explicitly with window w(n). However, when talking about a signal's average power rms, Uniform window is used in time domain. That is why I am trying to correct non Uniform window (Kaiser in this case) to match with Uniform window time domain rms. Please correct my post back to its original form. $\endgroup$
    – John L
    Commented Mar 18, 2023 at 18:59
  • $\begingroup$ Sorry, I just noticed that it was not modified by you. It was the first editor who incorrectly changed my equations. So my original post's equation was correct, and your equation is also correct. $\endgroup$
    – John L
    Commented Mar 19, 2023 at 1:18
  • $\begingroup$ @JohnL The original isn't correct either, $1/S_1$ should be $1/N$. $\endgroup$ Commented Mar 19, 2023 at 10:37
  • $\begingroup$ Thanks. You are correct, it should be $\frac{1}{N}$. Updated my post. $\endgroup$
    – John L
    Commented Mar 19, 2023 at 17:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.