请问双谱函数 bispecd 的参数设置问题

zx535777679

bispecd (y,  nfft, wind, segsamp, overlap)里面的各个参数怎么设置啊，原始信号是x1=A1*cos(2*pi*f1*t);%信号x2=A2*cos(2*pi*f2*t+C*cos(2*pi*fn*t));y=x1+x2

求大神指导

happy

应该是下面这个函数吧，函数前面的注释不是有很清楚的说明吗？

1. function [Bspec,waxis] = bispecd (y,  nfft, wind, nsamp, overlap)
   
2. %BISPECD Bispectrum estimation using the direct (fft-based) approach.
   
3. %[Bspec,waxis] = bispecd (y,  nfft, wind, segsamp, overlap)
   
4. %y    - data vector or time-series
   
5. %nfft - fft length [default = power of two > segsamp]
   
6. %wind - window specification for frequency-domain smoothing
   
7. %       if 'wind' is a scalar, it specifies the length of the side
   
8. %          of the square for the Rao-Gabr optimal window  [default=5]
   
9. %       if 'wind' is a vector, a 2D window will be calculated via
   
10. %          w2(i,j) = wind(i) * wind(j) * wind(i+j)
    
11. %       if 'wind' is a matrix, it specifies the 2-D filter directly
    
12. %segsamp - samples per segment [default: such that we have 8 segments]
    
13. %        - if y is a matrix, segsamp is set to the number of rows
    
14. %overlap - percentage overlap [default = 50]
    
15. %        - if y is a matrix, overlap is set to 0.
    
16. %
    
17. %Bspec   - estimated bispectrum: an nfft x nfft array, with origin
    
18. %          at the center, and axes pointing down and to the right.
    
19. %waxis   - vector of frequencies associated with the rows and columns
    
20. %          of Bspec;  sampling frequency is assumed to be 1.
    
21. 
    
22. %  Copyright (c) 1991-2001 by United Signals & Systems, Inc.
    
23. %       $Revision: 1.8 $
    
24. %  A. Swami   January 20, 1993.
    
25. 
    
26. %     RESTRICTED RIGHTS LEGEND
    
27. % Use, duplication, or disclosure by the Government is subject to
    
28. % restrictions as set forth in subparagraph (c) (1) (ii) of the
    
29. % Rights in Technical Data and Computer Software clause of DFARS
    
30. % 252.227-7013.
    
31. % Manufacturer: United Signals & Systems, Inc., P.O. Box 2374,
    
32. % Culver City, California 90231.
    
33. %
    
34. %  This material may be reproduced by or for the U.S. Government pursuant
    
35. %  to the copyright license under the clause at DFARS 252.227-7013.
    
36. 
    
37. % --------------------- parameter checks -----------------------------
    
38. 
    
39.     [ly, nrecs] = size(y);
    
40.     if (ly == 1) y = y(:);  ly = nrecs; nrecs = 1; end
    
41. 
    
42.     if (exist('nfft') ~= 1)            nfft = 128; end
    
43.     if (exist('overlap') ~= 1)      overlap = 50;  end
    
44.     overlap = min(99,max(overlap,0));
    
45.     if (nrecs > 1)                  overlap =  0;  end
    
46.     if (exist('nsamp') ~= 1)          nsamp = 0;   end
    
47.     if (nrecs > 1)                    nsamp = ly;  end
    
48. 
    
49.     if (nrecs == 1 & nsamp <= 0)
    
50.        nsamp = fix(ly/ (8 - 7 * overlap/100));
    
51.     end
    
52.     if (nfft  < nsamp)   nfft = 2^nextpow2(nsamp); end
    
53. 
    
54.     overlap  = fix(nsamp * overlap / 100);             % added 2/14
    
55.     nadvance = nsamp - overlap;
    
56.     nrecs    = fix ( (ly*nrecs - overlap) / nadvance);
    
57. 
    
58. 
    
59. % ------------------- create the 2-D window -------------------------
    
60.   if (exist('wind') ~= 1) wind = 5; end
    
61.   [m,n] = size(wind);
    
62.   window = wind;
    
63.   if (max(m,n) == 1)     % scalar: wind is size of Rao-Gabr window
    
64.      winsize = wind;
    
65.      if (winsize < 0) winsize = 5; end        % the window length L
    
