www.gusucode.com > signal 案例源码程序 matlab代码 > signal/BiasAndVariabilityInThePeriodogramExample.m
%% Bias and Variability in the Periodogram % This example shows how to reduce bias and variability in the periodogram. % Using a window can reduce the bias in the periodogram, and using windows % with averaging can reduce variability. % Copyright 2015 The MathWorks, Inc. %% % Use wide-sense stationary autoregressive (AR) processes to show the % effects of bias and variability in the periodogram. AR processes present % a convenient model because their PSDs have closed-form expressions. % Create an AR(2) model of the following form: % % $$y(n)-0.75y(n-1)+0.5y(n-2)=\varepsilon(n),$$ % % where $\varepsilon(n)$ is a zero mean white noise sequence with some % specified variance. In this example, assume the variance and the sampling % period to be 1. To simulate the preceding AR(2) process, create an % all-pole (IIR) filter. View the filter's magnitude response. B2 = 1; A2 = [1 -0.75 0.5]; fvtool(B2,A2) %% % This process is bandpass. The dynamic range of the PSD is approximately % 14.5 dB, as you can determine with the following code. [H2,W2] = freqz(B2,A2,1e3,1); dr2 = max(20*log10(abs(H2)))-min(20*log10(abs(H2))) %% % By examining the placement of the poles, you see that this AR(2) process % is stable. The two poles are inside the unit circle. fvtool(B2,A2,'Analysis','polezero') %% % Next, create an AR(4) process described by the following equation: % % $$y(n)-2.7607y(n-1)+3.8106y(n-2)-2.6535y(n-3)+ % 0.9238y(n-4)=\varepsilon(n).$$ % % Use the following code to view the magnitude response of this IIR system. B4 = 1; A4 = [1 -2.7607 3.8106 -2.6535 0.9238]; fvtool(B4,A4) %% % Examining the placement of the poles, you can see this AR(4) process is % also stable. The four poles are inside the unit circle. fvtool(B4,A4,'Analysis','polezero') %% % The dynamic range of this PSD is approximately 65 dB, much larger than % the AR(2) model. [H4,W4] = freqz(B4,A4,1e3,1); dr4 = max(20*log10(abs(H4)))-min(20*log10(abs(H4))) %% % To simulate realizations from these AR(_p_) processes, use |randn| and % |filter|. Set the random number generator to the default settings to % produce repeatable results. Plot the realizations. rng default x = randn(1e3,1); y2 = filter(B2,A2,x); y4 = filter(B4,A4,x); subplot(2,1,1) plot(y2) title('AR(2) Process') xlabel('Time') subplot(2,1,2) plot(y4) title('AR(4) Process') xlabel('Time') %% % Compute and plot the periodograms of the AR(2) and AR(4) realizations. % Compare the results against the true PSD. Note that |periodogram| % converts the frequencies to millihertz for plotting. Fs = 1; NFFT = length(y2); subplot(2,1,1) periodogram(y2,rectwin(NFFT),NFFT,Fs) hold on plot(1000*W2,20*log10(abs(H2)),'r','linewidth',2) title('AR(2) PSD and Periodogram') subplot(2,1,2) periodogram(y4,rectwin(NFFT),NFFT,Fs) hold on plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2) title('AR(4) PSD and Periodogram') text(350,20,'\downarrow Bias') %% % In the case of the AR(2) process, the periodogram estimate follows the % shape of the true PSD but exhibits considerable variability. This is due % to the low degrees of freedom. The pronounced negative deflections (in % dB) in the periodogram are explained by taking the log of a chi-square % random variable with two degrees of freedom. % % In the case of the AR(4) process, the periodogram follows the shape of % the true PSD at low frequencies but deviates from the PSD in the high % frequencies. This is the effect of the convolution with Fejer's kernel. % The large dynamic range of the AR(4) process compared to the AR(2) % process is what makes the bias more pronounced. % % Mitigate the bias demonstrated in the AR(4) process by using a taper, or % window. In this example, use a Hamming window to taper the AR(4) % realization before obtaining the periodogram. figure periodogram(y4,hamming(length(y4)),NFFT,Fs) hold on plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2) title('AR(4) PSD and Periodogram with Hamming Window') legend('Periodogram','AR(4) PSD') %% % Note that the periodogram estimate now follows the true AR(4) PSD over % the entire Nyquist frequency range. The periodogram estimates still only % have two degrees of freedom so the use of a window does not reduce the % variability of periodogram, but it does address bias. % % In nonparametric spectral estimation, two methods for increasing the % degrees of freedom and reducing the variability of the periodogram are % Welch's overlapped segment averaging and multitaper spectral estimation. % % Obtain a multitaper estimate of the AR(4) time series using a time half % bandwidth product of 3.5. Plot the result. NW = 3.5; figure pmtm(y4,NW,NFFT,Fs) hold on plot(1000*W4,20*log10(abs(H4)),'r','linewidth',2) legend('Multitaper Estimate','AR(4) PSD') %% % The multitaper method produces a PSD estimate with significantly less % variability than the periodogram. Because the multitaper method also uses % windows, you see that the bias of the periodogram is also addressed.