Transformada de Fourier

Para una función no periódica P-->∞

F(ω)= f(t) e iωt dt f(t)= 1 2π F(ω) e iωt dω

La primera integral que obtiene F(ω) se denomina transformada de Fourier de f(t), y la segunda se denomina transformada inversa de Fourier.

El cuadrado f2(t) nos da una idea de cómo la energía contenida en la onda se distribuye en el tiempo, mientas que F2(ω) nos da una idea de como la energía se distribuye en el espectro de frecuencias. Naturalmente,

| f(t) | 2 dt= | F(ω) | 2 dω

Pulso rectangular

Sea un pulso rectangular tal que f(t) es cero excepto en el intervalo [-a,a] que vale A, tal como se muestra en la figura

La transformada de Fourier de f(t) vale

F(ω)= f(t)exp(iωt)·dt F(ω)= a a Aexp(iωt)·dt = 2A ω sin(ωa)

>> syms a t w;
>> int('exp(-i*w*t)',t,-a,a)
ans =(2*sin(a*w))/w

Alternativamente, utilizamos la función fourier de MATLAB

syms a t;
ft=heaviside(t+a)-heaviside(t-a);
Fw=fourier(ft);
Fw=simplify(Fw)

ft=subs(ft,a,1);
subplot(2,1,1)
ezplot(ft,[-2,2]);
ylim([-0.2 1.2])
xlabel('t');
ylabel('f(t)')
title('Pulso rectangular')

Fw=subs(Fw,a,1);
subplot(2,1,2)
hg=ezplot(Fw,[-10,10]);
set(hg,'color','r')
ylim([-1 2.2])
xlabel('\omega');
ylabel('F(\omega)')
title('Transformada de Fourier')
grid on

En la figura vemos un pulso rectangular de semianchura a=1 y su transformada de Fourier

No es necesario utilizar Math Symbolic Toolbox para dibujar la función y su transformada de Fourier como vemos en este script.

t=[-2 -1 -1 1 1 2];
ft=[0 0 1 1 0 0];
subplot(2,1,1)
plot(t,ft,'b');
ylim([-0.2 1.2])
xlabel('t');
ylabel('f(t)')
title('Pulso rectangular')

fw=@(w) 2*sin(w+eps)./w;
subplot(2,1,2)
fplot(fw,[-10,10], 'color','r');
ylim([-1,2.2])
xlabel('\omega');
ylabel('F(\omega)')
grid on

Pulso triangular

