Skip to main content
1 de 3
ederwander
  • 6,6mil
  • 24
  • 32

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 é:

inserir a descrição da imagem aqui

É importante saber que também é possível efetuar autocorrelação no domínio do tempo:

inserir a descrição da imagem aqui

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:

inserir a descrição da imagem aqui

É 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...

ederwander
  • 6,6mil
  • 24
  • 32