3

Olá, pessoal!

Eu tenho um arquivo .txt de dados com duas colunas (em X tenho os valores em dias de observações e em Y o Fluxo Medido de uma amostra)

Pretendo calcular o Período com que esses dados se repetem e desejaria fazer isso usando a ACF (autocorrelation function) no Python, mas até agora não conseguir.

O que tenho feito é o seguinte:

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

dado = np.loadtxt(r'.../file.txt') #busco meu arquivo no meu PC

t = dado[:,0] #esta é apenas a coluna que corresponde ao tempo
y = dado[:,1] #esta é apenas a coluna que corresponde ao fluxo medido

# Fazer uma filtragem a partir da Transformada Rápida de Fourier (Fast Fourier Transformed - FFT)

fft_of_signal_with_noise = np.fft.fft(y)
# Normalizar o tempo
N = len(t)
f = np.fft.fftfreq(len(fft_of_signal_with_noise),1600)
frequency_of_signal  = 1.0 


def bandpass_filter(x, freq, frequency_of_signal=frequency_of_signal, band = 0.05):
    if (frequency_of_signal - band) < abs(freq) < (frequency_of_signal + band):
        return x
    else:
        return 0

F_filtered = np.asanyarray([bandpass_filter(x,freq) for x,freq in zip(fft_of_signal_with_noise, f)]);

filtered_signal = np.fft.ifft(F_filtered);

Porém ainda não obtive sucesso em encontrar o valor correto do Período desta amostra, por exemplo.

OBS: esta amostra apresenta cerca de 1600 dias, porém o arquivo tem na coluna X cerca de 64.000 (sessenta e quatro MIL) pontos!

Alguém pode ajudar-me a encontrar o valor do Período com que se reperte?

3
  • humm mas esse seu código não implementa autocorrelação, já conseguiu resolver o seu problema? se ainda não conseguiu eu posso tentar explicar como efetuar os cálculos ...
    – ederwander
    Commented 2/05/2018 às 19:12
  • Ainda não obtive sucesso, mas adoraria ver como funciona. Commented 21/05/2018 às 18:56
  • Espero que a minha resposta te ajuda
    – ederwander
    Commented 23/05/2018 às 13:20

1 Resposta 1

1

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 FFT.py
Periodo usando o maior componente espectral de Fourier
88.0860215054

Sem precisar de AutoCorrelação cheguei bem próximo, tendo em vista que a ordem de resolução para esse exemplo é de 44100/ 4096 = 10,7666015625hz

Claro além desses dois métodos você ainda pode utilizar Cepstrum ou autocorrelação no domínio do tempo da segunda equação...

Você deve fazer log-in para responder a esta pergunta.

Esta não é a resposta que você está procurando? Pesquise outras perguntas com a tag .