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. Alternativamente, se puede descargar en este enlace yearssn.dat. 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 ventana Wokspace 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])
[M,I]=max(power);
periodo=2*pi/w(I);
text(w(I)+.1,power(I), ['periodo',num2str(periodo)]);
grid on
xlabel('rad/año')
ylabel('Potencia')
title('Espectro de frecuencias')

calculamos el periodo de forma más exacta, utilizando la función max de MATLAB. La función max nos calcula el valor máximo M de un vector y nos devuelve I 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.

Concentración de CO2 en la atmósfera.

La curva de Keeling es una gráfica de la concentración de dióxido de carbono en la atmósfera desde 1958. Las medidas fueron tomadas en Mauna Loa (Hawái) bajo la supervisión de Charles David Keeling. Muestran los rápidos incrementos en los niveles de dióxido de carbono en la atmósfera

Los datos se encuentran en el fichero monthly_in_situ_co2_mlo.csv que se descarga de la página titulada Primary Mauna Loa CO2 Record. También se puede descargar en este enlace monthly_in_situ_co2_mlo.csv

Se guarda este fichero en carpeta de trabajo de MATLAB, se crea un script mediante New Script. En la pestaña Home se pulsa en Import data y se carga el fichero monthy_in_situ_co2_mlo.csv

En Output type seleccionamos Column vectors, indicando que no nos interesa toda la tabla de datos (diez columnas), solamente estamos interesados en la columna E que guarda las concentraciones de CO2 en ppm. Al principio de la columna y al final, hay celdas que contienen el dato -99.99, quiere decir que en ese mes no se ha medido la concentración de CO2. A partir de 1965/01 (enero) y hasta 2020/07 (julio), todas las celdas guardan medidas.

Seleccionamos mediante el ratón o en Range, en la columna E desde E142 hasta E808 ambos inclusive y pulsamos el botón Import selection. En la ventana Workspace cambiamos el nombre de la variable a co, un vector que contiene 667 datos de medidas de las concentraciones efectudas en los meses de esos años.

Representamos gráficamente estos datos, mediante un script. Dentro de la ventana gráfica representamos las medidas de la concentración en cada uno de los meses del año 1965.

plot(co)
grid on
xlim([1,667])
set(gca,'XTick',1:60:667)
set(gca,'XTickLabel',1965:5:2020)
xlabel('Año')
ylabel('Concentración CO_2, ppm')
title('Concentración CO_2, promedio mensual')

axes('Position',[0.21,0.58,2/7,3/10])
hold on
plot(co(1:12));
plot(co(1:12),'ro','markersize',3,'markerfacecolor','r');
hold off
grid on
title('variación mensual')
box off

La causa de las oscilaciones de la concentración se deben a la absorción y a la liberación del CO2 por la vegetación. En el hemisferio norte, en la primavera y verano crecen las plantas que absorben el CO2 de la atmósfera en el proceso de la fotosíntesis. Por otra parte, en el otoño y en el invierno las plantas pierden sus hojas y el el CO2 se libera a la atmósfera. El periodo de oscilación de la concentración de CO2 coincidirá con el tiempo que tarda la Tierra en dar una vuelta alrededor del Sol

Creamos otro script. Importamos el fichero de datos mencionado. En este caso, seleccionamos 28=256 datos de la columna E. Desde E538 hasta E793, desde enero (01) de 1998 hasta abril (04) de 2019. Guardamos los datos en el vector co_1.

El primer paso, es utilizar el comando detrend para eliminar la tendencia creciente de las concentraciones, el segundo paso es aplicar el comando fft, Transformada Rápida de Fourier

t=1:256;
hold on
Fs=1; %frecuencia de muestreo. Inversa del intervalo de tiempo, un mes
f=(0:(length(co_1)-1)/2)*(Fs/length(co_1));
D=detrend(co_1);
FD=fft(D);
Amp=abs(FD(1:length(FD)/2));
plot(f,Amp);
[M,I]=max(Amp); %máximo
plot(f(I),M,'ro','markersize',4,'markerfacecolor','r')
text(f(I)+0.01,M, ['periodo: ',num2str(1/f(I))]);
hold off
grid on
xlabel('frecuencia, 1/mes')
ylabel('Amplitud')
title('Espectro de frecuencias')

Obtenemos un periodo algo mayor de 12 meses.

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)

Referencias

Del apartado titulado, 'Concentración de CO2 en la atmósfera'.

Unofre B Pili, Renante R Violanda. Extracting Earth's orbital period from atmospheric CO2 concentrations using the Fourier transform based on Matlab. Phys. Educ. 55 (2020) 055001