001 function F=mctp2tc(freqPVR,utc,param,freqPVL)
002
003
004
005
006
007
008
009
010
011
012
013
014
015
016
017
018
019
020 if nargin<4
021 freqPVL=freqPVR;
022 end
023 if nargin<3
024 display('too few inputs parameters, stop')
025 break
026 end
027
028 u=levels(param);
029 udisc=fliplr(u);
030 ntc=sum(udisc>=utc);
031
032
033
034
035
036
037
038
039
040
041
042 if ntc < 2
043 display('the reference level out of range, stop')
044 break
045 end
046 if ntc > (param(3)-1)
047 display('the reference level out of range, stop')
048 break
049 end
050
051
052
053 nP=length(freqPVR);
054 for i=1:nP,
055 rowsum=sum(freqPVR(i,:));
056 if rowsum~=0
057 freqPVR(i,:)=freqPVR(i,:)/rowsum;
058 end
059 end
060 P=fliplr(freqPVR);
061
062 Ph=rot90(fliplr(freqPVL),-1);
063 for i=1:nP,
064 rowsum=sum(Ph(i,:));
065 if rowsum~=0
066 Ph(i,:)=Ph(i,:)/rowsum;
067 end
068 end
069 Ph=fliplr(Ph);
070
071
072 n=nP; F=zeros(n,n);
073
074 if ntc > n-1
075 display('index for mean-level out of range, stop')
076 break
077 end
078
079
080
081 F(1:ntc-1,1:(n-ntc))=freqPVL(1:ntc-1,1:(n-ntc));
082
083
084 F=cmat2nt(F);
085
086 for i=2:ntc,
087 for j=ntc:n-1,
088
089 if i<ntc
090
091 Ap=P(i:ntc-1,i+1:ntc); Bp=Ph(i+1:ntc,i:ntc-1);
092 dim_p=ntc-i;
093 tempp=zeros(dim_p,1);
094 I=eye(size(Ap));
095 if i==2
096 e=Ph(i+1:ntc,1);
097 else
098 e=sum(Ph(i+1:ntc,1:i-1)')';
099 end
100 if max(abs(e))>1e-10
101 if dim_p==1
102 tempp(1,1)=(Ap/(1-Bp*Ap)*e);
103 else
104 tempp=Ap*((I-Bp*Ap)\e);
105 end
106 end
107 end
108
109 if j>ntc
110
111 Am=P(ntc:j-1,ntc+1:j); Bm=Ph(ntc+1:j,ntc:j-1);
112 dim_m=j-ntc;
113 tempm=zeros(dim_m,1);
114 Im=eye(size(Am));
115 if j==n-1
116 em=P(ntc:j-1,n);
117 else
118 em=sum(P(ntc:j-1,j+1:n)')';
119 end
120 if max(abs(em))>1e-10
121 if dim_m==1
122 tempm(1,1)=(Bm/(1-Am*Bm)*em);
123 else
124 tempm=Bm*((Im-Am*Bm)\em);
125 end
126 end
127 end
128
129 if (j>ntc)*(i<ntc)
130 F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-ntc:-1:n-j+1)*tempm;
131 F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1);
132 F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm;
133 end
134 if (j==ntc)*(i<ntc)
135 F(i,n-j+1)=F(i,n-j+1)+tempp'*freqPVL(i:ntc-1,n-j:-1:1)*ones(n-j,1);
136 for k=1:ntc
137 F(i,n-k+1)=F(i,n-ntc+1);
138 end
139 end
140 if (j>ntc)*(i==ntc)
141 F(i,n-j+1)=F(i,n-j+1)+ones(1,i-1)*freqPVL(1:i-1,n-ntc:-1:n-j+1)*tempm;
142 for k=ntc:n
143 F(k,n-j+1)=F(ntc,n-j+1);
144 end
145 end
146
147 end
148 end
149 F;
150 F=nt2cmat(F);
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171