www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/bicoherx.m
function [bic,waxis] = bicoherx (w,x,y, nfft, wind, nsamp, overlap) %BICOHERX - Direct (FD) method for estimating cross-bicoherence % [bic,waxis] = bicoherx (w,x,y, nfft, wind, segsamp, overlap) % w,x,y - data vector or time-series % - should have identical dimensions % nfft - fft length [default = power of two > nsamp] % actual size used is power of two greater than 'nsamp' % wind - specifies the time-domain window to be applied to each % data segment; should be of length 'segsamp' (see below); % otherwise, the default Hanning window is used. % segsamp - samples per segment [default: such that we have 8 segments] % - if x is a matrix, segsamp is set to the number of rows % overlap - percentage overlap, 0 to 99 [default = 50] % - if y is a matrix, overlap is set to 0. % bic - estimated cross-bicoherence: an nfft x nfft array, with % origin at center, and axes pointing down and to the right. % waxis - vector of frequencies associated with the rows and columns % of bic; sampling frequency is assumed to be 1. % Copyright (c) 1991-1999 by United Signals & Systems, Inc. and The Mathworks, Inc. All Rights Reserved. % $Revision: 1.7 $ % A. Swami January 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 ----------------------------- if (size(w) ~= size(x) | size(x) ~= size(y) ) error(' w, x, and y should have identical dimensions') end [ly, nrecs] = size(y); if (ly == 1) ly = nrecs; nrecs = 1; w = w(:); x = x(:); y = y(:); end if (exist('nfft') ~= 1) nfft = 128; end if (exist('overlap') ~= 1) overlap = 50; end overlap = max(0,min(overlap,99)); if (nrecs > 1) overlap = 0; end if (exist('nsamp') ~= 1) nsamp = 0; end if (nrecs > 1) nsamp = ly; end if (nrecs == 1 & nsamp <= 0) nsamp = fix(ly/ (8 - 7 * overlap/100)); end if (nfft < nsamp) nfft = 2^nextpow2(nsamp); end overlap = fix(overlap/100 * nsamp); nadvance = nsamp - overlap; nrecs = fix ( (ly*nrecs - overlap) / nadvance); % ---------------------------------------------------------------------- if (exist('wind') ~= 1) wind = hanning(nsamp); end [rw,cw] = size(wind); if (min(rw,cw) ~= 1 | max(rw,cw) ~= nsamp) disp(['Segment size is ',int2str(nsamp)]) disp(['"wind" array is ',int2str(rw),' by ',int2str(cw)]) disp(['Using default window']) wind = hanning(nsamp); end wind = wind(:); % ---------------- accumulate triple products ---------------------- bic = zeros(nfft,nfft); Pyy = zeros(nfft,1); Pww = Pyy; Pxx = Pyy; mask = hankel([1:nfft],[nfft,1:nfft-1] ); % the hankel mask (faster) Yf12 = zeros(nfft,nfft); ind = [1:nsamp]'; for k = 1:nrecs ws = w(ind); ws = (ws-mean(ws)).* wind; Wf = fft(ws,nfft) / nsamp; CWf = conj(Wf); Pww = Pww + Wf .* CWf; xs = x(ind); xs = (xs-mean(xs)).* wind; Xf = fft(xs,nfft) / nsamp; CXf = conj(Xf); Pxx = Pxx + Xf .* CXf; ys = y(ind); ys = (ys(:) - mean(ys)) .* wind; Yf = fft(ys,nfft) / nsamp; CYf = conj(Yf); Pyy = Pyy + Yf .* CYf; Yf12(:) = CYf(mask); bic = bic + (Wf * Xf.') .* Yf12; ind = ind + nadvance; end bic = bic / nrecs; Pyy = Pyy / nrecs; Pww = Pww / nrecs; Pxx = Pxx / nrecs; mask(:) = Pyy(mask); bic = abs(bic).^2 ./ (Pww * Pxx.' .* mask); bic = fftshift(bic) ; % ------------ contour plot of magnitude bispectum -------------------- if (rem(nfft,2) == 0) waxis = [-nfft/2:(nfft/2-1)]'/nfft; else waxis = [-(nfft-1)/2:(nfft-1)/2]'/nfft; end hold off, clf % contour(bic,4,waxis,waxis),grid contour(waxis,waxis,bic,4), grid on title('Cross-Bicoherence') xlabel('f1'), ylabel('f2') set(gcf,'Name','Hosa BICOHERX') [colmax,row] = max(bic) ; [maxval,col] = max(colmax); row = row(col); disp(['Max: bic(',num2str(waxis(row)),',',num2str(waxis(col)),') = ', ... num2str(maxval) ]) return