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 Blackman-Tukey 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 w1] = freqz(b,a,512); 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); % 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 Blackman - Tukey Method: win = 1; M = 4; Pma(k,:) = per_smooth(x, win, M); end subplot(4,2,2*m-1) plot(w1/pi,10*log10(Pma(:,1:L/2)),'k','LineWidth',0.1) title(['Overlay of ',num2str(NumOfRuns),' MA(4) spectra by Blackman-Tukey and ', ... num2str(N(m)),' samples']) grid on; axis tight; ylim([-40 30]) xlabel('Frequency (units of pi)'); subplot(4,2,2*m) plot(w1/pi,1/NumOfRuns*sum(10*log10(Pma(:,1:L/2))),'r','LineWidth',1) title(['Average of ',num2str(NumOfRuns),' MA(4) spectra by Blackman-Tukey (red)']) % title(['Average of ',num2str(NumOfRuns),' MA(4) spectra by Blackman-Tukey and ', ... % num2str(N(m)),' samples (red)']) axis tight; ylim([-30 20]) 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'); ```