MATLAB Examples

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