home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 8.ddi / SIGNAL.DI$ / SPECTRUM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.5 KB  |  152 lines
  1. function P = spectrum(x,y,m,noverlap)
  2. %SPECTRUM Power spectrum estimate of one or two data sequences.
  3. %    P = SPECTRUM(X,Y,M) performs FFT analysis of the two sequences
  4. %    X and Y using the Welch method of power spectrum estimation. 
  5. %    The X and Y sequences of N points are divided into K sections of
  6. %    M points each (M must be a power of two). Using an M-point FFT,
  7. %    successive sections are Hanning windowed, FFT'd and accumulated.
  8. %    SPECTRUM returns the M/2 by 8 array
  9. %       P = [Pxx Pyy Pxy Txy Cxy Pxxc Pyyc Pxyc]
  10. %    where
  11. %       Pxx  = X-vector power spectral density
  12. %       Pyy  = Y-vector power spectral density
  13. %       Pxy  = Cross spectral density
  14. %       Txy  = Complex transfer function from X to Y
  15. %             (Use ABS and ANGLE for magnitude and phase)
  16. %       Cxy  = Coherence function between X and Y
  17. %       Pxxc,Pyyc,Pxyc = Confidence range (95 percent).
  18. %
  19. %    See SPECPLOT to plot these results.
  20. %    P = SPECTRUM(X,Y,M,NOVERLAP) specifies that the M-point sections
  21. %    should overlap NOVERLAP points.
  22. %    Pxx = SPECTRUM(X,M) and SPECTRUM(X,M,NOVERLAP) return the single
  23. %    sequence power spectrum and confidence range.
  24. %
  25. %       See also ETFE, SPA, and ARX in the Identification Toolbox.
  26.  
  27. %    J.N. Little 7-9-86
  28. %    Revised 4-25-88 CRD, 12-20-88 LS, 8-31-89 JNL, 8-11-92 LS
  29. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  30.  
  31. %    The units on the power spectra Pxx and Pyy are such that, using
  32. %    Parseval's theorem: 
  33. %
  34. %         SUM(Pxx)/LENGTH(Pxx) = SUM(X.^2)/LENGTH(X) = COV(X)
  35. %
  36. %    The RMS value of the signal is the square root of this.
  37. %    If the input signal is in Volts as a function of time, then
  38. %    the units on Pxx are Volts^2*seconds = Volt^2/Hz.
  39. %    To normalize Pxx so that a unit sine wave corresponds to
  40. %    one unit of Pxx, use Pn = 2*SQRT(Pxx/LENGTH(Pxx))
  41. %
  42. %    Here are the covariance, RMS, and spectral amplitude values of
  43. %    some common functions:
  44. %         Function   Cov=SUM(Pxx)/LENGTH(Pxx)   RMS        Pxx
  45. %         a*sin(w*t)        a^2/2            a/sqrt(2)   a^2*LENGTH(Pxx)/4
  46. %Normal:  a*rand(t)         a^2              a           a^2
  47. %Uniform: a*rand(t)         a^2/12           a/sqrt(12)  a^2/12
  48. %    
  49. %    For example, a pure sine wave with amplitude A has an RMS value
  50. %    of A/sqrt(2), so A = SQRT(2*SUM(Pxx)/LENGTH(Pxx)).
  51. %
  52. %    See Page 556, A.V. Oppenheim and R.W. Schafer, Digital Signal
  53. %    Processing, Prentice-Hall, 1975.
  54.  
  55. if (nargin == 2), m = y; noverlap = 0; end
  56. if (nargin == 3)
  57.     if (max(size(y)) == 1)
  58.         noverlap = m;
  59.         m = y;
  60.         nargin = 2;
  61.     else
  62.         noverlap = 0;
  63.     end
  64. end
  65.  
  66. x = x(:);        % Make sure x and y are column vectors
  67. y = y(:);
  68. n = max(size(x));    % Number of data points
  69. k = fix((n-noverlap)/(m-noverlap));    % Number of windows
  70.                     % (k = fix(n/m) for noverlap=0)
  71. index = 1:m;
  72. w = hanning(m);        % Window specification; change this if you want:
  73.             % (Try HAMMING, BLACKMAN, BARTLETT, or your own)
  74. KMU = k*norm(w)^2;    % Normalizing scale factor
  75.  
  76. if (nargin == 2)    % Single sequence case.
  77.     Pxx = zeros(m,1); Pxx2 = zeros(m,1);
  78.     for i=1:k
  79.         xw = w.*detrend(x(index));
  80.         index = index + (m - noverlap);
  81.         Xx = abs(fft(xw)).^2;
  82.         Pxx = Pxx + Xx;
  83.         Pxx2 = Pxx2 + abs(Xx).^2;
  84.     end
  85.     % Select first half
  86.     select = [1:m/2];
  87.     Pxx = Pxx(select);
  88.     Pxx2 = Pxx2(select);
  89.     cPxx = zeros(m/2,1);
  90.     if k > 1
  91.         c = (k.*Pxx2-abs(Pxx).^2)./(k-1);
  92.         c = max(c,zeros(m/2,1));
  93.         cPxx = sqrt(c);
  94.     end
  95.     pp = 0.95; % 95 percent confidence.
  96.     f = sqrt(2)*erfinv(pp);  % Equal-tails.
  97.     P = [Pxx f.*cPxx]/KMU;
  98.     return
  99. end
  100.  
  101. Pxx = zeros(m,1); % Dual sequence case.
  102. Pyy = Pxx; Pxy = Pxx; Pxx2 = Pxx; Pyy2 = Pxx; Pxy2 = Pxx;
  103.  
  104. for i=1:k
  105.     xw = w.*detrend(x(index));
  106.     yw = w.*detrend(y(index));
  107.     index = index + (m - noverlap);
  108.     Xx = fft(xw);
  109.     Yy = fft(yw);
  110.     Yy2 = abs(Yy).^2;
  111.     Xx2 = abs(Xx).^2;
  112.     Xy  = Yy .* conj(Xx);
  113.     Pxx = Pxx + Xx2;
  114.     Pyy = Pyy + Yy2;
  115.     Pxy = Pxy + Xy;
  116.     Pxx2 = Pxx2 + abs(Xx2).^2;
  117.     Pyy2 = Pyy2 + abs(Yy2).^2;
  118.     Pxy2 = Pxy2 + Xy .* conj(Xy);
  119. end
  120.  
  121. % Select first half
  122. select = [1:m/2];
  123.  
  124. Pxx = Pxx(select);
  125. Pyy = Pyy(select);
  126. Pxy = Pxy(select);
  127. Pxx2 = Pxx2(select);
  128. Pyy2 = Pyy2(select);
  129. Pxy2 = Pxy2(select);
  130.  
  131. cPxx = zeros(m/2,1);
  132. cPyy = cPxx;
  133. cPxy = cPxx;
  134. if k > 1
  135.    c = max((k.*Pxx2-abs(Pxx).^2)./(k-1),zeros(m/2,1));
  136.    cPxx = sqrt(c);
  137.    c = max((k.*Pyy2-abs(Pyy).^2)./(k-1),zeros(m/2,1));
  138.    cPyy = sqrt(c);
  139.    c = max((k.*Pxy2-abs(Pxy).^2)./(k-1),zeros(m/2,1));
  140.    cPxy = sqrt(c);
  141. end
  142.  
  143. Txy = Pxy./Pxx;
  144. Cxy = (abs(Pxy).^2)./(Pxx.*Pyy);
  145.  
  146. pp = 0.95; % 95 percent confidence.
  147. f = sqrt(2)*erfinv(pp);  % Equal-tails.
  148.  
  149. P = [ [Pxx Pyy Pxy]./KMU ...
  150.       Txy Cxy ...
  151.       f.*[cPxx cPyy cPxy]./KMU ];
  152.