A amostragem é definida como, $$ f_s = \frac{1}{dt} $$ em que $dt$ é a variação de tempo entre as observações
pelos dados abaixo, $dt$ médio das RRLyraes é $6.7884$, ou seja, a amostragem, ou sampling rate, é $$f_s = \frac{1}{6.7884} = 0.1473$$
Agora, vou gerar um sinal senoidal simples com essa amostragem,
In [1]:
'''
gerar um sinal senoidal com a amostragem acima e o periodo médio da RRLyraes
'''
tf = 4539.4593
ti = 2152.5019
n = 351.6
dt = (tf-ti)/n
fs = 1/dt
time = np.arange(ti, tf, 1/fs)
ttime = np.arange(ti,tf,dt/4)
p= 0.576
w = 2*np.pi/p
#Modificando a amostragem
Ny_f = fs/2
time2 = np.arange(ti, tf, 1/Ny_f)
#plot(time2, np.sin(w*time2), 'ko', time2, np.sin(w*time2), 'r-')
#xlim(2500,3000)
subplot(211) #dados sintéticos com a amostragem padrao
plot(time, np.sin(w*time), 'ro', ttime, np.sin(w*ttime), 'b-')
xlim(2500,3000)
subplot(212) # dadis sintéticos com a frequencia de Nyquist (metade da amostragem padrao)
plot(time2, np.sin(w*time2), 'ko', ttime, np.sin(w*ttime), 'r-')
xlim(2500,3000)
Out[1]:
De acordo com o artigo http://acta.astrouw.edu.pl/Vol59/n1/pdf/pap_59_1_1.pdf, 95% da RRLyr AB possuem periodos na faixa entre 0.45 ~ 0.75 dias com a média dos periodos sendo $<P_{ab}> 0.576$ dias.
A ideia é descobrir a quantidade média de pontos de observação para as estrelas para descobrir uma amostragem média.
A partir dos dados das estrelas, vou fazer os histogramas para determinar o ponto inicial e final de observações e quantidade de pontos entre essas observações.
In [71]:
#Pontos iniciais
%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
dados = np.array(np.loadtxt('/home/glauffer/Dropbox/FURG/final_project/data/rrab_ti_tf_pontos.txt'))
#a, b, c = plt.hist(dados.T[0], bins=30)
#plt.xlim(400, 2600)
fig,(ax,ax2) = plt.subplots(1, 2, sharey=True)
a, b, c = ax.hist(dados.T[0], bins=90)
ax2.hist(dados.T[0], bins=90)
ax.set_xlim(0,1000) # most of the data
ax2.set_xlim(2000,2500) # outliers only
# hide the spines between ax and ax2
ax.spines['right'].set_visible(False)
ax2.spines['left'].set_visible(False)
ax.yaxis.tick_left()
ax.tick_params(labeltop='off') # don't put tick labels at the top
ax2.yaxis.tick_right()
# Make the spacing between the two axes a bit smaller
plt.subplots_adjust(wspace=0.15)
d = .015 # how big to make the diagonal lines in axes coordinates
# arguments to pass plot, just so we don't keep repeating them
kwargs = dict(transform=ax.transAxes, color='k', clip_on=False)
ax.plot((1-d,1+d),(-d,+d), **kwargs) # top-left diagonal
ax.plot((1-d,1+d),(1-d,1+d), **kwargs) # bottom-left diagonal
kwargs.update(transform=ax2.transAxes) # switch to the bottom axes
ax2.plot((-d,d),(-d,+d), **kwargs) # top-right diagonal
ax2.plot((-d,d),(1-d,1+d), **kwargs) # bottom-right diagonal
#plt.title('Data Juliana $t_i$ RRab, max = %s'%b[np.argmax(a)])
#plt.suptitle('HJD $t_i$ max = %s'%b[np.argmax(a)], fontsize = 18)
plt.suptitle('Histograma - $t_i$', fontsize = 18)
#ax.set_ylabel('Frequência', fontsize = 14)
#ax2.set_xlabel('Tempo Inicial')
fig.text(0.5, 0.0, 'Tempo Inicial [JD]', ha='center', fontsize = 14)
fig.text(0.0, 0.5, 'Frequência', va='center', rotation='vertical', fontsize = 14)
b[np.argmax(a)]
#plt.show()
#plt.savefig('hist_ti.png', dpi=100)
Out[71]:
In [72]:
#Pontos finais
a, b, c = plt.hist(dados.T[1], bins=300)
plt.xlim(4300, 4600)
#plt.title('Histograma $t_f$ RRab, max = %s'%b[np.argmax(a)])
plt.title('Histograma - $t_f$', fontsize=18)
plt.xlabel('Tempo Final [JD]', fontsize=14)
plt.ylabel('Frequência', fontsize=14)
#plt.title('HJD $t_f$ max = %s'%b[np.argmax(a)])
#plt.savefig('hist_tf.png', dpi=100)
b[np.argmax(a)]
Out[72]:
In [73]:
#Quantidade de pontos
a, b, c = plt.hist(dados.T[2], bins=140)
plt.xlim(0,1400)
#plt.title('Quantidade de pontos RRab, max = %s'%b[np.argmax(a)])
#plt.title('Histograma $n$, max = %s'%b[np.argmax(a)])
plt.title('Histograma - $n$', fontsize=18)
plt.xlabel('Quantidade de pontos', fontsize=14)
plt.ylabel('Frequência', fontsize=14)
b[np.argmax(a)]
#plt.savefig('hist_n.png', dpi=100)
Out[73]:
A partir dos histogramas, temos que a maior frequencia dos dados são para $t_i = 2152.5019$, $t_f = 4539.4593$ e $n = 351.6$
De acordo com o artigo, o tempo de observação é $t_i = HJD - 2450000$, assim, o valor correto para os dados são, $t_i = 247152.5019$, $t_f = 249539.4593$
Com estes dados, podemos calcular a amostragem da seguinte forma, $$ f_s= \frac{1}{dt}$$
$$dt= \frac{t_f - t_i}{n}$$$$f_s = \frac{n}{t_f - t_i}$$Ou seja, quanto maior a quantidade de pontos maior a amostragem. Mantendo os pontos final e inicial fixos, posso lidar com a quantidade de pontos e a amostragem de forma simples.
In [5]:
tf = 4539.4593
ti = 2152.5019
n = 351.6
dt = (tf - ti) / n
p_med = 0.576
fs = 1 / dt
print('dt = ', dt)
print('amostragem = ', fs)
Os dados podem ser gerados da seguinte forma, \begin{align} m(t) &= A_0 + \sum_i^3 A_n \sin \Big( \frac{2 n \pi t}{P} \Big) + \sigma \end{align}
em que, \begin{align} A_0 = 15 \\ A_1 = -0.5 \\ A_2 = 0.15 \\ A_3 = -0.05 \\ P = \text{Período, } \\ \sigma = \text{ruído} \end{align}
In [6]:
def curva_de_luz(tempo, periodo, ruido, sigma,
a0=15, a1=-0.5, a2=0.15, a3=-0.05):
w = 2*np.pi*tempo/periodo
mag = a0 + a1*np.sin(w) + a2*np.sin(2*w) + \
a3*np.sin(3*w) + ruido*sigma
return mag
#Gerar curva com a amostragem necessaria
a = 0.5
amos = a*fs
tau = np.arange(2152.5019, 4539.4593, 1/amos)
n = np.arange(50, 350, 100)
b = np.linspace(0.1, 1.0, 10)
n = size(tau)
sigma = np.random.normal(0.0, 1.0, tau.size)
m = curva_de_luz(tau, 0.576, 0, sigma)
#Gerar curva original com maior quantidade de pontos
ttau = np.arange(2152.5019, 4539.4593, dt/4)
nn = size(ttau)
ssigma = np.random.normal(0.0, 1.0, ttau.size)
curva_original = curva_de_luz(ttau, 0.576, 0, ssigma)
#xlim(2150, 2160)
#subplot(211)
plot(tau, m, 'ko', label=str(a)+'*fs')
plot(ttau, curva_original, 'b-', label='full data')
plt.legend()
xlim(2500,3000)
#subplot(212)
#plot(ttau,
Out[6]:
In [6]:
In [6]: