plota.m

matlab数字信号处理工具箱

                                else
                                        area(X.f,abs(tmp));
                                        axis([0,max(X.f),0,1])
                                end;
                                set(s,'FontSize',6); 
                                %end
                        end;
                end;
                suptitle('PDC') 
        elseif (length(size(X.PDC))==4) & (size(X.PDC,4)==2),
                M = size(X.PDC,1);
                for k1 = 1:M;
                        for k2 = 1:M;
                                subplot(M,M,k2+(k1-1)*M);
                                tmp = abs(squeeze(X.PDC(k1,k2,:,1)));
                                plot(X.f,squeeze(X.PDC(k1,k2,:,1)),'b',X.f,squeeze(X.PDC(k1,k2,:,2)),'g');        
                                axis([0,max(X.f),0,1])        
                                if k1==1,
                                        ylabel(X.ElectrodeName{k2});
                                end;
                                if 1==k2,
                                        title(X.ElectrodeName{k1});
                                end;
                        end;
                end;
                suptitle('PDC')
        elseif length(size(X.PDC))==4,
                if nargin<2,
                        list1=1:size(X.PDC,1);
                        list2=1:size(X.PDC,2);
                        list3=1:size(X.PDC,3);
                else
                        list1=arg2;
                        list2=arg3;
                        list3=arg4';
                end;
                
                M = size(X.PDC,1);
                for k1 = 1:length(list1);
                        for k2 = 1:length(list2);
                                s = subplot(length(list2),length(list1),k2+(k1-1)*length(list1));
                                tmp = abs(squeeze(X.PDC(list1(k1),list2(k2),list3,:)));
                                [a,b] = meshgrid(X.t,list3);
                                mesh(a,b,abs(tmp));
                                v = axis; v(5:6)=[0,1]; axis(v);
                                xlabel('Time','FontSize',6);
                                ylabel('Freq','FontSize',6);
                                set(s,'FontSize',6);
                        end;
                end;
                suptitle('Time-varying PDC')
                figure;
                for k1 = 1:length(list1);
                        for k2 = 1:length(list2);
                                s = subplot(length(list2),length(list1),k2+(k1-1)*length(list1));
                                tmp = abs(squeeze(X.PDC(list1(k1),list2(k2),list3,:)));
                                imagesc(X.t,list3,tmp,[0 1]);
                                axis([min(X.t) max(X.t) min(list3) max(list3)])
                                xlabel('Time','FontSize',5);
                                ylabel('Freq','FontSize',5);
                                set(s,'YDir','normal','FontSize',5);
                                c = colorbar;
                                set(c,'FontSize',5);
                        end
                end
                suptitle('Time-varying PDC')
                
        elseif length(size(X.PDC))==5,
                [M1,M2,M3,M4,M5] = size(X.PDC);
                load cmap;
                for m1=1:M1
                        for m2=1:M2
                                if m1 ~= m2
                                        colormap(cmap);
                                        s = subplot(M1,M1,m1+(m2-1)*M1);
                                        imagesc(squeeze(X.DS(m1,m2,:,:)),[-1,1]);
                                        xlabel('Time','FontSize',5);
                                        ylabel('Freq','FontSize',5);
                                        set(s,'YDir','normal','FontSize',5);
                                        c = colorbar;
                                        set(c,'FontSize',5);
                                end
                        end
                end   
                suptitle('Significane difference of time-varying PDC')
        end; 
        
