# Example 8.6.5. A Comparison of Frequency Estimation Methods.

## Workspace Initialization.

```clc; clear; close all;
```

## Signal Definition

```N = [64 128 256 512]; % Number of samples used for each case.

NumOfRuns = 10;    % Number of frequency estimation functions to be calculated for every different data record.

% Open a new figure:
h1 = figure('NumberTitle', 'off','Name', ...
'Figure 8.37 Frequency Estimation for a Process Consisting of 4 Complex Exponentials in WGN', ...
'Visible','off','Position', [100 0 800 950]);

sigmaw = sqrt(0.5);
% Create frequency axis:
L = 1024;
w = 0:2*pi/L:2*pi*(1-1/L);

Pmusic = zeros(NumOfRuns,L);

for m=1:length(N)
for k=1:NumOfRuns
```
```           noise = sigmaw*randn(1,N(m));
n = 0:N(m)-1;
phi = 2*pi*(rand(4,1) - 0.5);
omegas = [0.2*pi 0.3*pi 0.8*pi 1.2*pi].';

% Compute N(m) samples of the given Harmonic Process:
x = sum(exp(1i*omegas*n + phi*ones(1,N(m)))) + noise;
```

## Estimate the Power Spectrum using the MUSIC algorithm:

```        p = 4;   % Number of complex exponentials present assumed known.
M = N(m); % Size of autocorrelation matrix to be used in the following methods:
Pmusic(k,:) = music(x,p,M).';
```
```    end

subplot(4,1,m)
plot(w/pi,Pmusic,'k','LineWidth',0.1)
title(['Overlay plot of ',num2str(NumOfRuns),' MUSIC Spectra Using ',num2str(N(m)),' samples'])
grid on;
axis tight;
ylim([-45 10])

end

xlabel('Frequency (units of pi)');

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