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 1024 2048]; % Number of samples used for each case. % N = [4*1024 8*1024 16*1024 32*1024 64*1024 128*1024]; % Try me !!! NumOfRuns = 50; % Number of periodograms to be calculated for every different data record. h1 = figure('NumberTitle', 'off','Name', ... 'Figure 8.25 (1 of 2) Spectrum Estimation of a AR(4) Process Using Yule-Walker method',... 'Visible','off','Position', [100 0 800 800]); % 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 the 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)/2 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 Autocorrelation Method: p = 4; [aest, err] = acm(x,p); % 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(length(N)/2,2,2*m-1) plot(w1/pi,10*log10(Par),'k','LineWidth',0.1) title(['Overlay of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra using ', ... num2str(N(m)),' samples']) grid on; axis tight; ylim([-30 50]) xlabel('Frequency (units of pi)'); subplot(length(N)/2,2,2*m) plot(w1/pi,1/NumOfRuns*sum(10*log10(Par)),'r','LineWidth',1) %title(['Average of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra using ', ... % num2str(N(m)),' samples (red)']) title(['Average of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra (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'); ```

## Creation of Second Figure.

```h2 = figure('NumberTitle', 'off','Name', ... 'Figure 8.25 (2 of 2) Spectrum Estimation of a AR(4) Process Using Yule-Walker method', ... 'Visible','off','Position', [100 0 800 800]); set(0,'CurrentFigure',h2); for m=length(N)/2+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 Autocorrelation Method: p = 4; [aest, err] = acm(x,p); % 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(length(N)/2,2,2*m-1-length(N)) plot(w1/pi,10*log10(Par),'k','LineWidth',0.1) title(['Overlay of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra using ', ... num2str(N(m)),' samples']) grid on; axis tight; ylim([-30 50]) xlabel('Frequency (units of pi)'); subplot(length(N)/2,2,2*m-length(N)) plot(w1/pi,1/NumOfRuns*sum(10*log10(Par)),'r','LineWidth',1) % title(['Average of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra using ', ... % num2str(N(m)),' samples (red)']) title(['Average of ',num2str(NumOfRuns),' Yule-Walker AR(4) spectra (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 first figure. set(h2, 'Visible','on'); ```