MATLAB Examples

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

Contents

Case II: Using the Covariance Method.

Workspace Initialization.

clc; clear; close all;

Signal Definition

N = [64 128 256 512 1024 2048]; % 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.26 (1 of 2) Spectrum Estimation of an AR(4) Process Using the Covariance method',...
                    'Visible','off','Position', [100 0 800 900]);

% 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)/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 Covariance Method:
        p = 4;
        [aest, err] = covm(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),' AR(4) Spectra by Covariance Method and ', ...
                                 num2str(N(m)),' samples'])
    grid on;
    axis tight;
    ylim([-30 60])
    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),' AR(4) Spectra by Covariance Method (red)']);
%     title(['Average of ',num2str(NumOfRuns),' AR(4) Spectra by Covariance 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');

Creation of Second Figure.

h2 = figure('NumberTitle', 'off','Name', ...
             'Figure 8.26 Spectrum Estimation of an AR(4) Process by the Covariance method', ...
             'Visible','off','Position', [100 0 800 900]);
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 Covariance Method:
        p = 4;
        [aest, err] = covm(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),' AR(4) Spectra by Covariance Method and ', ...
                                 num2str(N(m)),' samples'])
    grid on;
    axis tight;
    ylim([-30 60])
    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),' AR(4) Spectra by Covariance Method (red)']);
%     title(['Average of ',num2str(NumOfRuns),' AR(4) Spectra by Covariance 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(h2, 'Visible','on');