octave_packages/tsa-4.2.4/ar_spa.m

   1 function  [w,A,B,R,P,F,ip] = ar_spa(ARP,nhz,E);
   2 % AR_SPA decomposes an AR-spectrum into its compontents
   3 % [w,A,B,R,P,F,ip] = ar_spa(AR,fs,E);
   4 %
   5 %  INPUT:
   6 % AR   autoregressive parameters
   7 % fs    sampling rate, provide w and B in [Hz], if not given the result is in radians
   8 % E     noise level (mean square),  gives A and F in units of E, if not given as relative amplitude
   9 %
  10 %  OUTPUT
  11 % w     center frequency
  12 % A     Amplitude
  13 % B     bandwidth
  14 %       - less important output parameters -
  15 % R     residual
  16 % P     poles
  17 % ip    number of complex conjugate poles
  18 % real(F)       power, absolute values are obtained by multiplying with noise variance E(p+1)
  19 % imag(F)       assymetry, - " -
  20 %
  21 % All input and output parameters are organized in rows, one row
  22 % corresponds to the parameters of one channel
  23 %
  24 % see also ACOVF ACORF DURLEV IDURLEV PARCOR YUWA
  25 %
  26 % REFERENCES:
  27 % [1] Zetterberg L.H. (1969) Estimation of parameter for linear difference equation with application to EEG analysis. Math. Biosci., 5, 227-275.
  28 % [2] Isaksson A. and Wennberg, A. (1975) Visual evaluation and computer analysis of the EEG - A comparison. Electroenceph. clin. Neurophysiol., 38: 79-86.
  29 % [3] G. Florian and G. Pfurtscheller (1994) Autoregressive model based spectral analysis with application to EEG. IIG - Report Series, University of Technolgy Graz, Austria.
  30
  31 %       $Id: ar_spa.m 5090 2008-06-05 08:12:04Z schloegl $
  32 %       Copyright (C) 1996-2003 by Alois Schloegl <a.schloegl@ieee.org>
  33 %       This is part of the TSA-toolbox see also:
  34 %          http://hci.tugraz.at/schloegl/matlab/tsa/
  35 %          http://octave.sf.net/
  36 %
  37 %    This program is free software: you can redistribute it and/or modify
  38 %    it under the terms of the GNU General Public License as published by
  39 %    the Free Software Foundation, either version 3 of the License, or
  40 %    (at your option) any later version.
  41 %
  42 %    This program is distributed in the hope that it will be useful,
  43 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  44 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  45 %    GNU General Public License for more details.
  46 %
  47 %    You should have received a copy of the GNU General Public License
  48 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  49
  50
  51 [NTR,pp]=size(ARP);
  52
  53 R=zeros(size(ARP));
  54 P=zeros(size(ARP));
  55 w=zeros(size(ARP));
  56 A=zeros(size(ARP));
  57 B=zeros(size(ARP));
  58 F=zeros(size(ARP));
  59 F1 = F;
  60 for k = 1:NTR, %if ~mod(k,100),k, end;
  61         [r,p,tmp] = residue(1,[1 -ARP(k,:)]);
  62         [tmp,idx] = sort(-abs(r));
  63         R(k,:) = r(idx)';               % Residual,
  64         P(k,:) = p(idx)';               % Poles
  65         %r(k,:)=roots([1 -ARP(k,:)])';
  66         w(k,:) = angle(p(idx)');        % center frequency (in [radians])
  67         A(k,:) = 1./abs(polyval([1 -ARP(k,:)],exp(i*w(k,:))));  % Amplitude
  68         %A(k,:) = freqz(1,[1 -ARP(k,:)],w(k,:));        % Amplitude
  69         %A2(k,:) = abs(r)'./abs(exp(i*w(k,:))-r');   % Amplitude
  70         B(k,:) = -log(abs(p(idx)'));  % Bandwidth
  71
  72         if nargout < 6,
  73
  74         elseif 0,
  75                 F(k,:) = (1+sign(imag(r(idx)')))./(polyval([-ARP(k,pp-1:-1:1).*(1:pp-1) pp],1./p(idx).').*polyval([-ARP(k,pp:-1:1) 1],p(idx).'));
  76
  77         elseif 1;
  78                 a3 = polyval([-ARP(k,pp-1:-1:1).*(1:pp-1), pp],1./p(idx).');
  79                 a  = polyval([-ARP(k,pp:-1:1) 1],p(idx).');
  80                 %F(k,:) = (1+(imag(P(k,:))~=0))./(a.*a3);
  81                 F(k,:) = (1+sign(imag(P(k,:))))./(a.*a3);
  82         end;
  83 end;
  84
  85 A = A.*sqrt(E(:,ones(1,pp))/(2*pi*nhz));
  86 if nargin>1,
  87         if size(nhz,1)==1,
  88                 nhz = nhz(ones(NTR,1),:);
  89         end;
  90         w = w.*nhz(:,ones(1,pp))/(2*pi);
  91         B = B.*nhz(:,ones(1,pp))/(2*pi);
  92 end;
  93 if nargin>2,
  94         F = F.*E(:,ones(1,pp));
  95         F1 = F1.*E(:,ones(1,pp));
  96 end;
  97
  98 ip = sum(imag(P)~=0,2)/2;
  99 return;
 100
 101 np(:,1) = sum(imag(P')==0)';    % number of real poles
 102 np(:,2) = pp-np(:,1);           % number of imaginary poles
 103
 104