% DEMHARM  HOSA Toolbox Demo: Harmonic Retrieval in Colored Gaussian Noise
echo off

% demo of harmest (harmgen) 

% A. Swami April 15, 1993
% Copyright (c) 1991-2001 by United Signals & Systems, Inc. 
%       $Revision: 1.7 $

%     RESTRICTED RIGHTS LEGEND
% Use, duplication, or disclosure by the Government is subject to
% restrictions as set forth in subparagraph (c) (1) (ii) of the 
% Rights in Technical Data and Computer Software clause of DFARS
% 252.227-7013. 
% Manufacturer: United Signals & Systems, Inc., P.O. Box 2374, 
% Culver City, California 90231. 
%
%  This material may be reproduced by or for the U.S. Government pursuant 
%  to the copyright license under the clause at DFARS 252.227-7013. 

clear, clc, 
echo on

%           Harmonic Retrieval in Colored Gaussian Noise

% The problem of estimating the parameters of harmonic signals, observed
% in the presence of noise, is an ubiquitous one. 

% The HOSA Toolbox offers two routines:
% 1.  HARMGEN may be used to generate synthetics consisting of 
%     harmonics contaminated by colored Gaussian noise. 
%
% 2.  HARMEST may be used to estimate the number of harmonics,
%     and their frequencies.  This method uses either the correlation or 
%     fourth-order cumulants; in the latter case, the signal may be corrupted
%     by colored Gaussian noise, even if the noise autocorrelation is unknown.
%
% We will test HARMEST on some data generated by HARMGEN.
% The data consist of two unity amplitude harmonics, with frequencies 
% 0.1 Hz and 0.2 Hz, contaminated by colored Gaussian noise with a 
% variance of 0.5.  The broad-band SNR is 3dB.  The data matrix consists
% of 32 independent realizations, each with 128 samples;  the starting
% phases are chosen randomly from an uniform distribution. 
% Colored Gaussian noise was generated by passing white Gaussian noise
% through an AR  filter with coefficients [1, -1.058, 0.81];  the noise
% spectrum has a pole at 0.15 Hz with a damping factor of 0.9. 

% Hit any key to continue 
pause 

load harm

% HARMEST obtains the solution to a linear regression equation based
% on fourth-order cumulants;  the roots of the parameter vector so
% obtained yield the frequencies of the harmonics.  Additionally,
% the rank of the matrix equals twice the number of harmonics. 
% Thus, HARMEST may be used to determine the number of harmonics as well.
% We use the SVD of the cumulant matrix to determine the rank. 
%
% You will notice from the singular value plot that the singular values
% show a steep drop at p=4, indicating the presence of two real harmonics.
% We let p=4 in the demo: 

% Hit any key to continue 
pause 
 
 clf
 Pxx4 = harmest(zmat,12,4,'unbiased',256,4); 
 set (gcf, 'Name', 'HOSA cumulant-based HARMEST')

% HARMEST estimates "power spectra" 
% These "Power Spectra" are estimated using Eigenvector, MUSIC, Pisarenko 
% ML(Capon) and AR methods based on fourth-order cumulants. 
% The conventional periodogram estimate is also obtained. 

% Note that the conventional periodogram displays peaks at 0.1 Hz and 0.2 Hz, 
% as well as a broad peak at 0.15 Hz which is due to the Gaussian noise. 
% The fourth-order cumulant based estimates display two sharp peaks at the
% correct frequencies of 0.1 Hz and 0.2 Hz.  

% We can correctly estimate the number of harmonics, and their parameters
% using fourth-order cumulants, even when the data are corrupted by
% colored Gaussian noise.  Used in conjunction with second-order statistics
% (the periodogram in this case), we also obtain information about the 
% additive noise. 

% Now we will apply the corresponding correlation-based algorithms to 
% the same data.  You will notice that the singular value plot is now
% indicative of p=6, rather than p=4; hence, we use p=6 in the demo

% Hit any key to continue
pause

clf
Pxx2 = harmest(zmat,12,6,'unbiased',256,2); 
set (gcf, 'Name', 'HOSA correlation-based HARMEST')

% Hit any key to continue
pause

% Let us look at both sets of estimates together 
echo off 
[nfft,ncol] = size(Pxx4);
waxis = [-nfft/2:nfft/2-1]'/nfft; 
clf, 

for k=1:2 
   if   (k==1) Pxx = Pxx4; subplot(211)
   else        Pxx = Pxx2; subplot(212) 
   end 
   spmax = max(Pxx);
   spmax = ones(1,length(spmax)) ./ spmax; 
   semilogy(waxis, Pxx * diag(spmax),       ...
            waxis, Pxx(:,1)*spmax(1), '-',  ...
            waxis, Pxx(:,2)*spmax(2), '+',  ...
            waxis, Pxx(:,3)*spmax(3), '*',  ...
            waxis, Pxx(:,4)*spmax(4), 'o',  ...
            waxis, Pxx(:,5)*spmax(5), '--', ...
            waxis, Pxx(:,6)*spmax(6), 'x',  ... 
            waxis, Pxx(:,7)*spmax(7), '-.' )
   xlabel('Frequency')
   title('eig(-) music(+) pisar(*) ml(o) ar(--) pxx(x) min-norm(-.)')
   ylabel('"C4-Spectra"')
   if (k==2)    ylabel('"C2-Spectra"'), end 
   grid on 
end
echo on 

% Compare the correlation based estimates (bottom panel) 
% with the cumulant-based estimates (top panel). 
% 
% Hit any key to return to the main menu. 
pause 
echo off
clc