Series de Fourier

Funciones armónicas

En primer lugar, vamos a distinguir entre las magnitudes: amplitud A, frecuencia f y fase φ en la función armónica

x=Asin(2πf·t+φ).

subplot(2,1,1)
x=@(t) 10*sin(2*pi*0.1*t); %amplitud 10
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
title('Distinta amplitud')
ylim([-11,11])
grid on

subplot(2,1,2)
x=@(t) 5*sin(2*pi*0.1*t); %amplitud 5
fplot(x,[0,50])
ylim([-10,10])
xlabel('t(ms)')
ylabel('x')
ylim([-11,11])
grid on

subplot(2,1,1)
x=@(t) 10*sin(2*pi*0.1*t); %frecuencia, 100 Hz
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
title('Distinta frecuencia')
ylim([-11,11])
grid on

subplot(2,1,2)
x=@(t) 10*sin(2*pi*0.2*t);  %frecuencia, 200 Hz
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
ylim([-11,11])
grid on

subplot(4,1,1)
x=@(t) 10*sin(2*pi*0.1*t);
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
title('Distinta fase inicial')
ylim([-11,11])
grid on

subplot(4,1,2)
x=@(t) 10*sin(2*pi*0.1*t+pi/2);
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
ylim([-11,11])
grid on

subplot(4,1,3)
x=@(t) 10*sin(2*pi*0.1*t+pi);
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
ylim([-11,11])
grid on

subplot(4,1,4)
x=@(t) 10*sin(2*pi*0.1*t+3*pi/2);
fplot(x,[0,50])
xlabel('t(ms)')
ylabel('x')
ylim([-11,11])
grid on

x=200sin(2πf·100+π/2)+100sin(2πf·200+π)+100sin(2πf·400+3π/2)

f=[100,200,400]; %frecuencias
A=[200,100,100]; %amplitudes
phi=[90,180,270]; %fases

subplot(2,2,1)
stem(f,A)
axis([0,500,0,210])
xlabel('Frecuencia')
ylabel('Amplitud')

