MATLAB Examples

Example of par. 8.5.2. MA Spectrum Estimation Using Parametric Methods.

Contents

Case II: Using the Durbin's Method, section 4.7.3.

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.29 Spectrum Estimation of an MA(4) Process Using the Durbin''s method', ...
          'Visible','off','Position', [100 0 800 950]);

% First calculate and plot the true spectrum of the given process:
z = [0.98*exp(1i*0.2*pi)  0.98*exp(-1i*0.2*pi)  0.98*exp(1i*0.5*pi)  0.98*exp(-1i*0.5*pi)].';
p = [];
k = 1;
% Calculate the spectrum as the frequency response of given filter:
[b,a] = zp2tf(z,p,k);
[H omega] = freqz(b,a,1024);

L = 1024;
Pma = zeros(NumOfRuns, L);

for m=1:length(N)

     for k=1:NumOfRuns

          w = randn(1,N(m));
          x = zeros(1,N(m));
          n = 0:N(m)-1;
% Compute N samples of the MA(4) process x[n] = w[n] - 1.5857*w[n-1] + 1.9208*w[n-2] + ...
%                                                                      -1.5229*w[n-3] + 0.9224*w[n-4].
          x(1) = b(1)*w(1);
          x(2) = b(1)*w(2) + b(2)*w(1);
          x(3) = b(1)*w(3) + b(2)*w(2) + b(3)*w(1);
          x(4) = b(1)*w(4) + b(2)*w(3) + b(3)*w(2) + b(4)*w(1);
          for k1=5:N(m)
               x(k1) = b(1)*w(k1) + b(2)*w(k1-1) + b(3)*w(k1-2) + b(4)*w(k1-3) + b(5)*w(k1-4) ;
          end

% Estimate the MA(4) Power Spectrum using the Durbin's Method:
         q = 4; p = 32;
         b1 = durbin(x,p,q);

     % Form digital frequencies:
         omega1 = [ones(1,L); exp(-1i*omega.'); exp(-1i*2*omega.'); exp(-1i*3*omega.'); ...
                                               exp(-1i*4*omega.');];

     % Calculate MA(4) Power Spectrum:
         Pma(k,:) = abs(b1.'*omega1).^2;  % Eq. (8.133).

    end

    subplot(4,2,2*m-1)
    plot(omega/pi,10*log10(Pma),'k','LineWidth',0.1)
    title(['Overlay of ',num2str(NumOfRuns),' MA(4) spectra by Durbin''s Method and ', ...
                                 num2str(N(m)),' samples'])
    grid on;
    axis tight;
    ylim([-25 20])
    xlabel('Frequency (units of pi)');

    subplot(4,2,2*m)
    plot(omega/pi,1/NumOfRuns*sum(10*log10(Pma)),'r','LineWidth',1)
    title(['Average of ',num2str(NumOfRuns),' MA(4) spectra by Durbin''s Method (red)'])
    % title(['Average of ',num2str(NumOfRuns),' MA(4) spectra by Durbin''s Method and ', ...
    %                                num2str(N(m)),' samples (red)'])
    axis tight;
    ylim([-30 20])
    grid on;    hold on;
    plot(omega/pi,20*log10(abs(H)))
    xlabel('Frequency (units of pi)');

end

% Restore the visibility of the figure.
set(h1, 'Visible','on');