home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 2.ddi / MUTOOLS1.DI$ / VSPECT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  7.6 KB  |  265 lines
  1. % function P = vspect(x,y,m,noverlap,window)
  2. %
  3. %   Spectral analysis routine modeled after SPECTRUM in the
  4. %   Signal Processing Toolbox.  X (and optionally Y) are VARYING
  5. %   matrices representing time domain sequences.  M point ffts are
  6. %   used with NOVERLAP points of overlap.  WINDOW is an optional
  7. %   string specifying a window (default is 'hanning'). Currently
  8. %   supported options within the Signal Processing toolbox are:
  9. %   'hanning', 'hamming', 'boxcar','blackman','bartlett', and 'triang'.
  10. %
  11. %   The returned VARYING matrix, P, has the following columns
  12. %   [Pxx Pyy Pxy Txy Cxy].  Pxx (Pyy) is the power spectrum of X (Y),
  13. %   Pxy is the cross spectrum, Txy is the transfer function, and
  14. %   Cxy is the coherence.  For a single input sequence only Pxx is
  15. %   returned.  The returned independent variable is frequency in
  16. %   radians/seconds. The signal Y (or X in the single signal case) 
  17. %   can be a vector of signals.  The transfer function is calculated 
  18. %   for each output. This corresponds to single-input/multi-output
  19. %   transfer function estimation.  Each row of P corresponds to the 
  20. %   appropriate row of Y.
  21. %
  22. %   See also: FFT, IFFT, SPECTRUM, VFFT and VIFFT.
  23.  
  24. function P = vspect(x,y,m,noverlap,window)
  25.  
  26. if nargin == 0,
  27.     disp('usage: P = vspect(x,y,m,noverlap,window)')
  28.     disp('       P = vspect(x,m,noverlap,window)')
  29.     return
  30.     end
  31.  
  32. if nargin < 2 | nargin > 5
  33.     disp('incorrect number of input arguments')
  34.     disp('usage: P = vspect(x,y,m,noverlap,window)')
  35.     return
  36.     end
  37.  
  38. [xtype,xnr,xnc,xnpts] = minfo(x);
  39. if xtype ~= 'vary'
  40.     error('input signal must be VARYING matrix')
  41.     return
  42.     end
  43.  
  44. [ytype,ynr,ync,ynpts] = minfo(y);
  45.  
  46. %    Use the second argument to determine whether or not
  47. %    we are dealing with a single spectrum case.
  48.  
  49. if ytype == 'cons',            % single spectrum case.
  50.     if nargin > 4,
  51.         error('too many arguments for single input spectrum')
  52.         return
  53.         end
  54.     if nargin < 4,
  55.         window = 'hanning';
  56.     else
  57.         window = noverlap;
  58.         end
  59.     if nargin < 3,
  60.         noverlap = 0;
  61.     else
  62.         noverlap = m;
  63.         end
  64.     m = y;
  65.     if fix(m/2) ~= m/2,
  66.         error('fft length must be even')
  67.         return
  68.         end
  69.  
  70. %       now create the window of the appropriate length.
  71.  
  72.     if ~isstr(window),
  73.         error('window argument must be a string')
  74.         return
  75.     else
  76.         eval(['w = ' window '(m);']);
  77.         end
  78.  
  79. %      calculate the number of windows and the appropriate
  80. %      normalization factors.
  81.  
  82.     k = fix((xnpts - noverlap)/(m - noverlap));
  83.     disp([ int2str(k) ' ' window ' windows in averaging calculation']);
  84.     KMU = k*norm(w)^2;
  85.  
  86. %       calculate a frequency vector for creating the output.
  87.  
  88.     t = getiv(x);
  89.     tinc = t(2) - t(1);
  90.     finc = 2*pi/(tinc*m);
  91.     omega = [0:finc:finc*(m/2-1)];
  92.  
  93. %      go through the averaging procedure for each row
  94.  
  95.     select = [2 2:m/2];
  96.     Pxx = [];
  97.     for i = 1:xnr,
  98.         Pxxcol = [];
  99.         for j = 1:xnc,
  100.             sisoPxx = zeros(m,1);
  101.             index = [1:m];
  102.             xdat = vunpck(sel(x,i,j));
  103.             for l = 1:k,
  104.                 xw = w.*detrend(xdat(index));
  105.                 index = index + (m - noverlap);
  106.                 Xx = abs(fft(xw)).^2;
  107.                 sisoPxx = sisoPxx + Xx;
  108.                 end
  109.             sisoPxx = vpck((1/KMU)*sisoPxx(select),omega);
  110.             Pxxcol = sbs(Pxxcol,sisoPxx);
  111.             end
  112.         Pxx = abv(Pxx,Pxxcol);
  113.         end
  114.     P = Pxx;
  115.  
  116. elseif ytype == 'vary',         % cross spectrum case
  117.     if nargin < 5,
  118.         window = 'hanning';
  119.         end
  120.     if nargin < 4,
  121.         noverlap = 0;
  122.         end
  123.     if nargin < 3,
  124.         error('FFT length must be specified')
  125.         return
  126.         end
  127.     if fix(m/2) ~= m/2,
  128.         error('fft length must be even')
  129.         return
  130.         end
  131.  
  132. %       check that the x and y vectors are over the same
  133. %       independent variables.
  134.  
  135.     if indvcmp(x,y) > 1,
  136.         error('INDEPENDENT VARIABLEs of the signals do not match')
  137.         return
  138.         end
  139.  
  140. %       check that y is a row vector.  If necessary perform its
  141. %       transpose.  This is done so that selecting different
  142. %       elements is faster.
  143.  
  144.     [ytype,ynr,ync,ynpts] = minfo(y);
  145.     if ynr > 1 & ync > 1,
  146.         error('output signal (y) is a matrix')
  147.         return
  148.     elseif ynr > 1,
  149.         y = transp(y);
  150.         [ytype,ynr,ync,ynpts] = minfo(y);
  151.         end
  152.  
  153. %       make sure that x is a 1 x 1 VARYING matrix.  This means
  154. %       that we can only do SISO or SIMO spectrum cases.
  155.  
  156.     if xnr~= 1 | xnc ~= 1,
  157.         error('the input signal (x) cannot be a vector')
  158.         return
  159.         end
  160.  
  161. %       now create the window of the appropriate length.
  162.  
  163.     if ~isstr(window),
  164.         error('window argument must be a string')
  165.         return
  166.     else
  167.         eval(['w = ' window '(m);']);
  168.         end
  169.  
  170. %      calculate the number of windows and the appropriate
  171. %      normalization factors.
  172.  
  173.     k = fix((xnpts - noverlap)/(m - noverlap));
  174.     disp([ int2str(k) ' ' window ' windows in averaging calculation']);
  175.     KMU = k*norm(w)^2;
  176.  
  177.  
  178.     [xdat,xptr,t] = vunpck(x);
  179.     [ydat,yptr,t] = vunpck(y);
  180.  
  181. %      go through the averaging procedure x.  The results of the fft
  182. %      are stored side by side for each k in Xxstack.  This is so
  183. %      that they can be used to calculate Pxy when going through
  184. %      the averaging procedure for y.
  185.  
  186.     Pxx = zeros(m,1);
  187.     Xxstack = [];
  188.     index = [1:m];
  189.     for i = 1:k,
  190.         xw = w.*detrend(xdat(index));
  191.         index = index + (m - noverlap);
  192.         Xx = fft(xw);
  193.         Xxstack = [Xxstack Xx ];
  194.         Pxx = Pxx + abs(Xx).^2;
  195.         end
  196.     select = [2 2:m/2];
  197.     Pxx = Pxx(select);
  198.  
  199.     Pyy = [];
  200.     Pxy = [];
  201.     for j = 1:ync,
  202.         Pyycol = zeros(m,1);
  203.         Pxycol = zeros(m,1);
  204.         index = [1:m];
  205.         for i = 1:k,
  206.             yw = w.*detrend(ydat(index,j));
  207.             index = index + (m - noverlap);
  208.             Yy = fft(yw);
  209.             Xy = Yy .* conj(Xxstack(:,i));
  210.             Pyycol = Pyycol + abs(Yy).^2;
  211.             Pxycol = Pxycol + Xy;
  212.             end
  213.         Pyycol = Pyycol(select);
  214.         Pyy = [Pyy Pyycol];
  215.         Pxycol = Pxycol(select);
  216.         Pxy = [Pxy Pxycol];
  217.         end
  218.  
  219. %        calculate the transfer function and coherence
  220.  
  221.     Txy = [];
  222.     Cxy = [];
  223.     for j = 1:ync,
  224.         Txy = [Txy Pxy(:,j)./Pxx];
  225.         Cxy = [Cxy (abs(Pxy(:,j).^2)./(Pxx.*Pyy(:,j)))];
  226.         end
  227.  
  228. %      form a frequency vector for appending to the
  229. %      result.
  230.  
  231.     tinc = t(2) - t(1);
  232.     finc = 2*pi/(tinc*m);
  233.     omega = [0:finc:finc*(m/2-1)];
  234.  
  235. %      Put together the final spectrum and do a transpose to
  236. %      get it back in column vector form.  Pxx is padded with
  237. %      zeros to get everything dimensionally correct.
  238.  
  239.     if ync > 1,
  240.         Pxx = Pxx * [1 zeros(1,ync-1)];
  241.         end
  242.     vPxx = vpck((1/KMU)*Pxx,omega);
  243.     vPxx = transp(vPxx);
  244.     vPyy = vpck((1/KMU)*Pyy,omega);
  245.     vPyy = transp(vPyy);
  246.     vPxy = vpck((1/KMU)*Pxy,omega);
  247.     vPxy = transp(vPxy);
  248.     vTxy = vpck(Txy,omega);
  249.     vTxy = transp(vTxy);
  250.     vCxy = vpck(Cxy,omega);
  251.     vCxy = transp(vCxy);
  252.  
  253. %      stick all the bits together so each row of the output
  254. %      has
  255. %                [ Pxx Pyy Pxy Txy Cxy]
  256.  
  257. P = sbs(vPxx, sbs(vPyy, sbs(vPxy, sbs(vTxy,vCxy))));
  258.  
  259. else                            % error -  it must be a system matrix
  260.     error('SYSTEM matrix arguments are not allowed')
  261.     return
  262.     end
  263. %
  264. % Copyright MUSYN INC 1991,  All Rights Reserved
  265.