MATLAB Examples

# Example of par. 8.5.1. AR Spectrum Estimation Using Parametric Methods.

## Workspace Initialization.

```clc; clear; close all; ```

## Signal Definition

```N = [64 128 256 512]; % Number of samples used for each case. NumOfRuns = 50; % Number of periodograms to be calculated for every different data record. h1 = figure('NumberTitle', 'off','Name', ... 'Figure 8.28 Spectrum Estimation of an AR(4) Process Using the Burg''s method', ... 'Visible','off','Position', [100 0 800 950]); % First calculate and plot the true spectrum of the given process: z = []; p = [0.98*exp(1i*0.2*pi) 0.98*exp(-1i*0.2*pi) 0.98*exp(1i*0.3*pi) 0.98*exp(-1i*0.3*pi)]; k = 1; % Calculate the spectrum as the frequency response of given filter: [b,a] = zp2tf(z,p,k); [H w1] = freqz(b,a,1024); Par = zeros(NumOfRuns, 1024); % Calculate the Spectrum from Process Data using the Yule-Walker method: for m=1:length(N) for k=1:NumOfRuns w = randn(1,N(m)); x = zeros(1,N(m)); n = 0:N(m); % Compute N samples of the AR(4) process x[n] = 2.7377*x[n-1] - 3.7476*x[n-2] + ... % 2.6293*x[n-3] - 0.9224*x[n-4] + w[n]. x(1) = b(5)*w(1); x(2) = -a(2)*x(1) + b(5)*w(2); x(3) = -a(2)*x(2) - a(3)*x(1) + b(5)*w(3); x(4) = -a(2)*x(3) - a(3)*x(2) - a(4)*x(1) + b(5)*w(4); for k1=5:N(m) x(k1) = -a(2)*x(k1-1) - a(3)*x(k1-2) - a(4)*x(k1-3) - a(5)*x(k1-4) + b(5)*w(k1); end % Estimate the AR(4) Coefficients using the Burg's Method: p = 4; % Calculate the gamma coefficients... [gamma,err] = burg(x,p); % Then calculate the IIR filter coefficients using the Step-up Recursion: aest = gtoa(gamma); % Calculate a small part of the ACF of x to be used for b0 estimation: maxlag = 10; [rx0 lags] = xcorr(x,maxlag); Rx = rx0/N(m); % |b0|^2 = Rx(0) + a(1)*Rx(1) + a(2)*Rx(2) + % a(3)*Rx(3) + a(4)*Rx(4) b0est_squared = Rx(maxlag+1) + aest(2)*Rx(maxlag+2) + aest(3)*Rx(maxlag+3) + ... aest(4)*Rx(maxlag+4) + aest(5)*Rx(maxlag+5); % The estimated AR(4) spectrum is: Par(k,:) = b0est_squared./(abs(1 + aest(2)*exp(-1i*w1) + aest(3)*exp(-1i*2*w1) + ... aest(4)*exp(-1i*3*w1) + aest(5)*exp(-1i*4*w1))).^2; end subplot(4,2,2*m-1) plot(w1/pi,10*log10(abs(Par)),'k','LineWidth',0.1) title(['Overlay plot of ',num2str(NumOfRuns),' AR(4) spectra by Burg''s Method and ', ... num2str(N(m)),' samples']) grid on; axis tight; ylim([-30 60]) xlabel('Frequency (units of pi)'); subplot(4,2,2*m) plot(w1/pi,1/NumOfRuns*sum(10*log10(abs(Par))),'r','LineWidth',1) title(['Average of ',num2str(NumOfRuns),' AR(4) spectra by Burg''s Method (red)']) % title(['Average of ',num2str(NumOfRuns),' AR(4) spectra by Burg''s Method and ', ... % num2str(N(m)),' samples (red)']) axis tight; ylim([-30 50]) grid on; hold on; plot(w1/pi,20*log10(abs(H))) xlabel('Frequency (units of pi)'); end % Restore the visibility of the figure. set(h1, 'Visible','on'); ```