www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/harmest.m
function [Pxx,ar1,ar2] = harmest(y,maxlag,p_order,flag,nfft,norder)
%HARMEST Frequencies of harmonics in colored Gaussian noise.
% [Pxx,ar1,ar2] = harmest(y,maxlag,p_order,flag,nfft,norder)
% y - data matrix [nsamp x nrecs]
% maxlag - number of cumulant lags to compute [default = nsamp/12]
% p_order - order to use (dimension of signal subspace)
% user will be prompted if p_order <= 0
% flag - 'biased' or 'unbiased'
% nfft - fft length for spectra [default = 256]
% norder - cumulant order to use: 2 or 4 [default = 4]
%
% Pxx - a nfft x 7 matrix of spectral estimates
% ar1 - estimated parameters for the AR method.
% ar2 - estimated parameters for the min-norm method
%
% The seven columns of Pxx contain the spectral estimates based on
% the Eigenvector, Music, Pisarenko, ML, AR, periodogram methods,
% and the minimum-norm method.
%
% Copyright (c) 1991-98 by United Signals & Systems, Inc. and The MathWorks, Inc.
% $Revision: 1.6 $
% A. Swami January 20, 1993; August 8, 1994
% 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.
% --------- Parameter checks ---------------------------------
[nsamp,nrecs] = size(y);
if (nsamp == 1) nsamp = nrecs; nrecs = 1; y = y(:); end
if (exist('p_order') ~=1) p_order = 0; end
if (exist('flag') ~= 1) flag = 'biased'; end
if (flag(1) ~= 'b') flag = 'unbiased'; end
if (exist('maxlag') ~= 1) maxlag = nsamp/12; end
if (exist('nfft') ~= 1) nfft = 256; end
if (exist('norder') ~= 1) norder = 4; end
if (norder ~= 2 & norder ~= 4)
error('cumulant order - norder - should be 2 or 4')
end
Pxx = zeros(nfft,7);
% ------------ estimate second/fourth-order cumulants --------------
if (norder == 4)
M = maxlag;
cum_4y = cumest(y,4,M,nsamp,0,flag,0,0);
cum_4y = (cum_4y(M+1:2*M+1) + conj(cum_4y(M+1:-1:1)) )/2;
cum_4y = cum_4y/cum_4y(1);
Amat = toeplitz(conj(cum_4y),cum_4y);
else
Amat = zeros(maxlag+1,maxlag+1);
ind = [0:-1:-maxlag];
for n=maxlag+1:nsamp
xf = y(n+ind,:);
xb = conj(flipud(xf));
Amat = Amat + xf * xf' + xb * xb';
end
mu = mean(mean(y));
Amat = Amat - mu*mu'; % remove mean
Amat = Amat / ( nrecs*(nsamp-maxlag) );
end
[U,S,V] = svd(Amat); sval = diag(S);
% ----------- how many harmonics ? -------------------------
hold off, clf,
set(gcf,'Name', ...
['Hosa HARMEST - cum-order=', int2str(norder)] )
ls = length(sval);
plot(1:ls, sval,'-', 1:ls,sval,'go'), grid on
title(['sval: cum-',int2str(norder)])
drawnow
p = p_order;
while (p <= 0 | p > ls )
p = input([' enter order to use [1,' int2str(ls) '] ---> ']);
if (isempty(p)) p = 0; end
end
% ---------- Noise subspace methods --------------
M = maxlag+1;
nvec = M - p; mfft = nfft/2;
wte = ones(nvec,1) ./ sval(p+1:M);
Pvm = V(:,p+1:M) * V(:,p+1:M)';
Pve = V(:,p+1:M) * diag(wte) * V(:,p+1:M)';
psum = zeros(nfft,2);
psum(1,1) = sum(diag(Pve));
psum(1,2) = sum(diag(Pvm));
for k=1:maxlag
psum(k+1,1) = sum(diag(Pve,k));
psum(nfft-k+1,1) = sum(diag(Pve,-k));
psum(k+1,2) = sum(diag(Pvm,k));
psum(nfft-k+1,2) = sum(diag(Pvm,-k));
end
Pxx(:,1:2) = ones(nfft,2) ./ real(fft(psum,nfft));
[u1,s1,v1] = svd(Amat(1:p+1,1:p+1)); v1 = v1(:,p+1);
Pxx(:,3) = ones(nfft,1) ./ abs(fft(conj(v1),nfft)).^2;
% ------------ signal subspace method -----------------------------
% mlcapon
for k=1:p
Pxx(:,4) = Pxx(:,4) + abs(fft(conj(V(:,k))/sqrt(sval(k)),nfft)) .^ 2;
end
% --------- AR method -----------------------------------------------
% --------- SVD low-rank approximation (a la Cadzow) ----------------
Amat = U(:,1:p) * diag(sval(1:p)) * (V(:,1:p))';
Ahat = [];
[m,n] = size(U);
for k=p+1:n
Ahat = [Ahat; Amat(:,k-p:k)];
end
% force unity modulus solution
avec = [1; -Ahat(:,2:p+1) \ Ahat(:,1)];
avec = conj(avec);
if (exist('debug') == 1)
ncol = floor((p+2)/2);
nflip = floor((p+1)/2);
Ahat(:,1:nflip) = Ahat(:,1:nflip) + Ahat(:,p+1:-1:p+2-nflip);
avec = Ahat(:,2:ncol) \ Ahat(:,1);
avec = [1; -avec];
avec = [avec; avec(nflip:-1:1)];
end
ar1 = avec;
Pxx(:,5) = ones(nfft,1) ./ abs(fft(avec,nfft)) .^2;
% ----------- The periodogram estimate ----------------------
ind = [1:nsamp]';
for k=1:nrecs
ys = y(ind);
ys = ys - mean(ys);
Yf = fft(ys,nfft) / nsamp;
Pxx(:,6) = Pxx(:,6) + Yf .* conj(Yf);
ind = ind + nsamp;
end
Pxx(:,6) = Pxx(:,6) / nrecs;
Pxx(1,6) = Pxx(2,6); % dynamic range problems, else
% ----------------- Minimum-norm method -------------------
gvec = U(1,1:p).';
gmat = U(2:length(sval),1:p);
ar2 = [1; - conj(gmat) * gvec / (1 - gvec' * gvec)];
Pxx(:,7) = fft(ar2, nfft);
Pxx(:,7) = ones(nfft,1) ./ ( Pxx(:,7) .* conj(Pxx(:,7)) ) ;
% ----------- Display estimates ---------------------
waxis = [-mfft:mfft-1]/nfft;
Pxx = Pxx([mfft+1:nfft,1:mfft],:);
% - scale to max abs of unity for plots only
spmax = max(Pxx);
spmax = ones(1,length(spmax)) ./ spmax;
clf,
subplot(421), plot(1:ls, sval,'-', 1:ls,sval,'g.'), grid on
title(['sval: cum-',int2str(norder)])
subplot(422), plot(waxis,10*log10(Pxx(:,6)*spmax(6))), title('pxx'),grid on
subplot(423), plot(waxis,10*log10(Pxx(:,1)*spmax(1))), ylabel('eig'),grid on
subplot(424), plot(waxis,10*log10(Pxx(:,2)*spmax(2))), ylabel('music'),grid on
subplot(425), plot(waxis,10*log10(Pxx(:,3)*spmax(3))), ylabel('pisar'),grid on
subplot(426), plot(waxis,10*log10(Pxx(:,4)*spmax(4))), ylabel('ml'),grid on
subplot(427), plot(waxis,10*log10(Pxx(:,5)*spmax(5))), ylabel('ar'),grid on
subplot(428), plot(waxis,10*log10(Pxx(:,7)*spmax(7))), ylabel('min-norm'),
grid on
return