elseif strcmp(X.datatype,'MVAR'),
        if ~isfield(X,'A') | ~isfield(X,'B'),
                fprintf(2,'Error PLOTA: MVAR missing input data\n');
                return;
        end;
        
        [K1,K2] = size(X.A);
        p = K2/K1-1;
        %a=ones(1,p+1);
        [K1,K2] = size(X.B);
        q = K2/K1-1;
        %b=ones(1,q+1);
        if ~isfield(X,'C');
                X.C=ones(K1,K1);
        end;
        if nargin<2,
                Mode= 'DTF';
                Fs  = 1;
                f   = (0:128)/128*pi;
        else
                Mode = arg2;
        end;
        if nargin<3,
                N=512;
        else
                N=arg3;
        end;
        if nargin<4
                Fs=pi*2;        
        else
                Fs=arg4;
        end;
        if all(size(N)==1)	
                f = (1:N)/N/2*Fs;
        else
                f = N;
                N = length(N);
        end;
        if isfield(X,'SampleRate'); 
                Fs=X.SampleRate;
        end;
        if isfield(X,'ElectrodeName'); 
                if isstruct(X.ElectrodeName);
                        Label = X.ElectrodeName;
                else
                        for k=1:K1,
                                Label{k}=X.ElectrodeName(k,:);
                        end;
                end;
        else
                for k=1:K1,
                        Label{k}=sprintf('#%02i',k);
                end;
        end;
        
        s = exp(i*2*pi*f/Fs);
        z = i*2*pi/Fs;
        
        h=zeros(K1,K1,N);
        SP=zeros(K1,K1,N);
        DTF=zeros(K1,K1,N);
        COH=zeros(K1,K1,N);
        COH2=zeros(K1,K1,N);
        PDC=zeros(K1,K1,N);
        PDCF=zeros(K1,K1,N);
        invC=inv(X.C);
        tmp1=zeros(1,K1);
        tmp2=zeros(1,K1);
        for n=1:N,
                atmp = zeros(K1,K1);
                for k = 1:p+1,
                        atmp = atmp + X.A(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
                end;        
                btmp = zeros(K1,K2);
                for k = 1:q+1,
                        btmp = btmp + X.B(:,k*K1+(1-K1:0))*exp(z*(k-1)*f(n));
                end;        
                h(:,:,n) = atmp\btmp;        
                S(:,:,n) = h(:,:,n)*X.C*h(:,:,n)';        
                
                for k1 = 1:K1,
                        tmp = squeeze(atmp(:,k1));
                        tmp1(k1) = sqrt(tmp'*tmp);
                        tmp2(k1) = sqrt(tmp'*invC*tmp);
                end;
                
                %PDCF(:,:,n,kk) = abs(atmp)./tmp2(ones(1,K1),:);
                %PDC(:,:,n,kk)  = abs(atmp)./tmp1(ones(1,K1),:);
                PDCF(:,:,n) = abs(atmp)./tmp2(ones(1,K1),:);
                PDC(:,:,n)  = abs(atmp)./tmp1(ones(1,K1),:);
        end;
        
        if strcmpi(Mode,'Spectrum'),      
                maxS=max(abs(S(:)));
                minS=min(abs(S(:)));
                for k1=1:K1;
                        for k2=1:K2;
                                subplot(K1,K2,k2+(k1-1)*K1);
                                semilogy(f,squeeze(abs(S(k1,k2,:))));        
                                axis([0,max(f),minS,maxS]);
                                if k2==1;
                                        ylabel(Label{k1});
                                end;
                                if k1==1;
                                        title(Label{k2});
                                end;
                        end;
                end;
                suptitle('Power spectrum')
                return;
        elseif strcmpi(Mode,'Phase'),      
                maxS=max(angle(S(:)));
                minS=min(angle(S(:)));
                figure(2); clf;
                for k1=1:K1;
                        for k2=1:K2;
                                subplot(K1,K2,k2+(k1-1)*K1);
                                plot(f,unwrap(squeeze(angle(S(k1,k2,:))))*180/pi);        
                                axis([0,max(f),-360,360])        
                                if k2==1;
                                        ylabel(Label{k1});
                                end;
                                if k1==1;
                                        title(Label{k2});
                                end;
                        end;
                end;
                suptitle('phase')
                return;
        elseif strcmpi(Mode,'PDC'),      
                for k1=1:K1;
                        for k2=1:K2;
                                subplot(K1,K2,k2+(k1-1)*K1);
                                area(f,squeeze(PDC(k1,k2,:)));        
                                axis([0,max(f),0,1]);
                                if k2==1;
                                        ylabel(Label{k1});
                                end;
                                if k1==1;
                                        title(Label{k2});
                                end;
                        end;
                end;
                suptitle('partial directed coherence PDC');
                return;
        end;        
        
        DC = zeros(K1);
        for k = 1:p,
                DC = DC + X.A(:,k*K1+(1:K1)).^2;
        end;
        if strcmpi(Mode,'DC'),      
                fprintf(2,'Warning PLOTA: DC not implemented yet\n');
                return;
        end;        
        
        %%%%% directed transfer function
        for k1=1:K1;
                DEN=sqrt(sum(abs(h(k1,:,:)).^2,2));	        
                for k2=1:K2;
                        %COH2(k1,k2,:) = abs(S(k1,k2,:).^2)./(abs(S(k1,k1,:).*S(k2,k2,:)));
                        COH(k1,k2,:) = abs(S(k1,k2,:))./sqrt(abs(S(k1,k1,:).*S(k2,k2,:)));
                        %DTF(k1,k2,:) = sqrt(abs(h(k1,k2,:).^2))./DEN;	        
                        DTF(k1,k2,:) = abs(h(k1,k2,:))./DEN;
                end;
        end;
        
        if strcmpi(Mode,'Coherence'),      
                for k1=1:K1;
                        for k2=1:K2;
                                subplot(K1,K2,k2+(k1-1)*K1);
                                plot(f,squeeze(COH(k1,k2,:)));        
                                axis([0,max(f),0,1])
                                if k2==1;
                                        ylabel(Label{k1});
                                end;
                                if k1==1;
                                        title(Label{k2});
                                end;
                        end;
                end;
                suptitle('Ordinary coherence')
                return;
        elseif strcmpi(Mode,'DTF'),      
                for k1=1:K1;
                        for k2=1:K2;
                                subplot(K1,K2,k2+(k1-1)*K1);
                                area(f,squeeze(DTF(k1,k2,:)));        
                                axis([0,max(f),0,1]);
                                if k2==1;
                                        ylabel(Label{k1});
                                end;
                                if k1==1;
                                        title(Label{k2});
                                end;
                        end;
                end;
                suptitle('directed transfer function DTF');        
                return;
        end;        
        
elseif strcmp(X.datatype,'TF-MVAR') & (nargin>1) %& ~any(strmatch(arg2,{'S1','logS1',})),
        
        %GF = {'C','DC','AR','PDC','DTF','dDTF','ffDTF','COH','pCOH','pCOH2','S','h','phaseS','phaseh','coh','logh','logS'};
        
        if nargin<2,
                arg2 = 'S1';