Home > PR > constantes_sin_continua.m

constantes_sin_continua

PURPOSE ^

function [H,k1,k2]=constantes_sin_continua(fs,tam2,prediccion,fun,factor1,factor2);

SYNOPSIS ^

function [H,k1,k2]=constantes_sin_continua(fs,tam2,prediccion,fun,factor1,factor2);

DESCRIPTION ^

 function [H,k1,k2]=constantes_sin_continua(fs,tam2,prediccion,fun,factor1,factor2);
  
  Constantes necesarias para obtener la respuesta del modelo.
  Si esto se va a utilizar para predecir un coeficiente a partir
  de los vecinos la matriz con las interacciones deberia tener un 
  cero en la posicion correspondiente (en la diagonal)

 
   USO:

      [H,k1,k2]=constantes2(64,32,1);

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 % function [H,k1,k2]=constantes_sin_continua(fs,tam2,prediccion,fun,factor1,factor2);
0002 %
0003 %  Constantes necesarias para obtener la respuesta del modelo.
0004 %  Si esto se va a utilizar para predecir un coeficiente a partir
0005 %  de los vecinos la matriz con las interacciones deberia tener un
0006 %  cero en la posicion correspondiente (en la diagonal)
0007 %
0008 %
0009 %   USO:
0010 %
0011 %      [H,k1,k2]=constantes2(64,32,1);
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 % Dominio de contrastes
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 %              C
0032 % R(C) = Ro + int grad(R(C)) dC
0033 %              0
0034 %
0035 % Tomando Ro=0 :
0036 
0037 R=cumsum(IAFet')';
0038 
0039 
0040 mm = maxi(R);
0041 
0042 R=R/maxi(R);
0043 
0044 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0045 %                          %
0046 %  Respuesta Naka-Rushton  %
0047 %                          %
0048 %%%%%%%%%%%%%%%%%%%%%%%%%%%%
0049 
0050 
0051 
0052 % Dominio de frecuencias de Fourier
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 % Reordenacion 1D del dominio
0062 
0063 fx=zigzag(Fx);
0064 fy=zigzag(Fy);
0065 f=zigzag(F); 
0066 
0067 % Posicion en la que se encuentra la frecuencia 0
0068 
0069 ind=find(f==0);
0070 
0071    
0072 % Calculo del kernel de interaccion (que se aplicara sobre la ordenacion 1D de la señal).
0073 % Ahora no debe tener ceros
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 % Esto hay que hacerlo ya que el dominio es de Fourier
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 % CSF
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 %opciones=optimset('Display','final','MaxFunEvals',25000,'MaxIter',25000,'TolFun',1e-20,'TolX',1e-20)
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 % Esto es solo por si se quiere depurar
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 % Esto es solo por si se quiere depurar
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 % a partir de una fila central costruimos la matriz con simetria de revolución
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 % Esto es para interpolar a otros tamaños de bloque, simplemente por
0196 % no tocar lo anterior que ya esta ajustado para bloques de 32x32
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 % Reordenacion 1D del dominio
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 % Si es para precedir hay que poner la diagonal a cero para no usar
0249 % el que se quiere predecir
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

Generated on Wed 29-Nov-2006 16:19:19 by m2html © 2003