En esta entrada voy a hablar brevemente de la transformada de Hilbert.
Donde podemos querer aproximar la transformada de Hilbert de una determinada función periódica.
Para ello lo que hacemos es utilizar un spline cúbico y calcular las integrales usando la cancelación que viene de la integral en valor principal.
Dada la función seno (rojo en la imagen) nuestro código aproxima su transformada de Hilbert (azul en la imagen).
function splin3 = spline3(U,E,N)
Función que calcula el spline cúbico de la función U, con N puntos espaciados E.
A=zeros(3*N,3*N);
for k=2:N-1
A(k,k)=1;
A(k,N+k)=E(k);
A(k,2*N+k)=E(k)*E(k);
A(N+k,k)=1;
A(N+k,k+1)=-1;
A(N+k,N+k)=2*E(k);
A(N+k,2*N+k)=3*E(k)*E(k);
A(2*N+k,N+k)=1;
A(2*N+k,N+k+1)=-1;
A(2*N+k,2*N+k)=3*E(k);
end
los bordes
1
k=0;
A(1,1)=1;
A(1,N+1)=E(1);
A(1,2*N+1)=E(1)*E(1);
A(N+1,1)=1;
A(N+1,2)=-1;
A(N+1,N+1)=2*E(1);
A(N+1,2*N+1)=3*E(1)*E(1);
A(2*N+1,N+1)=1;
A(2*N+1,N+2)=-1;
A(2*N+1,2*N+1)=3*E(1);
A(N,N)=1;
A(N,2*N)=E(N);
A(N,3*N)=E(N)*E(N);
A(2*N,N)=1;
A(2*N,1)=-1;
A(2*N,2*N)=2*E(N);
A(2*N,3*N)=3*E(N)*E(N);
A(3*N,2*N)=1;
A(3*N,N+1)=-1;
A(3*N,3*N)=3*E(N);
A=inv(A);
B=zeros(3*N,1);
for k=2:N-1
B(k)=(U(k)-U(k-1))/E(k);
end
los bordes
1
k=0;
B(1)=(U(1)-U(N))/E(1);
B(N)=(U(N)-U(N-1))/E(N);
splin3=A*B;
Ahora el código que calcula la transformada de Hilbert:
function HH=hilbertt(Z,x,E,N,tol2)
S=spline3(Z,E,N);
HH=zeros(N,1);
Hilbert Transform
Corregido y bien
for i=1:N/2 %Izquierda del dominio
%ILFT int_{x(i-1)}^(x(i)) usando p(i); integral singular
j=i;
if i==1
a=Z(N);%periodicity Z(N)=Z(0)
XX=x(i)-E(i);
else
a=Z(j-1);
XX=x(j-1);
end
aa=Z(i);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(ILFTH,0,1,tol2);
for k=1:N/2-1
Parte izquierda de la integral
ban=i-1-k; integral en el dominio por debajo de x(i-1)
if ban<0
j=N+ban+1; p(j)
a=Z(j-1);
aa=Z(i);
XX=x(j-1)-N*E(1);Nodos uniformes
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);
else
j=ban+1;
if j==1
a=Z(N);periodicidad Z(x(0))
XX=x(1)-E(1);
else
a=Z(j-1);
XX=x(j-1);
end
aa=Z(i);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);
end
end
IRGH int_{x(i)}^(x(i+1)) using p(i+1)
j=i+1;
a=Z(j-1);
aa=Z(i);
XX=x(j-1);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(IRGHH,0,1,tol2);
for k=1:N/2-1 Parte de la derecha del dominio
j=i+1+k; integral en el dominio por encima de x(i+1)
a=Z(j-1);
aa=Z(i);
XX=x(j-1);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);
end
end
for i=N/2+1:NParte de la derecha del dominio
ILFT int_{x(i-1)}^(x(i)) usando p(i); integral singular
j=i;
a=Z(j-1);
XX=x(j-1);
aa=Z(i);
if i==N
b=S(1);partial Z(x(0))+periodicidad
else
b=S(i+1);partial_x Z(x(i)) usando p(i+1)
end
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(ILFTH,0,1,tol2);
for k=1:N/2-1
Parte de la izquierda de la integral no singular
j=i-k; integral en el dominio por debajo de x(i-1),
usando p(i-1),p(i-2)…
a=Z(j-1);
XX=x(j-1);
aa=Z(i);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);%
end
IRGH int_{x(i)}^(x(i+1)) usando p(i+1)
if i==N %p(i+1) en este caso es p(1) por la periodicidad
j=1;
a=Z(N);
aa=Z(i);
XX=x(N);
b=S(1);%partial_x Z(x(N)) por periodicidad
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(1);
X=x(i);
else
j=i+1;
a=Z(j-1);
aa=Z(i);
XX=x(j-1);
b=S(i+1);partial_x Z(x(i))
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
end
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+EE*quadl(IRGHH,0,1,tol2);
con EE evitamos E(j)=E(N+1).
for k=1:N/2-1 parte de la derecha de la integral no singular
ban=i+1+k;
if ban<=N
j=i+1+k; integral para el dominio por encima de x(i+1)
a=Z(j-1);
XX=x(j-1);
aa=Z(i);
b=S(i+1);partial_x Z(x(i)) usando p(i+1)
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);
else
j=ban-N;
if j==1
a=Z(N);periodicidad
XX=x(N);
else
a=Z(j-1);
XX=x(j-1)+N*E(1);Nodos espacialmente uniformes
end
aa=Z(i);
b=S(i+1);partial_x Z(x(i)) usando p(i+1)
cc=S(j);
d=S(N+j);
e=S(2*N+j);
EE=E(j);
X=x(i);
save(‘integrando’,'aa’,'a’,'b’,'cc’,'d’,'e’,'EE’,'X’,'XX’)
HH(i)=HH(i)+E(j)*quadl(INSH,0,1,tol2);%
end
end
end
Se habrán fijado que faltan 3 archivos más…
Aquí están:
function v=INSH(y)
load(‘integrando’)
zz=a+EE*y.*(cc+EE*y.*(d+e*EE*y));
x2=X-XX-EE*y;
v=(1/(2*pi))*(-aa+zz)./(tan(x2/2));
function v=ILFTH(y)
load(‘integrando’)
ff=cc+d*EE*(1+y)+e*(EE^2)*(1+y.*(1+y));
v=(1/(2*pi))*(-ff).*(2-(1/6).*(EE*(1-y)).^2);
function v=IRGHH(y)
load(‘integrando’)
ff=cc+(EE*y.*(d+e*EE*y));
v=(1/(2*pi))*(-ff).*(2-(1/6).*(EE*y).^2);
–Nota: Este código está pensado como una caja negra… es decir,
úsalo si lo necesitas pero no te preocupes mucho de cómo funciona.
Eso es así porque tratar de explicar un código tan largo (pero no difícil)
se sale del objetivo de la entrada.
