MATLAB Examples

Example 8.6.5. A Comparison of Frequency Estimation Methods.

Contents

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 maximized figure:
scrsz = get(0,'ScreenSize');
hfig1 = figure('NumberTitle', 'off','Name', ...
          'Figure 8.37 (1 of 2) Frequency estimation for a process made of four complex exponentials in WGN', ...
                     'Position', [0 0 800 900]);
hfig2 = figure('NumberTitle', 'off','Name', ...
          'Figure 8.37 (2 of 2) Frequency estimation for a process made of four complex exponentials in WGN', ...
                     'Position', [800 0 800 900]);

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

Pphd            = zeros(NumOfRuns,L);
Pmusic        = zeros(NumOfRuns,L);
Pmin_norm = zeros(NumOfRuns,L);
Pev              = zeros(NumOfRuns,L);

for m=1:length(N)
     for k=1:NumOfRuns
          noise = sigmaw*randn(1,N(m));
   %            x = zeros(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,:));

Estimate the Power Spectrum using the MUSIC algorithm:

        M = N(m); % Size of autocorrelation matrix to be used in the following methods:
        Pmusic(k,:) = music(x,p,M).';

Estimate the Power Spectrum using the Minimum Norm algorithm:

        Pmin_norm(k,:) = min_norm(x,p,M).';

Estimate the Power Spectrum using the Eigenvector method:

        Pev(k,:) = ev(x,p,M).';
     end

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

     subplot(4,2,2*m)
     plot(w/pi,Pmusic,'k','LineWidth',0.1)
     title(['Overlay plot of ',num2str(NumOfRuns),' MUSIC spectra and ',num2str(N(m)),' samples'])
     grid on;
     axis tight;
     ylim([-45 10])
    % xlabel('Frequency (units of pi)');

    figure(hfig2);
    subplot(4,2,2*m-1)
    plot(w/pi,Pmin_norm,'k','LineWidth',0.1)
    title(['Overlay plot of ',num2str(NumOfRuns),' Minimum Norm spectra and ',num2str(N(m)), ' samples'])
    grid on;
    axis tight;
    ylim([-5 30])
    % xlabel('Frequency (units of pi)');

    subplot(4,2,2*m)
    plot(w/pi,Pev,'k','LineWidth',0.1)
    title(['Overlay plot of ',num2str(NumOfRuns),' Eigenvector spectra and ',num2str(N(m)), ' samples'])
    grid on;
    axis tight;
    ylim([-40 40])
    xlabel('Frequency (units of pi)');

end


% If you want to plot on a single figure:
%
%      subplot(4,4,4*m-1)
%      plot(w/pi,Pphd,'k','LineWidth',0.1)
%      title(['Overlay plot of ',num2str(NumOfRuns),' PHD spectra and ',num2str(N(m)),' samples'])
%      grid on;
%      axis tight;
%      ylim([-10 100])
%      xlabel('Frequency (units of pi)');
%
%      subplot(4,4,4*m-2)
%      plot(w/pi,Pmusic,'k','LineWidth',0.1)
%      title(['Overlay plot of ',num2str(NumOfRuns),' MUSIC spectra and ',num2str(N(m)),' samples'])
%      grid on;
%      axis tight;
%      ylim([-45 10])
%      xlabel('Frequency (units of pi)');
%
%     subplot(4,4,4*m-3)
%     plot(w/pi,Pmin_norm,'k','LineWidth',0.1)
%     title(['Overlay plot of ',num2str(NumOfRuns),' Minimum Norm spectra and ',num2str(N(m)), ...
%             ' samples'])
%     grid on;
%     axis tight;
%     ylim([-5 30])
%     xlabel('Frequency (units of pi)');
%
%     subplot(4,4,4*m)
%     plot(w/pi,Pev,'k','LineWidth',0.1)
%     title(['Overlay plot of ',num2str(NumOfRuns),' Eigenvector spectra and ',num2str(N(m)), ...
%             ' samples'])
%     grid on;
%     axis tight;
%     ylim([-40 40])
%     xlabel('Frequency (units of pi)');
%     tightfig(hfig);