MATLAB Examples

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

Contents

Case I: Using the Autocorrelation (Yule-Walker) Method.

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');