66.      winsize = winsize - rem(winsize,2) + 1;  % make it odd
    
67.      if (winsize > 1)
    
68.         mwind   = fix (nfft/winsize);            % the scale parameter M
    
69.         lby2    = (winsize - 1)/2;
    
70. 
    
71.         theta  = -lby2:lby2;
    
72.         opwind = ones(winsize,1) * (theta .^2);       % w(m,n)=m^2
    
73.         opwind = opwind + opwind' + theta' * theta;   % m^2 + n^2 + mn
    
74.         opwind = 1 - (2*mwind/nfft)^2 * opwind;       %
    
75.         hex    = ones(winsize,1) * theta;             % m
    
76.         hex    = abs(hex) + abs(hex') + abs(hex+hex');
    
77.         hex    = (hex < winsize);
    
78.         opwind = opwind .* hex;
    
79.         opwind = opwind * (4 * mwind^2) / (7 * pi^2) ;
    
80.      else
    
81.         opwind = 1;
    
82.      end
    
83. 
    
84.   elseif (min(m,n) == 1)  % 1-D window passed: convert to 2-D
    
85.      window = window(:);
    
86.      if (any(imag(window) ~= 0))
    
87.         disp(['1-D window has imaginary components: window ignored'])
    
88.         window = 1;
    
89.      end
    
90.      if (any(window < 0))
    
91.         disp(['1-D window has negative components: window ignored'])
    
92.         window = 1;
    
93.      end
    
94.      lwind  = length(window);
    
95.      windf  = [window(lwind:-1:2); window];    % the full symmetric 1-D
    
96.      window = [window; zeros(lwind-1,1)];
    
97.      opwind = (windf * windf')      ...
    
98.               .* hankel(flipud(window), window); % w(m)w(n)w(m+n)
    
99.      winsize = length(window);
    
100. 
     
101.   else                    % 2-D window passed: use directly
     
102.     winsize = m;
     
103.     if (m ~= n)
     
104.        disp('2-D window is not square: window ignored')
     
105.        window = 1;
     
106.        winsize = m;
     
107.     end
     
108.     if (rem(m,2) == 0)
     
109.        disp('2-D window does not have odd length: window ignored')
     
110.        window = 1;
     
111.        winsize = m;
     
112.     end
     
113.     opwind  = window;
     
114.   end
     
115. 
     
116. % ---------------- accumulate triple products ----------------------
     
117. 
     
118.     Bspec    = zeros(nfft,nfft);
     
119. 
     
120.     mask = hankel([1:nfft],[nfft,1:nfft-1] );   % the hankel mask (faster)
     
121.     locseg = [1:nsamp]';
     
122.     for krec = 1:nrecs
     
123.         xseg   = y(locseg);
     
124.         Xf     = fft(xseg-mean(xseg), nfft)/nsamp;
     
125.         CXf    = conj(Xf);
     
126.         Bspec  = Bspec + (Xf * Xf.') .* ...
     
127.          reshape(CXf(mask), nfft, nfft);
     
128.         locseg = locseg + nadvance;
     
129.     end
     
130. 
     
131.     Bspec = fftshift(Bspec)/(nrecs);
     
132. 
     
133. 
     
134. 
     
135. % ----------------- frequency-domain smoothing ------------------------
     
136. 
     
137.   if (winsize > 1)
     
138.       lby2 = (winsize-1)/2;
     
139.       Bspec = conv2(Bspec,opwind);
     
140.       Bspec = Bspec(lby2+1:lby2+nfft,lby2+1:lby2+nfft);
     
141.   end
     
142. % ------------ contour plot of magnitude bispectum --------------------
     
143. 
     
144.    if (rem(nfft,2) == 0)
     
145.        waxis = [-nfft/2:(nfft/2-1)]'/nfft;
     
146.    else
     
147.        waxis = [-(nfft-1)/2:(nfft-1)/2]'/nfft;
     
148.    end
     
149. 
     
150.    hold off, clf
     
151. %  contour(abs(Bspec),4,waxis,waxis),grid
     
152.    contour(waxis,waxis,abs(Bspec),4),grid on
     
153.    title('Bispectrum estimated via the direct (FFT) method')
     
154.    xlabel('f1'), ylabel('f2')
     
155.    set(gcf,'Name','Hosa BISPECD')
     
156. return

复制代码