MATLAB Examples

Example 8.6.5. A Comparison of Frequency Estimation Methods.

Contents

Case I: Pisarenko Harmonic Decomposition (PHD).

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

Pphd = 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 Pisarenko Harmonic Decomposition (PHD):

           p = 4; % Number of complex exponentials present assumed known.
           [vmin, sigma] = phd(x,p);

           for k1=1:L
                e = exp(-1i*w(k1)*(0:p));
                Pphd(k,k1) = 1./abs(e*vmin).^2;
           end
           Pphd(k,:) = 10*log10(Pphd(k,:));
     end

     subplot(4,1,m)
     plot(w/pi,Pphd,'k','LineWidth',0.1)
     title(['Overlay of ',num2str(NumOfRuns),' PHD Spectra Using',num2str(N(m)),' samples'])
     grid on;
     axis tight;
     ylim([-10 100])
     xlabel('Frequency (units of pi)');

end

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