f(t)={ 1+ t a at<0 1 t a 0t<a 0otros t

La transformada de Fourier es

F(ω)= f(t)exp(iωt)dt= a 0 ( 1+ t a )exp(iωt)dt + 0 a ( 1 t a )exp(iωt)dt= 4 sin 2 ( ωa /2 ) a ω 2

>> syms t a w;
>> f1=(1+t/a)*(heaviside(t+a)-heaviside(t))*exp(-i*w*t);
>> f2=(1-t/a)*(heaviside(t)-heaviside(t-a))*exp(-i*w*t);
>> Fw=int(f1,t,-a,0)+int(f2,t,0,a)
Fw =1/(a*w^2) - (1/exp(a*w*i))/(a*w^2) - (exp(a*w*i) - 1)/(a*w^2)
>> Fw=simplify(Fw)
Fw =(4*sin((a*w)/2)^2)/(a*w^2)
syms t;
ft=(1-t)*(heaviside(t)-heaviside(t-1))+(1+t)*(heaviside(t+1)-heaviside(t));
Fw=fourier(ft);
Fw=simplify(Fw)

subplot(2,1,1)
ezplot(ft,[-2,2]);
ylim([-0.2 1.2])
xlabel('t');
ylabel('f(t)')
title('Pulso triangular')

subplot(2,1,2)
hg=ezplot(Fw,[-10,10]);
set(hg,'color','r')
ylim([-0.1 1.1])
xlabel('\omega');
ylabel('F(\omega)')
title('Transformada de Fourier')
grid on

En la ventana de comandos aparece la transformada de Fourier, que es la misma que hemos deducido con a=1

>> Fw=-(2*cos(w) - 2)/w^2

Función coseno limitado

Sea la función f(t) definida en el intervalo [-b/2, b/2]

f(t)={ Acos( ω 0 t ), b 2 t b 2 0

Su transformada de Fourier es

F( ω )= f(t)exp( iωt )dt=A b/2 b/2 cos( ω 0 t )exp( iωt )dt A{ b/2 b/2 cos( ω 0 t )cos( ωt )dt+i b/2 b/2 cos( ω 0 t )sin( ωt )dt }

La segunda integral es nula, el intervalo de integración es simétrico y la función es impar

>> syms w w0 b t;
>> int(cos(w0*t)*sin(w*t), t,-b/2, b/2)
ans =0

El resultado de la primera integral es

A b/2 b/2 cos( ω 0 t )cos( ωt )dt= A 2 b/2 b/2 ( cos( ( ω 0 ω )t )+cos( ( ω 0 +ω )t ) )dt= A 2 { sin( ( ω 0 ω )t ) ω 0 ω + sin( ( ω 0 +ω )t ) ω 0 +ω } | b/2 b/2 F( ω )=A sin( ( ω 0 ω ) b 2 ) ω 0 ω + sin( ( ω 0 +ω ) b 2 ) ω 0 +ω

b=pi;
w0=6;
subplot(2,1,1)
fplot(@(t) cos(w0*t),[-b/2,b/2])
line([-4,-b/2],[0,0])
line([b/2,4],[0,0])
grid on
xlabel('t')
ylabel('f(t)')
title('Función, f(t)')
subplot(2,1,2)
f=@(w) sin((w0+w)*b/2)./(w0+w)+sin((w0-w)*b/2)./(w0-w);
fplot(f,[-20,20])
grid on
xlabel('\omega')
ylabel(' F(\omega)')
title('Transformada, F(\omega)')

Función exponencial

Transformada de Fourier de la función f(t)=Aexp(-γ|t|)

Primero, calculamos la transformada de Fourier de la función f(t)=Aexp(-γtu(t). La integración de f(t) se extiende entre 0 e ∞

F(ω)= 0 A e γt e iωt dt= 0 Aexp( (γ+iω)t )dt= A γ+iω exp( (γ+iω)t ) | 0 = A γ+iω γ>0

A continuación, tenemos en cuenta la propiedad de la transformada de Fourier de la función f(-t) es F(-ω), la transformada de las dos exponenciales es la suma

F(ω)= A γ+iω + A γiω =A 2γ γ 2 + ω 2

g=1;
ft=@(t) exp(-g*abs(t));
subplot(2,1,1)
fplot(ft,[-5,5],'color','b');
ylim([-0.1 1.1])
xlabel('t');
ylabel('f(t)')
title('Exponencial')

Fw=@(w) 2*g./(g^2+w.^2);
subplot(2,1,2)
fplot(Fw,[-10,10],'color','r');
ylim([-0.1 2.1])
xlabel('\omega');
ylabel('F(\omega)')
grid on

Oscilación amortiguada

En este ejemplo calculamos la trasnformada de Fourier de la función que describe una oscilación amortiguada que parte del instante t=0.

f(t)=exp(γt)cos( ω 0 t)

Escribimos f(t) en forma equivalente y calculamos su transformada de Fourier

f(t)= 1 2 exp(γt)( exp(i ω 0 t)+exp(i ω 0 t) ) F(ω)= f(t)exp(iωt)dt= 1 2 ( 0 exp( γti(ω ω 0 )t )dt+ 0 exp( γti(ω+ ω 0 )t )dt )= 1 2 ( exp( γti(ω ω 0 )t ) γi(ω ω 0 ) + exp( γti(ω+ ω 0 )t ) γi(ω+ ω 0 ) ) 0 = 1 2 ( 1 γ+i(ω ω 0 ) + 1 γ+i(ω+ ω 0 ) )= 1 2 ( γi(ω ω 0 ) γ 2 + (ω ω 0 ) 2 + γi(ω+ ω 0 ) γ 2 + (ω+ ω 0 ) 2 )

g=1; %amortiguamiento
w0=6; %frecuencia angular
ft=@(t) exp(-g*t).*cos(w0*t);
subplot(2,1,1)
fplot(ft,[0,5],'color','b');
ylim([-1 1.1])
xlabel('t');
ylabel('f(t)')
title('Oscilación amortiguada')

Fw=@(w) abs(((g-1i*(w-w0))./(g^2+(w-w0).^2)+(g-1i*(w+w0))./(g^2+(w+w0).^2))/2);
subplot(2,1,2)
fplot(Fw,[-10,10],'color','r');
ylim([-0.1 0.6])
xlabel('\omega');
ylabel('F(\omega)')
title('Transformada')
grid on

Como F(ω) es complejo, representamos el módulo

Función de Gauss

La función de Gauss es una de las funciones más importantes, se define

f(x)= 1 σ 2π exp( (xμ) 2 2 σ 2 )

donde μ es la media y σ es la desviación estándar de acuerdo con su interpretación estadística.

x=-5:0.05:5;
s2=[0.2,1,5,0.5]; %cuadrado de la desviación estándar
mu=[0,0,0,-2]; %media 
col=['b','g','r','k'];
hold on
for k=1:length(s2)
    y=exp(-(x-mu(k)).^2/(2*s2(k)))/sqrt(s2(k)*2*pi);
    str=sprintf('\\mu=%1.1f, \\sigma^2=%1.1f',mu(k),s2(k));
    plot(x,y,col(k),'displayName',str)
end
hold off
xlabel('x')
ylabel('f(x)')
legend('-DynamicLegend','location','NorthEast')
title('Función de Gauss')
grid on

El máximo de la función se produce para x=μ. Denominamos Δx a la anchura de la función de Gauss para una ordenada igual a la mitad del máximo. Por ejemplo, en la curva de color verde el máximo vale y=0.4, aproximadamente. Determinamos el intervalo de valores de x, Δx cuya ordenada y≥0.2.

1 2 1 σ 2π = 1 σ 2π exp( x 1 2 2 σ 2 ) x 1 =σ 2ln2 Δx=2σ 2ln2

Sea f(t) una función de Gauss

f(t)= 1 σ 2π exp( t 2 2 σ 2 ) f(t)·dt =1

>> syms t;
>> syms s2 positive;
>> f=exp(-t^2/(2*s2))/(sqrt(2*pi*s2));
>> int(f,t,-inf,inf)
ans =1

Para calcular la transformada de Fourier tenemos que partir del resultado de la integral

I= 0 exp(α x 2 )dx = 1 2 π α

que se justifica del siguiente modo:

I= 0 exp(α x 2 )dx I= 0 exp(α y 2 )dx I 2 = 0 0 exp( α( x 2 + y 2 ) )dx·dy x 2 + y 2 = r 2 dx·dy=r·dr·dθ I 2 = 0 π/2 0 exp( α r 2 )r·dr·dθ= π 2 0 exp( α r 2 )r·dr= π 4α I= 0 exp(α x 2 )dx = 1 2 π α

Con ayuda de este resultado calculamos la transformada de Fourier de un pulso descrito por una función de Gauss

F(ω)= f(t)·exp(iωt)·dt = 1 σ 2π exp( t 2 2 σ 2 ) exp(iωt)·dt= 1 σ 2π exp( t 2 2 σ 2 iωt ) ·dt F(ω)= 1 σ 2π exp( σ 2 ω 2 2 ) exp( (t+i σ 2 ω) 2 2 σ 2 ) ·dt F(ω)=exp( σ 2 ω 2 2 )

La tansformada de Fourier de una función de Gauss es otra función de Gauss

>> syms t;
>> syms s positive;
>> ft=exp(-t^2/(2*s^2))/(s*sym('sqrt(2*pi)'));
>> Fw=fourier(ft)
Fw =1/exp((s^2*w^2)/2)

Se sugiere comprobar que cuanto es mayor σ, más ancha es la curva f(t) y más estrecha es la curva F(ω) y viceversa.

sigma=2;
ft=@(t) exp(-t.^2/(2*sigma^2))/(sigma*sqrt(2*pi));
subplot(2,1,1)
fplot(ft,[-7,7],'color','b');
ylim([-0.05 0.25])
xlabel('t');
ylabel('f(t)')
title('Función de Gauss')
grid on

fw=@(w) exp(-w.^2*sigma^2/2);
subplot(2,1,2)
fplot(fw,[-7,7],'color','r');
ylim([-0.2,1.2])
xlabel('w');
ylabel('F(w)')
grid on

Función delta de Dirac δ(t)

La función delta de Dirac tiene la siguiente definición

δ(t)={ 0t0 t=0 δ(t)·dt=1 f(t)δ(ta)·dt=f(a)

Se obtiene a partir de una función de Gauss cuyo parámetro σ tiende a cero

f(t)= 1 σ 2π exp( t 2 2 σ 2 ) lim f(t,σ) σ0 =δ(t)

Calculamos la transformada inversa de Fourier de la función delta de Dirac.

F(ω)=2π·δ( ω ω 0 ) f(t)= δ( ω ω 0 )exp(iωt)dt =exp(i ω 0 t) ω 0 + ω 0 δ( ω ω 0 )dt =exp(i ω 0 t)

F(ω) es cero para todos los valores de ω, excepto para ω=ω0. El área bajo la curva δ(ω-ω0) es la unidad.

>> syms w w0;
>> Fw=2*pi*dirac(w-w0);
>> ft=ifourier(Fw)
ft =exp(w0*x*i)

Teniendo en cuenta que

cos( ω 0 t )= 1 2 ( exp(i ω 0 t)+exp(i ω 0 t) )

la transformada de Fourier F(ω)=π[δ(ω-ω0)+δ(ω+ω0)]

>> syms w w0;
>> Fw=pi*(dirac(w-w0)+dirac(w+w0));
>> ft=ifourier(Fw)
ft =(1/exp(w0*x*i))/2 + exp(w0*x*i)/2
>> simplify(ft) 
ans =cos(w0*x)

De otro modo

f(x)=cos( k 0 x)= 1 2 ( exp(i k 0 x)+exp(i k 0 x) ) F(k)= f(x)exp(ikx)dx= 1 2 ( exp( i(k k 0 )x )dx+ exp( i(k+ k 0 )x )dx )= π( δ(k k 0 )+δ(k+ k 0 ) )

Dos funciones delta de Dirac situadas en +k0 y en -k0.

>> syms k0 x;
>> fx=cos(k0*x);
>> Fk=fourier(fx)
Fk =pi*(dirac(k0 - w) + dirac(k0 + w))

Combinación de las funciones coseno y Gauss

En primer lugar, vamos a obtener el resultado de la integral

I= exp( u 2 2 σ 2 +iuy ) du= exp( ( u 2 σ i yσ 2 ) 2 y 2 σ 2 2 ) du= exp( y 2 σ 2 2 ) exp( ( u 2 σ i yσ 2 ) 2 ) du

Haciendo el cambio de variable

x= u 2 σ i yσ 2 dx= 1 2 σ du

Obtenemos

I= 2 σexp( y 2 σ 2 2 ) exp( x 2 ) dx

Sabiendo el resultado de la integral

exp( x 2 ) dx= π

Obtenemos el resultado

I= 2π σexp( y 2 σ 2 2 )

que utlizaremos más adelante.

Transformada de Fourier

Vamos a obtener la transformada de Fourier de la función

f(x)=exp( ( x x 0 ) 2 2 σ x 2 )cos( k 0 x ) = 1 2 exp( ( x x 0 ) 2 2 σ x 2 )( exp(i k 0 x)+exp(i k 0 x) )

En primer lugar, obtenemos la transformada de Fourier F1(k) de la función f1(x)

f 1 (x)=exp( i k 0 x ( x x 0 ) 2 2 σ x 2 ) F 1 (k)= f(x)exp(ikx)dx = exp( ikx+i k 0 x ) exp( ( x x 0 ) 2 2 σ x 2 )dx F 1 (k)=exp( i(k k 0 ) x 0 ) exp( i( k 0 k)(x x 0 ) ( x x 0 ) 2 2 σ x 2 ) dx

Tenemos una integral del tipo descrito al principio de este apartado

u=x x 0 du=dx y= k 0 kσ= σ x F 1 (k)=exp( i(k k 0 ) x 0 ) exp( iyu u 2 2 σ 2 ) du F 1 (k)= 2π σ x exp( i k 0 x 0 )exp( ik x 0 )exp( ( k 0 k ) 2 σ x 2 2 )

De modo análogo, obtenemos la transformada de Fourier de la función f2(x)

f 2 (x)=exp( i k 0 x ( x x 0 ) 2 2 σ x 2 ) F 2 (k)= f(x)exp(ikx)dx = exp( ikxi k 0 x ) exp( ( x x 0 ) 2 2 σ x 2 )dx F 2 (k)=exp( i(k+ k 0 ) x 0 ) exp( i( k 0 k)(x x 0 ) ( x x 0 ) 2 2 σ x 2 ) dx u=x x 0 du=dx y=( k 0 +k)σ= σ x F 2 (k)= 2π σ x exp( i k 0 x 0 )exp( ik x 0 )exp( ( k 0 +k ) 2 σ x 2 2 )

La transformada de Fourier de la función f(x)= f1(x)+ f2(x) es F(k)=(F1(k)+F2(k))/2

x0=0;
k0=2;
s2=5;
x=-10:0.05:10;
fx=exp(-(x-x0).^2/(2*s2)).*cos(k0*x);
subplot(2,1,1)
plot(x,fx,'b')
grid on
xlabel('x')
ylabel('f(x)')
title('Función')

k=-5:0.05:5;
Fk=sqrt(2*pi*s2)*exp(-1i*k*x0).*(exp(1i*k0*x0)*exp(-(k-k0).^2*s2/2)...
+exp(-1i*k0*x0)*exp(-(k+k0).^2*s2/2))/2;
subplot(2,1,2)
plot(k,abs(Fk),'r')
grid on
xlabel('k')
ylabel('F(k)')
title('Transformada')

En la parte superior de la figura, vemos la función coseno, cuya envolvente (amplitud) es la función de Gauss.

No hay cambio apreciable en el valor absoluto de F(k) cuando se modifica la posición del centro de la función de Gauss, x0 a un valor distinto de cero. La transformada de Fourier de esta función consiste en dos picos centrados en ±k0.

Si incrementamos el parámetro σ2 a una valor grande por ejemplo, 100, vemos que la amplitud de la función f(x) es casi constante, y su transformada de Fourier son dos picos puntiagudos centrados en ±k0 que se van aproximando a una función delta de Dirac a medida que se incrementa σ2

Tabla de transformadas de Fourier

f(t) F(ω)
exp(-atu(t) 1 a+iω
t·exp(-atu(t) 1 ( a+iω ) 2
δ(t) 1
1 δ(ω)
exp(-0t) δ(ω-ω0)
cos(ω0t) π(δ(ω+ω0)+δ(ω-ω0))
sin(ω0t) iπ(δ(ω+ω0)-δ(ω-ω0))
u(t) πδ(ω)+1/()
exp(-at)·sin(ω0tu(t) ω 0 ( a+iω ) 2 + ω 0 2
exp(-at)·cos(ω0tu(t) a+iω ( a+iω ) 2 + ω 0 2
>> ft=sym('1');
>> Fw=fourier(ft)
Fw =2*pi*dirac(w)

>> syms t;
>> ft=heaviside(t);
>> Fw=fourier(ft)
Fw =pi*dirac(w) - i/w
 
>> ft=dirac(t);
>> Fw=fourier(ft)
Fw =1

>> syms a positive;
>> ft=exp(-a*t)*heaviside(t);
>> Fw=fourier(ft)
Fw =1/(a + w*i)

>> ft=t*exp(-a*t)*heaviside(t);
>> Fw=fourier(ft)
Fw =1/(a + w*i)^2

Transformada de Fourier de una función periódica f(t) de semiperiodo P.

f(t)= c k exp( ikπ t P ) F(ω)=2π c k δ( ωk ω 0 ) ω 0 = π P

Transformada de f(t-T)

f(tT)exp( iωt )dt = f(ξ)exp( iω(ξ+T) )dξ = exp( iωT ) f(ξ)exp( iωξ )dξ= exp( iωT )F(ω)

Transformada de exp(0t)f(t)

f(t)exp( i ω 0 t )exp( iωt )dt = f(t)exp( i(ω ω 0 )t )dt =F(ω ω 0 )

Transformada de f(at) para a>0

f(at)exp( iωt )dt = 1 a f(ξ)exp( iωξ/a )dξ = 1 a F( ω a )