0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 function [H,k1,k2]=constantes_sin_continua(fs,tam2,prediccion,fun,factor1,factor2);
0015
0016 fN=fs/2;
0017
0018 tam=32;
0019
0020 lcuan=32;
0021
0022
0023
0024 c=linspace(0,1,100);
0025
0026 [fx,fy]=freqspace(tam);
0027
0028 [IAFet,CSFet]=iafet(abs(fx*fN)',1,0,c,[1 1 0],1);
0029
0030
0031
0032
0033
0034
0035
0036
0037 R=cumsum(IAFet')';
0038
0039
0040 mm = maxi(R);
0041
0042 R=R/maxi(R);
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055 [fx,fy]=freqspace(tam,'meshgrid');
0056 Fx=fx*fN;
0057 Fy=fy*fN;
0058
0059 F=sqrt(Fx.^2+Fy.^2);
0060
0061
0062
0063 fx=zigzag(Fx);
0064 fy=zigzag(Fy);
0065 f=zigzag(F);
0066
0067
0068
0069 ind=find(f==0);
0070
0071
0072
0073
0074
0075
0076 H=zeros(lcuan^2,lcuan^2);
0077
0078 cont = 1;
0079 for i=1:lcuan
0080 for j=1:lcuan
0081 s = feval(fun,f(cont),factor1,factor2);
0082 hh=exp(-((Fx-fx(cont)).^2+(Fy-fy(cont)).^2)/s.^2);
0083 H(cont,:)=(zigzag(hh))';
0084 cont = cont + 1;
0085 end
0086 end
0087
0088
0089 H = H./ repmat(sum(H,2),1,lcuan^2);
0090 H(ind,ind)=1;
0091
0092
0093
0094
0095 h11f=dezigzag(diag(H));
0096
0097
0098 h11f=h11f(tam/2+1,:);
0099
0100
0101
0102
0103 rhs = interp1([1:13],h11f(17:29),[14 15 16],'linear','extrap');
0104
0105 rhs =[h11f(17:29) rhs];
0106 r=reverse(rhs);
0107
0108 h11f=[rhs(end) r(1:end-1) rhs];
0109
0110
0111
0112
0113
0114 [csf,logciaf]=ngan(fN,tam/2);
0115 r=reverse(csf);
0116 k1=[csf(end) r(1:end-1) csf];
0117
0118
0119 kk2=[ 0.1800 0.1800 0.1550 0.1100 0.0800 0.0500 0.0306 0.0219 0.0146 0.0102 0.0078 0.0068 0.0067 0.0070 0.0067 0.0083 0.3292 0.0083 0.0067 0.0070 0.0067 0.0068 0.0078 0.0102 0.0146 0.0219 0.0306 0.0500 0.0800 0.1100 0.1550 0.1800];
0120
0121
0122 opciones=optimset('Diagnostics','off','Display','off','TolX',1e-15,'TolFun',1e-15,'MaxIter',10000,'MaxFunEvals',10000);
0123
0124
0125
0126 for i=1:tam
0127
0128 x0 = fminsearch('error_naka9_1',k1(i),opciones,kk2(i),h11f(i),c,R(i,:));
0129
0130 kk1(i)=x0;
0131 end
0132
0133
0134
0135 for i=1:32
0136 RR(i,:)=kk1(i)/10+ kk1(i)*(c.^2)./(kk2(i)+h11f(i)*c.^2);
0137 end
0138
0139 x0 = fminsearch('error_naka9_1_ajustefino',[1 0.3],opciones,kk1,kk2,h11f,c(1:30),R(:,1:30));
0140
0141
0142 valnum = x0(1);
0143 valden = x0(2);
0144
0145
0146 kk1 = valnum*kk1;
0147 kk2 = valden*kk2;
0148
0149
0150 kk1mitadr = kk1(18:end);
0151 vmin = min(kk1mitadr);
0152 pos = find(kk1mitadr==vmin);
0153 kk1mitadr(pos(1):end)=vmin;
0154
0155 kk1mitadl = kk1(1:16);
0156 vmin = min(kk1mitadl);
0157 pos = find(kk1mitadl==vmin);
0158 kk1mitadl(1:pos(end))=vmin;
0159
0160 kk1 = [kk1mitadl kk1(17) kk1mitadr];
0161
0162
0163
0164
0165
0166 for i=1:32
0167 RRR(i,:)=kk1(i)/10 + kk1(i)*(c.^2)./(kk2(i)+h11f(i)*c.^2);
0168 end
0169
0170
0171 kkk1=zeros(lcuan);
0172 kkk2=kk1;
0173
0174
0175
0176 [fx,fy]=freqspace(lcuan);
0177
0178 fx=fx*fN;
0179
0180 for i=1:lcuan
0181 for j=1:lcuan
0182 fpunto=F(i,j);
0183 if (fpunto>=max(abs(fx)))
0184 kkk1(i,j)=kk1(lcuan);
0185 kkk2(i,j)=kk2(lcuan);
0186 else
0187
0188 kkk1(i,j)=interp1(fx,kk1,fpunto,'cubic');
0189 kkk2(i,j)=interp1(fx,kk2,fpunto,'cubic');
0190 end
0191 end
0192 end
0193
0194
0195
0196
0197
0198
0199 [fx,fy]=freqspace(tam,'meshgrid');
0200 Fx=fx*fN;
0201 Fy=fy*fN;
0202
0203 [fxa,fya]=freqspace(tam2,'meshgrid');
0204 Fxa=fxa*fN;
0205 Fya=fya*fN;
0206
0207 filtro = interp2(Fx,Fy,kkk1,Fxa,Fya);
0208
0209 beta = interp2(Fx,Fy,kkk2,Fxa,Fya);
0210
0211
0212 [fx,fy]=freqspace(tam2,'meshgrid');
0213 Fx=fx*fN;
0214 Fy=fy*fN;
0215
0216 F=sqrt(Fx.^2+Fy.^2);
0217
0218
0219
0220 fx=zigzag(Fx);
0221 fy=zigzag(Fy);
0222 f=zigzag(F);
0223
0224
0225 ind=find(f==0);
0226
0227 lcuan = tam2;
0228
0229
0230 H=zeros(lcuan^2,lcuan^2);
0231
0232 cont = 1;
0233 for i=1:lcuan
0234 for j=1:lcuan
0235 s = feval(fun,f(cont),factor1,factor2);
0236 hh=exp(-((Fx-fx(cont)).^2+(Fy-fy(cont)).^2)/s.^2);
0237 H(cont,:)=(zigzag(hh))';
0238 cont = cont + 1;
0239 end
0240 end
0241
0242
0243 H(:,ind) = 0;
0244 H(ind,ind)=1;
0245
0246 H = H./ repmat(sum(H,2),1,lcuan^2);
0247
0248
0249
0250
0251 if (prediccion==1)
0252 HH = diag(H);
0253 HH = diag(HH);
0254 H = H - HH;
0255 H(ind,ind)=1;
0256 H = H./ repmat(sum(H,2),1,lcuan^2);
0257 end
0258
0259
0260
0261
0262 k1=zigzag(filtro);
0263
0264 k2=zigzag(beta);
0265
0266