subplot(2,2,2)
stem(f,phi)
axis([0,500,0,360])
xlabel('Frecuencia')
set(gca,'YTick',0:90:360)
set(gca,'YTickLabel',{'0',90','180','270','360'})
ylabel('Fase')

subplot(2,2,3:4) %resultante
t=(0:0.1:30)/1000; %milisegundos
x=zeros(1,length(t));
for i=1:length(f)
    x=x+A(i)*sin(2*pi*f(i)*t+phi(i)*pi/180);
end
plot(t,x,'r')
xlabel('t(ms)')
ylabel('x')
title('Resultante')
ylim([-410,410])
set(gca,'XTick',(0:5:30)/1000)
set(gca,'XTickLabel',{'0','5','10','15','20','25','30'})
grid on

Función periódica

Una función es periódica de periodo P si hay un número P>0 tal que f(t+P)=f(t). Cualquier múltiplo n entero de P es también periodo f(t+nP)=f(t)

La función f(t)=cos(2πt)+cos(4πt)/2, es la suma de dos funciones periódicas de periodos 1 y 0.5, respectivamente. Como vemos en la gráfica f(t) es periódica con periodo P=1.

Las funciones cos(t) y cos( 2 t ) son periódicas de periodo 2π y 2π/ 2 respectivamente, pero la suma

f(t)=cos( t )+cos( 2 t )

no es periódica.

x=@(t) cos(2*pi*t)+cos(4*pi*t)/2;
subplot(2,1,1)
fplot(x,[0,10]);
xlabel('t')
ylabel('x')

subplot(2,1,2)
x=@(t) cos(2*pi*t)+cos(2*pi*sqrt(2)*t);
fplot(x,[0,10]);
xlabel('t')
ylabel('x')

Serie de Fourier

A primera vista, parece que el problema de analizar formas de ondas complejas representa una tarea formidable. Sin embargo, si la forma de la onda es periódica, se puede representar con una precisión arbitraria, mediante la superposición de un número suficientemente grande de ondas senoidales que forman una serie armónica.

Toda función f(t) periódica de periodo 2P, se puede representar en forma de una suma infinita de funciones armónicas, es decir,

f(t)= a 0 2 + k=1 ( a k cos kπt P + b k sin kπt P )

donde a0 a1 ...ak ... y b1 b2 .... bk .... son los denominados coeficientes de Fourier.

Teniendo en cuenta los resultados de las integrales

P P cos mπt P sin nπt P dt= P π π π cos(mz) sin(nz)dz=0

>> syms t m n;
>> assume(n,'integer');
>> assume(m,'integer');
>> int('sin(m*t)*cos(n*t)',t,-pi,pi)
ans =0

P P cos mπt P cos nπt P dt= P π π π cos(mz) cos(nz)dt P 2π π π ( cos( (m+n)z )+cos( (mn)z ) ) ·dz={ 0mn Pm=n

>> y=int('cos(m*t)*cos(n*t)',t,-pi,pi)
y=(2*(m*cos(pi*n)*sin(pi*m) - n*cos(pi*m)*sin(pi*n)))/(m^2 - n^2)
>> limit(y,m,n)
ans =pi*cos(pi*n)^2

Como n es entero, la respuesta es pi (π)

P P sin mπt P sin nπt P dt= P π π π sin(mz) sin(nz)dz P 2π π π ( cos( (m+n)z )+cos( (mn)z ) ) ·dz={ 0mn Pm=n

>> y=int('sin(m*t)*sin(n*t)',t,-pi,pi)
y =-(2*(m*cos(pi*m)*sin(pi*n) - n*cos(pi*n)*sin(pi*m)))/(m^2 - n^2)
>> limit(y,m,n)
ans =pi*cos(pi*n)^2

Los coeficientes del desarrollo en serie valen

a k = 1 P P P f(t) cos kπt P dt b k = 1 P P P f(t) sin kπt P dt

La suma parcial de las series de Fourier es

s n (t)= a 0 2 + k=1 n ( a k cos kπt P + b k sin kπt P )

Si la función f(t) tiene simetría, algunos de los coeficientes resultan nulos.

Función par

Por ejemplo, para el pulso rectangular simétrico de anchura 1, y periodo 2 se obtienen los siguientes coeficientes.

a 0 = 1 1 0.5 0.5 dt =1 a k = 1 1 0.5 0.5 cos( kπt )dt= 2 kπ ( sin( kπ 2 ) ){ 0kpar 2 kπ ( 1 ) (k1)/2 kimpar

Vamos a reconstruir la función f(t) a partir del desarrollo en serie de Fourier.

s n (t)= 1 2 + 2cos(πt) π 2cos(3πt) 3π + 2cos(5πt) 5π 2cos(7πt) 7π +...

Elaboramos un script en el que

n=7;
hold on
x=[-1 -0.5 -0.5 0.5 0.5 1];
y=[0 0 1 1 0 0];
plot(x,y,'b','linewidth',1.5)
x=linspace(-1,1,100);
y=zeros(length(x),1);
for i=1:length(x)
    y(i)=1/2;
    for k=1:2:n
        y(i)=y(i)+(-1)^((k-1)/2)*2*cos(k*pi*x(i))/(k*pi);
     end
end
plot(x,y, 'r');
title(sprintf('Aproximación de Fourier: %i términos',n))
xlabel('t'); 
ylabel('f(t)')
grid on
hold off

Fenómeno de Gibbs

Aproximamos el pulso rectangular con 50 términos y con 200 términos. Cabría esperar que la serie convergiese a la función cuando se incrementa el número de términos, esto ocurre en casi todos los puntos excepto en las discontinuidades, cuando la función pasa de 0 a 1 en x=-0.5 y de 1 a 0 en x=0.5.

Se ha utilizado, las herramientas Zoom y Pan de la ventana gráfica para mostrar un trozo de la función próximo a x=-0.5

Función impar

Sea ahora la función

Es una función impar, los coeficientes ak son nulos

b k = 1 0 sin( kπt )dt 0 1 sin( kπt )dt = 1 kπ ( 2+2cos(kπ) )={ 0kpar 4 kπ kimpar

El desarrollo en serie es

s n (t)= 4sin(πt) π 4sin(3πt) 3π 4sin(5πt) 5π 4sin(7πt) 7π +...

Escribimos un script en el que

n=7;
hold on
x=[-1 -1 0 0 1 1];
y=[0 1 1 -1 -1 0];
plot(x,y,'b','linewidth',1.5)
x=linspace(-1,1,100);
y=zeros(length(x),1);
for i=1:length(x)
    y(i)=0;
    for k=1:2:n
        y(i)=y(i)-4*sin(k*pi*x(i))/(k*pi);
     end
end
plot(x,y, 'r');
title(sprintf('Aproximación de Fourier: %i términos',n))
xlabel('t'); 
ylabel('f(t)')
grid on
hold off

Desarrollo en serie de Fourier con MATLAB

La función heaviside de MATLAB

En MATLAB la función heaviside(t) se define del siguiente modo:

u ( t ) = { 0 t < 0 0.5 t = 0 1 t > 0

Vamos a elaborar un script que nos permita obtener el desarrollo en serie de Fourier de una función f(t) periódica de periodo 2P.

Definimos primero, la función escalón

f(t)={ 0t<a 1ta

Describimos la función escalón mediante la llamada a la función heaviside(t-a)

La función pulso rectangular del primer ejemplo,

f(t)={ 0t<a 1ata 0t>a

la definimos en MATLAB como heaviside(t+a)-heaviside(t-a), tal como vemos en la figura

La función f(t).

f(t)={ 0t<1 11t<0 10t<1 0t1

Se escribe en MATLAB como, (heaviside(t+1)-heaviside(t))+(-1)*(heaviside(t)-heaviside(t-1))= heaviside(t+1)-2*heaviside(t)+heaviside(t-1)

>> syms t;
>> x=heaviside(t+1)-2*heaviside(t)+heaviside(t-1);
>> ezplot(x,[-2,2]);
>> xlabel('t'),
>> ylabel('f(t)')
>> title('Función heaviside')
>> grid on

La función f(t) cualesquiera que toma valores distintos de cero en el intervalo [a,b] se escribe

( u ( t a ) u ( t b ) ) f ( t ) = { 0 t < a f ( t ) a t b 0 t > b

Función par

Desarrollo en serie de Fourier de la primera función, el pulso rectangular

syms t;
P=1; %semiperiodo
f=heaviside(t+0.5)-heaviside(t-0.5);
a=@(k) int(f*cos(k*pi*t/P),t,-P,P)/P;
b=@(k) int(f*sin(k*pi*t/P),t,-P,P)/P;
fs=a(0)/2;
for k=1:10 %número de términos
    fs=fs+a(k)*cos(k*pi*t/P)+b(k)*sin(k*pi*t/P);
end
pretty(fs)

Otra forma alternativa, utilizando la función symsum de MATLAB

syms t k P n;
assume(k,'Integer')
a=@(f,t,k,P) int(f*cos(k*pi*t/P),x,-P,P)/P;
b=@(f,t,k,P) int(f*sin(k*pi*t/P),x,-P,P)/P;
fs=@(f,t,n,P) a(f,t,0,P)/2+symsum(a(f,t,k,P)*cos(k*pi*t/P)
+b(f,t,k,P)*sin(k*pi*t/P),k,1,n);
f=heaviside(t+0.5)-heaviside(t-0.5);
P=1; %semiperiodo
pretty(fs(f,t,10,P))

En la ventana de comandos vemos el desarrollo en serie de Fourier de esta función periódica para n=10 términos

2cos(pi t)   2cos(3pi t)   2cos(5pi t)   2cos(7pi t)   2cos(9pi t) 
--------- - ---------- + ------------- - ---------- + --------- + 1/2 
    pi        3pi           5pi            7pi            9pi

Función impar

Desarrollo en serie de Fourier de la segunda función

....
f=heaviside(t+1)-2*heaviside(t)+heaviside(t-1);
P=1; %semiperiodo
pretty(fs(f,x,10,P))

En la ventana de comandos vemos el desarrollo en serie de Fourier de esta función periódica para n=10 términos

4sin(pi t)   4sin(3pi t)   4sin(5pi t)   4sin(7pi t)   4 sin(9pi t) 
- -------- - ---------- - ------------- - ----------- - ------------- 
    pi            3 pi         5 pi        7 pi            9 pi

Ejemplo 1

Sea la función periódica f(t)=t entre -1 y 1

Para estudiar otra función periódica basta cambiar la definición de la función f y el valor de su semiperiodo P. Cambiando el valor de N se muestran más o menos términos del desarrollo en serie.

syms t k P n;
assume(k,'Integer')
a = @(f,t,k,P) int(f*cos(k*pi*t/P),t,-P,P)/P;
b = @(f,t,k,P) int(f*sin(k*pi*t/P),t,-P,P)/P;
fs=@(f,t,n,P) a(f,t,0,P)/2+symsum(a(f,t,k,P)*cos(k*pi*t/P)
+b(f,t,k,P)*sin(k*pi*t/P),k,1,n);

%definición de la fuerza y su semiperiodo P
f=t;
P=1;

N=6; %términos del desarrollo en serie
pretty(fs(f,t,N,P))

hold on
ezplot(f,[-P P])
hg=ezplot(fs(f,t,N,P),[-P P]);
set(hg,'color','r')
hold off
xlabel('t')
ylabel('f(t)')
title('Serie de Fourier')

%armónicos
figure
k=1:N;
ak=a(f,t,k,P);
bk=b(f,t,k,P);
subplot(2,1,1)
stem(ak)
xlabel('k');
ylabel('a(k)')

subplot(2,1,2)
stem(bk)
xlabel('k');
ylabel('b(k)')

2sin(pi t)   sin(2pi t)   2 sin(3pi t)   sin(4pi t)  2sin(5pi t) sin(6pi t) 
----------- - ----------- + ---------- - -------- + ---------- - ----------- 
    pi           pi            3 pi         2 pi        5 pi        3 pi

En la gráfica se representa la contribución de los N=6 primeros armónicos cada uno de los armónicos, en la parte superior los coeficientes ak, y en la parte inferior los coeficientes bk. Por ser f(t)=t una función impar, ak=0,

Ejemplo 2

Sea la función

f(t)={ 0π<t<0 t0t< π 2 π 2 π 2 tπ

syms t k P n;
assume(k,'Integer')
a = @(f,t,k,P) int(f*cos(k*pi*t/P),t,-P,P)/P;
b = @(f,t,k,P) int(f*sin(k*pi*t/P),t,-P,P)/P;
fs=@(f,t,n,P) a(f,t,0,P)/2+symsum(a(f,t,k,P)*
cos(k*pi*t/P)+b(f,t,k,P)*sin(k*pi*t/P),k,1,n);

%definición de la fuerza y su semiperiodo P
f=t*(heaviside(t)-heaviside(t-pi/2))+
(heaviside(t-pi/2)-heaviside(t-pi))*pi/2;
P=pi;

N=6; %términos del desarrollo en serie
latex(fs(f,t,N,P))
hold on
ezplot(f,[-P P])
hg=ezplot(fs(f,t,N,P),[-P P]);
set(hg,'color','r')
set(gca,'XTick',-pi:pi/2:pi)
set(gca,'XTickLabel',{'-\pi','-\pi/2','0','\pi/2','\pi'})
set(gca,'YTick',0:pi/2:pi)
set(gca,'YTickLabel',{'0','\pi/2','\pi'})
hold off
grid on
xlabel('t')
ylabel('f(t)')
title('Serie de Fourier')

%armónicos
figure
k=1:N;
ak=a(f,t,k,P);
bk=b(f,t,k,P);
subplot(2,1,1)
stem(ak)
xlabel('k');
ylabel('a(k)')

subplot(2,1,2)
stem(bk)
xlabel('k');
ylabel('b(k)')

Utilizamos la función latex en combinación con el programa MathType para presentar los resultados del desarrollo en serie de Fourier

3π 16 sin( 2t ) 4 sin( 4t ) 8 sin( 6t ) 12 cos( 2t ) 2π cos( 3t ) 9π cos( 5t ) 25π cos( 6t ) 18π cos( t ) π + sin( t )( π 2 +1 ) π + sin( 3t )( 3π 2 1 ) 9π + sin( 5t )( 5π 2 +1 ) 25π

La forma compleja de las series de Fourier

Teniendo en cuenta las relaciones

e iωt =cosωtisinωt e iωt =cosωt+isinωt

El desarrollo en serie de Fourier

f(t)= a 0 2 + k=1 ( a k cos( kπt P )+ b k sin( kπt P ) )

se expresa de la siguiente forma alternativa

f(t)= a 0 2 + k=1 ( a k 2 ( e ikπt/P + e ikπt/P )+ b k 2i ( e ikπt/P e ikπt/P ) ) c 0 = a 0 2 c k = a k i b k 2 c k = a k +i b k 2 f(t)= c 0 + k=1 ( c k e ikπt/P + c k e ikπt/P ) f(t)= c k exp( ikπ t P ) c k = 1 2P P P f(t)exp( ikπ t P )dt

donde P es el semiperiodo de la función periódica de periodo 2P

Pulso rectangular

P=π c k = 1 2π a a A e ikt dt= A kπ sin( ka )={ k=0 c 0 = A π a k0 c k = A kπ sin( ka )

Escribimos un script en el que

a=1; %Semianchura del pulso rectangular a<pi
n=7; %número de terminos
hold on
x=[-pi -a -a a a pi];
y=[0 0 1 1 0 0];
plot(x,y,'b','linewidth',1.5)
x=linspace(-pi,pi,100);
y=zeros(length(x),1);

for j=1:length(x)
    y(j)=0;
    for k=-n:n
        if k==0 
            y(j)=y(j)+a/pi;
        else
            y(j)=y(j)+sin(k*a)*exp(1i*k*x(j))/(k*pi);
        end
     end
end
%ejes
plot([-4 4],[0 0],'k')
plot([0 0],[-0.2 1.2],'k')
%serie de Fourier
plot(x,real(y),'r');
title(sprintf('Aproximación de Fourier: %i términos',n))
xlabel('t'); 
ylabel('f(t)')
grid on
hold off

Pulso diente de sierra

P=π c k = 1 2π π π A π t e ikt dt= A kπ icos( kπ )= A kπ ( 1 ) k i

Escribir un script en el que

n=7; %Número de términos
hold on
plot([-pi,pi],[-1 1],'b','linewidth',1.5)
x=linspace(-pi,pi,100);
y=zeros(length(x),1);
for j=1:length(x)
    y(j)=0;
    for k=-n:n
        y(j)=y(j)+(-1)^k*1i*exp(1i*k*x(j))/(k*pi);
    end
end
%ejes
plot([-4 4],[0 0],'k')
plot([0 0],[-1.5 1.5],'k')
%serie de Fourier
plot(x,real(y), 'r');
title(sprintf('Aproximación de Fourier: %i términos',n))
grid on
xlabel('t'); 
ylabel('f(t)')
hold off