www.gusucode.com > MATLAB高阶累积量工具箱 > MATLAB高阶累积量工具箱/MATLAB高阶累积量工具箱/hosa/bispecd.m
function [Bspec,waxis] = bispecd (y, nfft, wind, nsamp, overlap) %BISPECD Bispectrum estimation using the direct (fft-based) approach. % [Bspec,waxis] = bispecd (y, nfft, wind, segsamp, overlap) % y - data vector or time-series % nfft - fft length [default = power of two > segsamp] % wind - window specification for frequency-domain smoothing % if 'wind' is a scalar, it specifies the length of the side % of the square for the Rao-Gabr optimal window [default=5] % if 'wind' is a vector, a 2D window will be calculated via % w2(i,j) = wind(i) * wind(j) * wind(i+j) % if 'wind' is a matrix, it specifies the 2-D filter directly % segsamp - samples per segment [default: such that we have 8 segments] % - if y is a matrix, segsamp is set to the number of rows % overlap - percentage overlap [default = 50] % - if y is a matrix, overlap is set to 0. % % Bspec - estimated bispectrum: an nfft x nfft array, with origin % at the center, and axes pointing down and to the right. % waxis - vector of frequencies associated with the rows and columns % of Bspec; 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.8 $ % A. Swami January 20, 1993. % 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, nrecs] = size(y); if (ly == 1) y = y(:); ly = nrecs; nrecs = 1; end if (exist('nfft') ~= 1) nfft = 128; end if (exist('overlap') ~= 1) overlap = 50; end overlap = min(99,max(overlap,0)); 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(nsamp * overlap / 100); % added 2/14 nadvance = nsamp - overlap; nrecs = fix ( (ly*nrecs - overlap) / nadvance); % ------------------- create the 2-D window ------------------------- if (exist('wind') ~= 1) wind = 5; end [m,n] = size(wind); window = wind; if (max(m,n) == 1) % scalar: wind is size of Rao-Gabr window winsize = wind; if (winsize < 0) winsize = 5; end % the window length L winsize = winsize - rem(winsize,2) + 1; % make it odd if (winsize > 1) mwind = fix (nfft/winsize); % the scale parameter M lby2 = (winsize - 1)/2; theta = -lby2:lby2; opwind = ones(winsize,1) * (theta .^2); % w(m,n)=m^2 opwind = opwind + opwind' + theta' * theta; % m^2 + n^2 + mn opwind = 1 - (2*mwind/nfft)^2 * opwind; % hex = ones(winsize,1) * theta; % m hex = abs(hex) + abs(hex') + abs(hex+hex'); hex = (hex < winsize); opwind = opwind .* hex; opwind = opwind * (4 * mwind^2) / (7 * pi^2) ; else opwind = 1; end elseif (min(m,n) == 1) % 1-D window passed: convert to 2-D window = window(:); if (any(imag(window) ~= 0)) disp(['1-D window has imaginary components: window ignored']) window = 1; end if (any(window < 0)) disp(['1-D window has negative components: window ignored']) window = 1; end lwind = length(window); windf = [window(lwind:-1:2); window]; % the full symmetric 1-D window = [window; zeros(lwind-1,1)]; opwind = (windf * windf') ... .* hankel(flipud(window), window); % w(m)w(n)w(m+n) winsize = length(window); else % 2-D window passed: use directly winsize = m; if (m ~= n) disp('2-D window is not square: window ignored') window = 1; winsize = m; end if (rem(m,2) == 0) disp('2-D window does not have odd length: window ignored') window = 1; winsize = m; end opwind = window; end % ---------------- accumulate triple products ---------------------- Bspec = zeros(nfft,nfft); mask = hankel([1:nfft],[nfft,1:nfft-1] ); % the hankel mask (faster) locseg = [1:nsamp]'; for krec = 1:nrecs xseg = y(locseg); Xf = fft(xseg-mean(xseg), nfft)/nsamp; CXf = conj(Xf); Bspec = Bspec + (Xf * Xf.') .* ... reshape(CXf(mask), nfft, nfft); locseg = locseg + nadvance; end Bspec = fftshift(Bspec)/(nrecs); % ----------------- frequency-domain smoothing ------------------------ if (winsize > 1) lby2 = (winsize-1)/2; Bspec = conv2(Bspec,opwind); Bspec = Bspec(lby2+1:lby2+nfft,lby2+1:lby2+nfft); end % ------------ 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(abs(Bspec),4,waxis,waxis),grid contour(waxis,waxis,abs(Bspec),4),grid on title('Bispectrum estimated via the direct (FFT) method') xlabel('f1'), ylabel('f2') set(gcf,'Name','Hosa BISPECD') return