应该是下面这个函数吧,函数前面的注释不是有很清楚的说明吗?
- 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 <= 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
复制代码 |