Transformada rápida de Fourier (II)

Posición del asteriode Pallas

Gauss estaba interesado en calcular de un modo preciso la órbita del asteroide Pallas, a partir de las observaciones de su posición.

Ascensión 0 30 60 90 120 150 180 210 240 270 300 330
Declinación 408 89 -66 10 338 807 1238 1511 1583 1462 1183 804

Disponía de 12 pares de datos (xi,yi) que deseaba interpolar con seis funciones armónicas.

y= A0+A1cos(ωx)+ B1sin(ωx)+ A2cos(2ωx)+ B2sin(2ωx)+...+ A6cos(6ωx)+ B6sin(6ωx)

Al final de la página Ajuste de datos con funciones armónicas hemos tratado este mismo problema. Ahora, vamos a calcular los coeficientes A0...A6 y B1...B6 a partir de la Transformada Rápida de Fourier

x=0:30:330; 
y=[408 89 -66 10 338 807 1238 1511 1583 1462 1183 804];
hold on
plot(x,y,'ro','markersize',3,'markerfacecolor','r')

n=length(y);
w=2*pi/360;
d=fft(y);
M=floor((n+1)/2);
a0=d(1)/n; 
a=2*real(d(2:M))/n;
a(6)=d(M+1)/n;
b=-2*imag(d(2:M))/n;
b(6)=0;

x=0:360;
k=1:6;
y=a0+a*cos(w*k'*x)+b*sin(w*k'*x);
plot(x,y,'b')
hold off
xlim([0 360])
xlabel('Ascensión (grados)')
ylabel('Declinación (minutos)')
title('Posición del asteroide Pallas')
grid on

La solución encontrada por Gauss es equivalente a la Transformada Rápida de Fourier que fué descubierta mucho más tarde en 1965 por John Turkey y John Cooley

Manchas solares

(USET data/image, Royal Observatory of Belgium, Brussels, http://sidc.oma.be/uset/)

Durante siglos los astrónomos han notado que la superficie del Sol no es constante o uniforme en apariencia y que ciertas regiones oscuras aparecen en posiciones aleatorias con una base cíclica. En 1848, Rudolf Wolfer propuso una regla que combinaba el número y el tamaño de estas manchas solares en un único índice. Usando archivos históricos, los astrónomos aplicaron la regla de Wolfer para determinar la actividad de las manchas solares desde el año 1700.

De esta web y desde la dirección http://sidc.oma.be/sunspot-data/ se pueden descargar los datos de las manchas solares en un fichero ASCII. Desde el enlace yearly sunspot number se ha descargado el fichero yearssn.dat que contiene los datos del número de manchas solares desde el año 1700 hasta el año 2012. Una porción del largo fichero es el siguiente:

Año Indice Wolfer
1700.5  
1701.5  
1702.5 
1703.5 
1704.5  
1705.5  
1706.5 
1707.5  
.......
2007.5 
2008.5  
2009.5  
2010.5  
2011.5 
2012.5
5.0  
11.0  
16.0  
23.0  
36.0  
58.0  
29.0  
20.0  
......
7.5  
2.9  
3.1  
16.5  
55.7  
57

Situamos este fichero en la carpeta por defecto (Current folder). Creamos un script para cargar con el comando load los datos del fichero y representar gráficamente el número de las manchas solares en función del tiempo. Los datos del fichero se guardan en la matriz yearssn de 313×2. Una matriz que tiene el mismo nombre que el fichero. Estas dimensiones las podemos obtener llamando a la función size(yearssn) en la ventana de comandos o en la ventanaWokspace de las variables.

load yearssn.dat
t=yearssn(:,1); %año
m_s=yearssn(:,2); %número de manchas solares

hold on
plot(t,m_s,t,m_s,'bo','markersize',3,'markerfacecolor','b')
p=polyfit(t,m_s,1);
trend=polyval(p,t);
plot(t,trend)
hold off
xlim([1700 2020])
xlabel('Año')
ylabel('Número manchas solares')
grid on
title('Actividad solar')

Se puede observarse fácilmente la tendencia cíclica del fenómeno: los picos y los valles se encuentran separados entre sí algo más de 10 años (podemos seleccionar una parte de la imagen y ampliarla eligiendo los iconos Zoom y Pan debajo del menú de Figure Window). Con Data cursor señalamos dos picos consecutivos y medimos su diferencia en años (10)

La línea de color roja indica la tendencia media de los datos obtenida mediante ajuste por mínimos cuadrados llamando a la función MATLAB polyfit que nos devuelve la pendiente y la ordenda en el origen de la recta de regresión en el vector p. Antes de aplicar transformada de Fourier tenemos que substraer la tendencia lineal de los datos para lo que podemos utilizar alternativamente, la función detrend de MATLAB.

Tomamos como intervalo de tiempo Δt=1 año. Creamos la serie de frecuencias angulares del mismo modo que en el ejemplo anterior, siendo n es el número de datos y la frecuencia crítica o de Nyquist ωc=π/Δt=π rad/año.

Escribimos un script que primero crea la serie temporal (t,f) a partir de los datos originales y después la serie (ω,g) mediante la Transformada Rápida de Fourier fft. Representemos gráficamente la potencia (cuadrado de la amplitud) en una segunda ventana y ampliamos la parte de la imagen alrededor de las frecuencias cuya potencia es máxima con los iconos Zoom y Pan y vemos que el máximo se sitúa un poco más allá de 0.6, aproximadamente ω=0.62 rad/año lo que da un periodo P=2π/ω=10.1 años tal como se ve en la figura.

load yearssn.dat
t=yearssn(:,1); %año
m_s=yearssn(:,2); %número de manchas solares

f=detrend(m_s);
plot(t,f)
xlim([1700 2020])
xlabel('Año')
ylabel('Número manchas solares')
title('Serie temporal')

n=length(f);
dt=1; %año
dw=2*pi/(n*dt);
w=(0:n-1)*dw;
wc=pi/dt; %frecuencia angular crítica
g=fft(f);
power=abs(g).^2;
figure
stem(w,power,'markersize',4,'markerfacecolor','r')
xlim([0 wc])
indice=find(power==max(power));
periodo=2*pi/w(indice(1))
text(w(indice(1))+.1,power(indice(1)), ['periodo',num2str(periodo)]);
grid on
xlabel('rad/año')
ylabel('Potencia')
title('Espectro de frecuencias')

Podemos calcular el periodo de forma más exacta, utilizando las funciones max y find de MATLAB. La función max nos calcula el valor máximo de un vector y find nos devuelve el índice del elemento del vector cuyo valor es máximo. Como hemos visto al principio de esta página, las frecuencias se repiten, por lo que hay dos índices que dan el mismo valor máximo 32 y 283. Solamente es visible el máximo que corresponde al índice 32 que está por debajo de la frecuencia crítica ωc. Con el dato de la frecuencia angular correspondiente a este índice, calculamos el periodo y lo mostramos en la ventana gráfica mediante la función gráfica text.

Ejemplos en el curso de Física

El sonido

Análisis armónico de las mareas

Análisis de Fourier de los desplazamiento de las olas

Resumen

dt intervalo de tiempo dw=2π/(n·dt) intervalo de frecuencias angulares
t=(0:n-1)·dt vector de tiempos w=(0:n-1)·dw vector de frecuencias angulares
(tj,xj) serie temporal g=fft(x) Transformada Rápida de Fourier. Los elementos del vector g son números complejos
n=length(x) número de datos abs(g).^2 potencia
fs=1/dt frecuencia de muestreo ωc=π/dt frecuencia angular crítica (Nyquist)