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-2001 by United Signals & Systems, Inc.% $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 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
评论
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
查看更多评论
添加红包
