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);
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
14 % - less important output parameters -
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, - " -
21 % All input and output parameters are organized in rows, one row
22 % corresponds to the parameters of one channel
24 % see also ACOVF ACORF DURLEV IDURLEV PARCOR YUWA
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.
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/
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.
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.
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/>.
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
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).'));
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);
85 A = A.*sqrt(E(:,ones(1,pp))/(2*pi*nhz));
88 nhz = nhz(ones(NTR,1),:);
90 w = w.*nhz(:,ones(1,pp))/(2*pi);
91 B = B.*nhz(:,ones(1,pp))/(2*pi);
94 F = F.*E(:,ones(1,pp));
95 F1 = F1.*E(:,ones(1,pp));
98 ip = sum(imag(P)~=0,2)/2;
101 np(:,1) = sum(imag(P')==0)'; % number of real poles
102 np(:,2) = pp-np(:,1); % number of imaginary poles