MATLAB Examples

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

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