www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/biceps.m
function [hest,ceps,A,B,minh,maxh] = biceps(y,p,q,nsamp,overlap,flag,lh)
%BICEPS Non-parametric system identification using the bicepstrum.
% BICEPS estimates the complex cepstrum, using third-order
% cumulants, and then reconstructs the impulse response.
%
% [hest,ceps,A,B,minh,maxh] = biceps(y,p,q,sampseg,overlap,flag,lh)
% y - data vector or matrix
% p,q - cepstral orders (truncation points) [MUST be specified]
% sampseg - samples per record [default = row size]
% overlap - percentage overlap of records [default = 0; max=99]
% flag - 'biased' or 'unbiased' [default = 'biased']
% lh - hest indices will range from -lh to lh [default = 2(p+q)]
%
% if y is a matrix, each column is assumed to be a different realization
% or record; in this case, overlap is set to 0, and record size is set
% to the row length of the matrix.
%
% hest - estimated impulse response, h(n): n= -lh, ... ,lh
% ceps - estimated complex cepstrum
% A - estimated A(k)'s, k=1,...,p
% B - estimated B(k)'s, k=1,...,q
% minh - minimum phase component of hest; length lh+1
% maxh - maximum phase component of hest; length lh+1
% Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved.
% $Revision: 1.7 $
% A. Swami January 20, 1993. Revised Jan 20, 1995.
% 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 -----------------------
[ly, nrecord] = size(y);
if (ly == 1) y = y(:); ly = nrecord; nrecord = 1; end
if (exist('p') ~=1) error('order p must be specified'); end
if (exist('q') ~=1) error('order q must be specified'); end
if (nrecord > 1) nsamp = ly; overlap = 0; end
if (exist('nsamp') ~= 1) nsamp = 0; end
if (nrecord > 1) nsamp = ly; end
if (exist('overlap') ~= 1) overlap = 0; end
if (nrecord > 1) overlap = 0; end
if (exist('flag') ~= 1) flag = 'biased'; end
if (exist('lh') ~= 1) lh = 2*(p+q); end
overlap = min(99,(max(overlap,0)));
nsamp = min(ly,max(nsamp,0));
% ------ estimate cumulants ---------------------------------------
w = max(p,q);
w = w + rem(w,2);
wby2 = w/2;
cum_y = zeros(w+p+q+1,4*w+1);
M = 2 * w;
n_loc = wby2+w+1;
for k = -wby2-w:wby2+w
cum_y(n_loc+k,:) = cumest(y(:),3,M,nsamp,overlap,flag,k).';
end
% -------- set up system of linear equations ----------------------
% \sum_{i=-q}^{p} i ch(i) [ C(m-i,n) - C(m+1,n+i) ] = m C(m,n)
% with, |m| <= max(p,q) and |n| < max(p,q)/2
xorg = M + 1;
yorg = wby2 + w + 1; % c3(0,0) is in cum_y(xorg,yorg)
ind1 = (-w+q:w+q) + xorg;
ind2 = (-w+q:-1:-w-p) + xorg;
ind3 = yorg + (q:-1:-p);
ind4 = xorg + (-w:w);
Amat = [];
for n = -wby2:wby2
alpha = cum_y(yorg+n,:); % picks up (.,n) slice
tmp = toeplitz(alpha(ind1), alpha(ind2)); % I term b(m-k,n) m=-w:w
tmp1 = cum_y(ind3-n,ind4-n).';
Amat = [Amat; tmp-tmp1];
end
% ------ solve and convert to complex cepstral coefficients -----------
Amat = Amat(:,[1:q,q+2:p+q+1]);
rvec = cum_y(yorg+ (-wby2:wby2), ind4) * diag(-w:w);
rvec = rvec';
rvec = rvec(:);
dcepstrum = Amat\rvec;
cepstrum = dcepstrum ./ [-q:-1,1:p]';
ceps = [cepstrum(1:q); 0; cepstrum(q+1:q+p)];
% ------- estimate causal/anti-causal IR's directly:
maxh = real( ifft( exp( fft( [1;cepstrum(q:-1:1) ], lh+1 ) )));
minh = real( ifft( exp( fft( [1;cepstrum(q+1:q+p)], lh+1 ) )));
hest = conv(flipud(maxh),minh);
A = - dcepstrum(q+1:q+p);
B = dcepstrum(q:-1:1);
% ------- or estimate impulse response using Oppenheim and Schafer's method
opschaf = 0;
if (opschaf == 1)
cpos = cepstrum(q+1:q+p) .* [1:p]';
cneg = cepstrum(q:-1:1) .* [1:q]';
ir_len = lh;
cpos = [cpos; zeros(ir_len-p,1)];
cneg = [cneg; zeros(ir_len-q,1)];
ha = zeros(ir_len+1,1); hc = ha; ha(1) = 1; hc(1) = 1;
for n=1:ir_len
ha(n+1) = cneg(1:n).' * ha(n:-1:1) /n;
hc(n+1) = cpos(1:n).' * hc(n:-1:1) /n;
end
ha = ha(ir_len+1:-1:1);
hest = conv(ha,hc);
minh = hc; maxh = ha;
hest = hest/hest(2*(p+q)+1);
A = -cpos(1:p);
B = -cneg(1:q);
end
% --------- display estimates --------------------
clf, subplot(211)
plot([-q:p]', ceps), grid on
title('complex cepstrum')
xlabel('sample number')
subplot(212)
taxis = [-lh:lh]';
plot(taxis, hest), grid on
title('impulse response ')
xlabel('sample number')
set (gcf, 'Name','Hosa BICEPS')
return