www.gusucode.com > signal 案例源码程序 matlab代码 > signal/PerformanceOfThePeriodogramExample.m
%% Performance of the Periodogram
% The following sections discuss the performance of the periodogram with
% regard to the issues of leakage, resolution, bias, and variance.
%
% *Spectral Leakage*
%
% Consider the PSD of a finite-length (length $L$) signal $x_L(n)$. It is
% frequently useful to interpret $x_L(n)$ as the result of multiplying an
% infinite signal, $x(n)$, by a finite-length rectangular window, $w_R(n)$:
%
% $$x_L(n)=x(n)w_R(n).$$
%
% Because multiplication in the time domain corresponds to convolution in
% the frequency domain, the expected value of the periodogram in the
% frequency domain is
%
% $$E\{ {\hat P_{xx}}(f)\} = {1 \over {{F_s}}}\int_{ - {F_s}/2}^{{F_s}/2}
% {{{{{\sin }^2}(L\pi (f - f')/{F_s})} \over {L{{\sin }^2}(\pi (f -
% f')/{F_s})}}} \,{P_{xx}}(f')\,df',$$
%
% showing that the expected value of the periodogram is the convolution of
% the true PSD with the square of the Dirichlet kernel.
%
% The effect of the convolution is best understood for sinusoidal data.
% Suppose that $x(n)$ is composed of a sum of $M$ complex sinusoids:
%
% $$x(n) = \sum\limits_{k = 1}^N {{A_k}} {e^{j{\omega _k}n}}.$$
%
% Its spectrum is
%
% $$X(\omega) = \sum\limits_{k = 1}^N {A_k}\delta (\omega-{\omega_k}),$$
%
% which for a finite-length sequence becomes
%
% $$X(\omega ) = \int_{ - \pi }^\pi {\left( {\sum\limits_{k = 1}^N
% {{A_k}\delta (\varepsilon - {\omega _k})} } \right){W_R}(\omega -
% \varepsilon )\,d\varepsilon } .$$
%
% The preceding equation is equal to
%
% $$X(\omega ) = \sum\limits_{k = 1}^N {{A_k}} {W_R}(\omega - {\omega
% _k}).$$
%
% So in the spectrum of the finite-length signal, the Dirac deltas have
% been replaced by terms of the form $W_R(\omega-\omega_k)$, which
% corresponds to the frequency response of a rectangular window centered on
% the frequency $\omega_k$.
%
% The frequency response of a rectangular window has the shape of a
% periodic sinc:
%%
L = 32;
[h,w] = freqz(rectwin(L)/L,1);
y = diric(w,L);
plot(w/pi,20*log10(abs(h)))
hold on
plot(w/pi,20*log10(abs(y)),'--')
hold off
ylim([-40,0])
legend('Frequency Response','Periodic Sinc')
xlabel('\omega / \pi')
%%
% The plot displays a mainlobe and several sidelobes, the largest of which
% is approximately 13.5 dB below the mainlobe peak. These lobes account for
% the effect known as spectral leakage. While the infinite-length signal
% has its power concentrated exactly at the discrete frequency points
% $f_k$, the windowed (or truncated) signal has a continuum of power
% "leaked" around the discrete frequency points $f_k$.
%
% Because the frequency response of a short rectangular window is a much
% poorer approximation to the Dirac delta function than that of a longer
% window, spectral leakage is especially evident when data records are
% short. Consider the following sequence of 100 samples:
fs = 1000; % Sampling frequency
t = (0:fs/10)/fs; % One-tenth second worth of samples
A = [1 2]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)
%%
% It is important to note that the effect of spectral leakage is contingent
% solely on the length of the data record. It is not a consequence of the
% fact that the periodogram is computed at a finite number of frequency
% samples.
%%
% *Resolution*
%
% _Resolution_ refers to the ability to discriminate spectral features, and
% is a key concept on the analysis of spectral estimator performance.
%
% In order to resolve two sinusoids that are relatively close together in
% frequency, it is necessary for the difference between the two frequencies
% to be greater than the width of the mainlobe of the leaked spectra for
% either one of these sinusoids. The mainlobe width is defined to be the
% width of the mainlobe at the point where the power is half the peak
% mainlobe power (i.e., 3 dB width). This width is approximately equal to
% $f_s/L$.
%
% In other words, for two sinusoids of frequencies $f_1$ and $f_2$, the
% resolvability condition requires that
%
% $$f_2-f_1>{{F_s}\over L}.$$
%
% In the example above, where two sinusoids are separated by only 10 Hz,
% the data record must be greater than 100 samples to allow resolution of
% two distinct sinusoids by a periodogram.
%
% Consider a case where this criterion is not met, as for the sequence of
% 67 samples below:
fs = 1000; % Sampling frequency
t = (0:fs/15)/fs; % 67 samples
A = [1 2]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 0.1*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)
%%
% The above discussion about resolution did not consider the effects of
% noise since the signal-to-noise ratio (SNR) has been relatively high thus
% far. When the SNR is low, true spectral features are much harder to
% distinguish, and noise artifacts appear in spectral estimates based on
% the periodogram. The example below illustrates this:
fs = 1000; % Sampling frequency
t = (0:fs/10)/fs; % One-tenth second worth of samples
A = [1 2]; % Sinusoid amplitudes
f = [150;140]; % Sinusoid frequencies
xn = A*sin(2*pi*f*t) + 2*randn(size(t));
periodogram(xn,rectwin(length(xn)),1024,fs)
%%
% *Bias of the Periodogram*
%
% The periodogram is a biased estimator of the PSD. Its expected value was
% previously shown to be
%
% $$E\{ {\hat P_{xx}}(f)\} = {1 \over {{F_s}}}\int_{ - {F_s}/2}^{{F_s}/2}
% {{{{{\sin }^2}(L\pi (f - f')/{F_s})} \over {L{{\sin }^2}(\pi (f -
% f')/{F_s})}}} \,{P_{xx}}(f')\,df'.$$
%
% The periodogram is asymptotically unbiased, which is evident from the
% earlier observation that as the data record length tends to infinity, the
% frequency response of the rectangular window more closely approximates
% the Dirac delta function. However, in some cases the periodogram is a
% poor estimator of the PSD even when the data record is long. This is due
% to the variance of the periodogram, as explained below.
%%
% *Variance of the Periodogram*
%
% The variance of the periodogram can be shown to be
%
% $${\rm{Var}}({\hat P_{xx}}(f)) = \left\{ {\matrix{
% {P_{xx}^2(f),} & {0 < f < {{{F_s}} \mathord{\left/
% {\vphantom {{{F_s}} 2}} \right.
% \kern-\nulldelimiterspace} 2},} \cr
% {2P_{xx}^2(f),} & {f = 0\vee f = {{{F_s}} \mathord{\left/
% {\vphantom {{{F_s}} 2}} \right.
% \kern-\nulldelimiterspace} 2},} \cr
% } } \right.$$
%
% which indicates that the variance does not tend to zero as the data
% length $L$ tends to infinity. In statistical terms, the periodogram is
% not a consistent estimator of the PSD. Nevertheless, the periodogram can
% be a useful tool for spectral estimation in situations where the SNR is
% high, and especially if the data record is long.