Antes de iniciar qualquer coisa, preciso te passar algumas informações...
Existem duas maneiras de efetuar autocorrelação, no seu exemplo pude observar que tentou efetuar autocorrelação no domínio da frequência, a equação correta para efetuar Autocorrelação no domínio da frequência é:
É importante saber que também é possível efetuar autocorrelação no domínio do tempo:
Este último funciona meio que na Força Bruta, um pseudo código para autocorrelação da equação acima seria:
AC=zeros( 4096 , 1 );
for k=1:4096/2,
sum = 0;
for n=1:4096/2,
sum = sum + (x(n)*x(n + k));
end
AC(k) = AC(k) + sum;
end
A autocorrelação vai retornar os pontos mais parecidos dentro do sinal, o retorno dessa função(tanto no domínio do tempo quanto do domínio da frequência) vai demonstrar picos indicando as posições onde o sinal possui maior semelhança, estes picos podem ser traduzidos como periodicidade/frequência, geralmente eu acho a diferença/subtraio a posição de um pico pelo outro para encontrar a periodicidade do sinal analisado.
De cara eu consigo saber a taxa de amostragem dos seus dados, pelos valores apresentados:
esta amostra apresenta cerca de 1600 dias, porém o arquivo tem na
coluna X cerca de 64.000 (sessenta e quatro MIL) pontos!
Isso quer dizer que você possui uma taxa de amostragem de 40
amostras por dia = 64000/1600=40
, ou seja em 1 dia é recolhido 40 amostras que compõem os seus dados, se você quiser saber a periodicidade de 30 dias, vai precisar de 30*40=1200
amostras, quando você utiliza fourier é preciso se preocupar com a ordem de resolução, dá uma lida nessa minha resposta.
Posso demonstrar como encontrar a periodicidade de um conjunto de dados usando Fourier, como eu não tenho os seus dados em mãos, para demonstrar vou criar uma senoide com frequência em 500Hz
e com uma taxa de amostragem em 44100hz
, isso me diz que estamos procurando uma periodicidade de 44100/500=88.2000
, guarde esse número é o que estamos procurando como resultado para essas minhas entradas de dados, para demonstrar segue o plot de uma senoide à 500hz com taxa de amostragem em 44100hz(essa taxa de amostragem diz que a cada 1 segundo eu estou coletando 44100 amostras), veja o plot:
É um sinal/dado simples, consigo a olho nú encontrar repetições, deste modo eu marquei o período onde a olho parece ocorrer a primeira repetição, olhá lá x=89
, lembra do valor que estamos perseguindo lá encima né ? 88.2000
, ou seja a olho nú eu consigo definir a periodicidade desses dados, mas e como fazer isso utilizando a primeira equação ?
Segue meu código:
import numpy as np
#inicio gerando o sinal mostrado no plot
Fs = 44100
freq = 500
nsamples = 4096
sinal = np.arange(nsamples)
sinal = np.sin(2 * np.pi * freq * sinal / Fs)
#fim, sinal criado
#aplicando primeira equação
Mag = np.abs(np.fft.rfft(sinal,nsamples*2))**2
AC=np.fft.irfft(Mag[0:nsamples])
#fim da equação
#encontrar onde estão os picos dentro da AutoCorrelação
peaks = []
k=3
while(k < nsamples - 1):
y1 = (AC[k - 1])
y2 = (AC[k])
y3 = (AC[k + 1])
if (y2 > y1 and y2 >= y3):
peaks.append(k)
k=k+1
periodo = np.mean(np.diff(peaks))
print("Periodo usando a diferenca de picos da Autocorrelacao")
print(periodo)
O Resultado desse algoritmo:
C:\Python33>python.exe AC.py
Periodo usando a diferenca de picos da Autocorrelacao
88.1555555556
Bingo 88.1555555556
praticamente o resultado esperado...
Existe uma pegadinha nessa parte do código np.fft.rfft(sinal,nsamples*2))
, para efetuar Autocorrelação e não Convolução(são parecidos, mas não iguais), você tem que computar o dobro de samples nsamples*2
isso aplica zeropad
na segunda metade do sinal, ou seja metade são os seus dados e a outra metade é composto de zeros...
Existe outra maneira de achar a Frequência/Período que é somente aplicar fourier
e observar onde existem picos, a Série de Fourier existe pra esse propósito a transformada de Fourier de um sinal periódico gera um espectro discreto no domínio da frequência
, somente utilizando Fourier:
import numpy as np
Fs = 44100
freq = 500
nsamples = 4096
sinal = np.arange(nsamples)
sinal = np.sin(2 * np.pi * freq * sinal / Fs)
Mag = np.abs(np.fft.rfft(sinal,nsamples*2))**2
index = Mag[0:nsamples/2].argmax()
print("Periodo usando o maior componente espectral de Fourier")
print(1 / (index / nsamples/2))
Resultado:
C:\Python33>python.exe AC.py
Periodo usando o maior componente espectral de Fourier
88.0860215054
Também cheguei bem próximo, tendo em vista que a ordem de resolução para esse exemplo é de 44100/ 4096 = 10,7666015625hz
Claro além desse dois métodos você ainda pode utilizar Cepstrum
ou autocorrelação no domínio do tempo da segunda equação...