声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 1990|回复: 1

[编程技巧] 请问双谱函数 bispecd 的参数设置问题

[复制链接]
发表于 2015-7-29 22:09 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有账号?我要加入

x
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

求大神指导
回复
分享到:

使用道具 举报

发表于 2015-9-23 09:55 | 显示全部楼层
应该是下面这个函数吧,函数前面的注释不是有很清楚的说明吗?

  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. %  Copyright (c) 1991-2001 by United Signals & Systems, Inc.
  22. %       $Revision: 1.8 $
  23. %  A. Swami   January 20, 1993.

  24. %     RESTRICTED RIGHTS LEGEND
  25. % Use, duplication, or disclosure by the Government is subject to
  26. % restrictions as set forth in subparagraph (c) (1) (ii) of the
  27. % Rights in Technical Data and Computer Software clause of DFARS
  28. % 252.227-7013.
  29. % Manufacturer: United Signals & Systems, Inc., P.O. Box 2374,
  30. % Culver City, California 90231.
  31. %
  32. %  This material may be reproduced by or for the U.S. Government pursuant
  33. %  to the copyright license under the clause at DFARS 252.227-7013.

  34. % --------------------- parameter checks -----------------------------

  35.     [ly, nrecs] = size(y);
  36.     if (ly == 1) y = y(:);  ly = nrecs; nrecs = 1; end

  37.     if (exist('nfft') ~= 1)            nfft = 128; end
  38.     if (exist('overlap') ~= 1)      overlap = 50;  end
  39.     overlap = min(99,max(overlap,0));
  40.     if (nrecs > 1)                  overlap =  0;  end
  41.     if (exist('nsamp') ~= 1)          nsamp = 0;   end
  42.     if (nrecs > 1)                    nsamp = ly;  end

  43.     if (nrecs == 1 & nsamp <= 0)
  44.        nsamp = fix(ly/ (8 - 7 * overlap/100));
  45.     end
  46.     if (nfft  < nsamp)   nfft = 2^nextpow2(nsamp); end

  47.     overlap  = fix(nsamp * overlap / 100);             % added 2/14
  48.     nadvance = nsamp - overlap;
  49.     nrecs    = fix ( (ly*nrecs - overlap) / nadvance);


  50. % ------------------- create the 2-D window -------------------------
  51.   if (exist('wind') ~= 1) wind = 5; end
  52.   [m,n] = size(wind);
  53.   window = wind;
  54.   if (max(m,n) == 1)     % scalar: wind is size of Rao-Gabr window
  55.      winsize = wind;
  56.      if (winsize < 0) winsize = 5; end        % the window length L
  57.      winsize = winsize - rem(winsize,2) + 1;  % make it odd
  58.      if (winsize > 1)
  59.         mwind   = fix (nfft/winsize);            % the scale parameter M
  60.         lby2    = (winsize - 1)/2;

  61.         theta  = -lby2:lby2;
  62.         opwind = ones(winsize,1) * (theta .^2);       % w(m,n)=m^2
  63.         opwind = opwind + opwind' + theta' * theta;   % m^2 + n^2 + mn
  64.         opwind = 1 - (2*mwind/nfft)^2 * opwind;       %
  65.         hex    = ones(winsize,1) * theta;             % m
  66.         hex    = abs(hex) + abs(hex') + abs(hex+hex');
  67.         hex    = (hex < winsize);
  68.         opwind = opwind .* hex;
  69.         opwind = opwind * (4 * mwind^2) / (7 * pi^2) ;
  70.      else
  71.         opwind = 1;
  72.      end

  73.   elseif (min(m,n) == 1)  % 1-D window passed: convert to 2-D
  74.      window = window(:);
  75.      if (any(imag(window) ~= 0))
  76.         disp(['1-D window has imaginary components: window ignored'])
  77.         window = 1;
  78.      end
  79.      if (any(window < 0))
  80.         disp(['1-D window has negative components: window ignored'])
  81.         window = 1;
  82.      end
  83.      lwind  = length(window);
  84.      windf  = [window(lwind:-1:2); window];    % the full symmetric 1-D
  85.      window = [window; zeros(lwind-1,1)];
  86.      opwind = (windf * windf')      ...
  87.               .* hankel(flipud(window), window); % w(m)w(n)w(m+n)
  88.      winsize = length(window);

  89.   else                    % 2-D window passed: use directly
  90.     winsize = m;
  91.     if (m ~= n)
  92.        disp('2-D window is not square: window ignored')
  93.        window = 1;
  94.        winsize = m;
  95.     end
  96.     if (rem(m,2) == 0)
  97.        disp('2-D window does not have odd length: window ignored')
  98.        window = 1;
  99.        winsize = m;
  100.     end
  101.     opwind  = window;
  102.   end

  103. % ---------------- accumulate triple products ----------------------

  104.     Bspec    = zeros(nfft,nfft);

  105.     mask = hankel([1:nfft],[nfft,1:nfft-1] );   % the hankel mask (faster)
  106.     locseg = [1:nsamp]';
  107.     for krec = 1:nrecs
  108.         xseg   = y(locseg);
  109.         Xf     = fft(xseg-mean(xseg), nfft)/nsamp;
  110.         CXf    = conj(Xf);
  111.         Bspec  = Bspec + (Xf * Xf.') .* ...
  112.          reshape(CXf(mask), nfft, nfft);
  113.         locseg = locseg + nadvance;
  114.     end

  115.     Bspec = fftshift(Bspec)/(nrecs);



  116. % ----------------- frequency-domain smoothing ------------------------

  117.   if (winsize > 1)
  118.       lby2 = (winsize-1)/2;
  119.       Bspec = conv2(Bspec,opwind);
  120.       Bspec = Bspec(lby2+1:lby2+nfft,lby2+1:lby2+nfft);
  121.   end
  122. % ------------ contour plot of magnitude bispectum --------------------

  123.    if (rem(nfft,2) == 0)
  124.        waxis = [-nfft/2:(nfft/2-1)]'/nfft;
  125.    else
  126.        waxis = [-(nfft-1)/2:(nfft-1)/2]'/nfft;
  127.    end

  128.    hold off, clf
  129. %  contour(abs(Bspec),4,waxis,waxis),grid
  130.    contour(waxis,waxis,abs(Bspec),4),grid on
  131.    title('Bispectrum estimated via the direct (FFT) method')
  132.    xlabel('f1'), ylabel('f2')
  133.    set(gcf,'Name','Hosa BISPECD')
  134. return
复制代码
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 16:23 , Processed in 0.084678 second(s), 18 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

快速回复 返回顶部 返回列表