MATLAB Examples

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

Contents

Case III: Using the Modified Covariance Method.

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.27 Spectrum Estimation of an AR(4) Process by the Modified Covariance 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 Modified Covariance Method:
        p = 4;
        [aest, err] = mcov(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(4,2,2*m-1)
    plot(w1/pi,10*log10(abs(Par)),'k','LineWidth',0.1)
    title(['Overlay of ',num2str(NumOfRuns),' AR(4) spectra by Mod. Covariance 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 Mod. Covariance (red)'])
%     title(['Average of ',num2str(NumOfRuns),' AR(4) spectra by Mod. Covariance 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');