In [13]:
import numpy
import tensorflow as tf
import pylab
maxFreq = 120
minFreq = 80
tFft = 8192
stepFreq = 1/tFft
nFreqs = numpy.int((maxFreq-minFreq)/stepFreq)
sessione = tf.Session()
detector1 = tf.random_normal([nFreqs], dtype=tf.float32)
detector2 = tf.random_normal([nFreqs], dtype=tf.float32)
spettroProdTF = tf.constant(spettroProd,dtype=tf.float32)
frequenze = tf.range(minFreq,maxFreq, delta = stepFreq)
alpha = tf.range(-4,4, delta = 0.1)
nAlphas = tf.size(alpha)
alpha = tf.reshape(alpha, shape=[nAlphas,1] )
alpha = 0
spettroSegnale = tf.pow(frequenze, alpha-3)
#detector1 = detector1+spettroSegnale*10
#detector2 = detector2+spettroSegnale*10
#detector1NP = sessione.run(detector1)
#detector2NP = sessione.run(detector2)
funzione = spettroSegnale/(spettroProdTF)
funzNP = sessione.run(funzione)
prodotto = detector2*funzione
estimatoreSum = detector1*prodotto
#righe = tf.shape(inputConv)[0]
#colonne = tf.shape(inputConv)[1]
inputConv = prodotto
inputConv = tf.reshape(inputConv,(1,1,nFreqs,1))
kernel = tf.reshape(detector1, (1,nFreqs,1,1))
#kernel = tf.reshape(detector1, (tf.size(detector1),1,1,1))
estimatore = tf.nn.convolution(inputConv,kernel, padding = 'VALID')
#swipe = [0,0,0,0]
#estimatore = tf.nn.conv2d(inputConv,kernel,strides=swipe, padding = 'SAME')
#inputConv = tf.reshape(detector1,(1,101,1))
#kernel = tf.reshape(detector2, (101,1,1))
#crosscorr = tf.nn.conv1d(inputConv,kernel,stride = 1,padding='VALID')/101
#cross = sessione.run(crosscorr)
print(estimatoreNP, autoestimatoreNP)
#inputConv = prodotto
#inputConv = tf.reshape(inputConv,(1,nFreqs,nAlphas,1))
#kernel = tf.reshape(detector1, (nFreqs,1,1,1))
#estimatore = tf.nn.convolution(inputConv,kernel, padding = 'VALID')
#estimatoreMulti = sessione.run(estimatore)
#print(estimatoreMulti.shape)
#alphaNP = numpy.arange(-4,4, 0.1)
freqNP = numpy.arange(80,120, 1/8192)
##estimatoreMulti = numpy.zeros(alphaNP.size)
##for i in numpy.arange(alphaNP.size-1):
##inputConv = prodotto[0]
##inputConv = tf.reshape(inputConv,(1,1,nFreqs,1))
##estimatore = tf.nn.convolution(inputConv,kernel, padding = 'VALID')
###start = time.time()
##estimatoreMulti[i] = sessione.run(estimatore)
###stop = time.time()
###print(stop-start)
#spettroSegnale = sessione.run(spettroSegnale)
##%matplotlib notebook
#a = pyplot.scatter(alphaNP,estimatoreMulti, s = 5)
#pyplot.show()
%time estimatore = sessione.run(estimatore)
%time estimatoreSum = sessione.run(estimatoreSum)
In [14]:
print(estimatore)
In [40]:
%matplotlib notebook
pylab.plot(freqNP, funzNP)
#pylab.plot(L, 1e-8*F*L/c)
#pylab.tight_layout()
#pylab.loglog()
pylab.show()
In [24]:
vettore1 = tf.constant([1,2,3,4])
vettore2 = tf.constant([[2],[2],[2]])
matrice = tf.multiply(vettore1,vettore2)
print(sessione.run(matrice))
In [6]:
import numpy
han = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2016-12-08_C01_H1_O2_Sensitivity_strain_asd.txt')
liv = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2017-07-20_C01_L1_O2_Sensitivity_strain_asd.txt')
maxFreq = 120
minFreq = 80
tFft = 8192
stepFreq = 1/tFft
nFreqs = numpy.int((maxFreq-minFreq)/stepFreq)
filtrohan = numpy.where(han[:,0] > 80)
filtroliv = numpy.where(liv[:,0] > 80)
hansel = han[filtrohan]
livsel = liv[filtroliv]
filtrohan = numpy.where(hansel[:,0] < 120)
filtroliv = numpy.where(livsel[:,0] < 120)
hansel = hansel[filtrohan]
livsel = livsel[filtroliv]
spettroCoord = hansel[:,0]
spettroHan = hansel[:,1]
spettroLiv = livsel[:,1]
spettroHanExp = numpy.zeros(nFreqs)
spettroLivExp = numpy.zeros(nFreqs)
freqNP = numpy.arange(minFreq,maxFreq,stepFreq)
print(freqNP,spettroCoord.size)
filtro = numpy.where(freqNP<=spettroCoord[0])[0]
fine = filtro[filtro.size-1]
spettroHanExp[:fine] = spettroHan[0]
spettroLivExp[:fine] = spettroLiv[0]
for i in numpy.arange(spettroCoord.size-1):
filtro1 = numpy.where(freqNP<=spettroCoord[i+1])[0]
filtro2 = numpy.where(freqNP>spettroCoord[i])[0]
inizio = filtro2[0]-1
fine = filtro1[filtro1.size-1]
#print(inizio,fine)
spettroHanExp[inizio:fine] = spettroHan[i]
spettroLivExp[inizio:fine] = spettroLiv[i]
filtro = numpy.where(freqNP>spettroCoord[-1])[0]
inizio = filtro[0]-1
spettroHanExp[inizio:] = spettroHan[-1]
spettroLivExp[inizio:] = spettroLiv[-1]
print(spettroHanExp, spettroLivExp)
spettroHanExp = spettroHanExp*1e23
spettroLivExp = spettroLivExp*1e23
print(spettroHanExp, spettroLivExp)
spettroProd = numpy.multiply(spettroHanExp,spettroLivExp)
print(spettroProd)
In [7]:
import numpy
import pylab
han = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2016-12-08_C01_H1_O2_Sensitivity_strain_asd.txt')
liv = numpy.loadtxt('/home/protoss/Documenti/TESI/DATI/strainsensib/2017-07-20_C01_L1_O2_Sensitivity_strain_asd.txt')
%matplotlib qt
pylab.figure(figsize=(10,7.5))
pylab.loglog()
pylab.plot(han[:,0],han[:,1],label = "LIGO H",linewidth=1.0)
pylab.plot(liv[:,0],liv[:,1],label = "LIGO L",linewidth=1.0)
pylab.xlabel('Frequencies ($Hz$)')
pylab.ylabel('Strain amplitude ($1/\sqrt{Hz}$)')
pylab.title('Noise spectral amplitude')
pylab.legend()
Out[7]:
In [19]:
from matplotlib import pyplot
%matplotlib notebook
pyplot.figure(figsize=(10, 8))
#pyplot.loglog()
a = pyplot.scatter(freqNP,spettroHanExp, s = 0.5)
a = pyplot.scatter(freqNP,spettroLivExp, s = 0.5)
a = pyplot.scatter(freqNP,spettroProd, s = 0.5)
In [41]:
from matplotlib import pyplot
%matplotlib notebook
pyplot.figure(figsize=(10, 8))
pyplot.loglog()
a = pyplot.scatter(han[:,0],han[:,1], s = 0.5)
b = pyplot.scatter(liv[:,0],liv[:,1], s = 0.5)
In [45]:
In [62]:
print(livsel[:,0]- hansel[:,0])
In [90]:
from matplotlib import pyplot
%matplotlib notebook
pyplot.figure(figsize=(10, 8))
pyplot.loglog()
a = pyplot.scatter(hansel[:,0],hansel[:,1], s = 0.5)
b = pyplot.scatter(livsel[:,0],livsel[:,1], s = 0.5)
In [ ]:
In [47]:
import numpy
import tensorflow as tf
from matplotlib import pyplot
import time
maxFreq = 120
minFreq = 80
tFft = 8192
stepFreq = 1/tFft
nFreqs = numpy.int((maxFreq-minFreq)/stepFreq)
sessione = tf.Session()
detector1 = tf.random_normal([nFreqs], dtype=tf.float32)
detector2 = tf.random_normal([nFreqs], dtype=tf.float32)
#spettroProdTF = tf.constant(spettroProd,dtype=tf.float32)
frequenze = tf.range(minFreq,maxFreq, delta = stepFreq)
alpha = tf.range(0,3, delta = 0.01)
nAlphas = tf.size(alpha)
alpha = tf.reshape(alpha, shape=[nAlphas,1] )
espansoreFreq = tf.ones([nAlphas,1])
espansoreAlpha = tf.ones([1,nFreqs])
freqExp = tf.multiply(frequenze,espansoreFreq)
alphExp = tf.multiply(espansoreAlpha,alpha)
spettroSegnale = tf.pow(freqExp, alphExp-3)
#detector2 = detector2/spettroProdTF
prodotto = detector2*spettroSegnale
inputConv = prodotto#[0]
inputConv = tf.reshape(inputConv,(1,nAlphas,nFreqs,1))
kernel = tf.reshape(detector1, (1,nFreqs,1,1))
estimatore = tf.nn.convolution(inputConv,kernel, padding = 'VALID')
start = time.time()
estimatoreMulti = sessione.run(estimatore)
stop = time.time()
print(stop-start)
print(estimatoreMulti[0,:,0,0])
alphaNP = numpy.arange(0,3, 0.01)
%matplotlib notebook
#pyplot.loglog()
a = pyplot.scatter(alphaNP,estimatoreMulti[0,:,0,0], s = 5)
pyplot.show()
In [ ]:
In [30]:
inputConv = detector2
inputConv = tf.reshape(inputConv,(1,1,nFreqs,1))
kernel = funzione[0]
kernel = tf.reshape(kernel, (1,nFreqs,1,1))
matchingfilter = tf.nn.convolution(inputConv,kernel, padding = 'SAME')
inputConv = detec1
inputConv = tf.reshape(inputConv,(1,1,nFreqs,1))
kernel = matchingfilter
kernel = tf.reshape(kernel, (1,nFreqs,1,1))
estimatore = tf.nn.convolution(inputConv,kernel, padding = 'SAME')
In [67]:
import numpy
import tensorflow as tf
from matplotlib import pyplot
import time
maxFreq = 120
minFreq = 80
tFft = 8192
stepFreq = 1/tFft
nFreqs = numpy.int((maxFreq-minFreq)/stepFreq)
#nFreqs = 101
sessione = tf.Session()
def preparDati(N,sign):
valori = tf.random_normal([N], dtype=tf.float32)
ampiezza = tf.pow(tf.abs(valori),2)
dati = ampiezza + sign
return dati
frequenze = tf.range(minFreq,maxFreq, delta = stepFreq)
alpha = tf.range(-4,4, delta = 0.1)
nAlphas = tf.size(alpha)
alpha = tf.reshape(alpha, shape=[nAlphas,1] )
alpha = 3
refreq = minFreq
spettroSegnale = tf.pow(frequenze/refreq, alpha)
detector1 = preparDati(nFreqs,spettroSegnale)
detector2 = preparDati(nFreqs,spettroSegnale)
#detector1 = preparDati(nFreqs,0)
#detector2 = preparDati(nFreqs,0)
def crosscorrela(funz1,funz2, N):
media1 = tf.reduce_mean(funz1)
media2 = tf.reduce_mean(funz2)
funz1= funz1-media1
funz2 = funz2-media2
sigma1 = tf.sqrt(tf.reduce_sum(funz1*funz1)/N)
sigma2 = tf.sqrt(tf.reduce_sum(funz2*funz2)/N)
prodotto = funz1*funz2
crosscorr = (tf.reduce_sum(prodotto)/N)/(sigma1*sigma2)
return crosscorr
def matchedfiltra(funz,filtro, N):
media = tf.reduce_mean(funz)
funz= funz-media
filtro = filtro-media
sigma = tf.sqrt(tf.reduce_sum(funz*funz)/N)
prodotto = filtro*funz/tf.pow(sigma,2)
normalizzatore = tf.reduce_sum(tf.pow(filtro/sigma,2))
filtraggio = tf.reduce_sum(prodotto)/normalizzatore
return filtraggio
#def crossfiltra(funz1,funz2,sign,N)
filtrato1 = matchedfiltra(detector1,spettroSegnale,nFreqs)
filtrato2 = matchedfiltra(detector2,spettroSegnale,nFreqs)
estimatore = crosscorrela(detector1, detector2, nFreqs)
autoestimatore = crosscorrela(detector1,detector1, nFreqs)
fil1NP = sessione.run(filtrato1)#detectorsNP[:,2]
fil2NP = sessione.run(filtrato2)#detectorsNP[:,3]
print(fil1NP,fil2NP)
detectorsNP = sessione.run(tf.stack([detector1,detector2],1))
det1NP = detectorsNP[:,0]
det2NP = detectorsNP[:,1]
estimatoreNP = sessione.run(estimatore)
autoestimatoreNP = sessione.run(autoestimatore)
In [72]:
import numpy
from matplotlib import pyplot
import tensorflow as tf
spettroHan = numpy.load('/home/protoss/Documenti/Gravit_sper/dati/spettroHanMediato.npy')
spettroLiv = numpy.load('/home/protoss/Documenti/Gravit_sper/dati/spettroLivMediato.npy')
han = tf.constant(spettroHan)+spettroSegnale
liv = tf.constant(spettroLiv)+spettroSegnale
crossSpettri = sessione.run(crosscorrela(han,liv, nFreqs))
autoCorrHan = sessione.run(crosscorrela(han,han,nFreqs))
autoCorrLiv = sessione.run(crosscorrela(liv,liv,nFreqs))
filtroHan = sessione.run(matchedfiltra(han, spettroSegnale,nFreqs))
filtroLiv = sessione.run(matchedfiltra(liv, spettroSegnale,nFreqs))
print(crossSpettri, autoCorrHan, autoCorrLiv,filtroHan,filtroLiv)
freqNP = numpy.arange(80,120, 1/8192)
%matplotlib notebook
pyplot.loglog()
a = pyplot.scatter(freqNP,spettroHan, s = 5)
a = pyplot.scatter(freqNP,spettroLiv, s = 5)
#pyplot.show()
In [52]:
spettroHan
Out[52]:
In [ ]:
In [48]:
import numpy
import tensorflow as tf
from matplotlib import pyplot
import time
maxFreq = 120
minFreq = 80
samplerate = 16384
tFft = 192
sizetSeries = 16384*tFft
stepFreq = 1/tFft
nFreqs = numpy.int((maxFreq-minFreq)/stepFreq)
sessione = tf.Session()
def preparDati(N,sign=0):
valori = tf.random_normal([N], dtype=tf.float32)
mediaAmp = tf.reduce_mean(valori)
dati = valori - mediaAmp
dati = dati + sign
return dati
timeseries1 = preparDati(sizetSeries)
#timeseries1 = tf.cast(timeseries1, dtype=tf.complex64)
timeseries2 = preparDati(sizetSeries)
#timeseries2 = tf.cast(timeseries2, dtype=tf.complex64)
#che succede se crosscorrelo sui tempi e faccio trasf di fourier
#invece di fare media delle componenti di fourier come maggiore & romano?
#fft1 = tf.fft(timeseries1)
#fft2 = tf.fft(timeseries2)
#sizefft = tf.size(fft1)
times1NP = sessione.run(timeseries1)
times2NP = sessione.run(timeseries2)
#fftNP = sessione.run(fft1)
fft1NP = numpy.fft.fft(times1NP)
fft2NP = numpy.fft.fft(times2NP)
fft1 = tf.constant(fft1NP, dtype=tf.complex64)
fft2 = tf.constant(fft2NP, dtype=tf.complex64)
prod = tf.conj(fft1)*fft2
print(sessione.run(prod))
def crossfiltra(dati1,dati2,sign=1):
prodComplesso = tf.conj(dati1)*dati2*sign
estimatore = tf.reduce_sum(prodComplesso)
return estimatore
estimatore = crossfiltra(fft1,fft2)
print(sessione.run(estimatore))
sizefft = tf.size(fft1)
frequenze = tf.range(sizefft)
alpha = tf.range(-4,4, delta = 0.1)
nAlphas = tf.size(alpha)
alpha = tf.reshape(alpha, shape=[nAlphas,1] )
alpha = 3
refreq = minFreq
spettroSegnale = tf.pow(frequenze/refreq, alpha)
In [43]:
#non ho capito come funziona ora il range di frequenza. che intervallo di frequenze copro con questo tfft
#e questo sampling rate?
#da 0 a 192 hz con 16384 valori per ogni hertz? ma se sensib dovrebbe essere 1/192, dovrei mediare i valori
#in mezzo?
fft1NP.size/192
Out[43]:
In [26]:
a = numpy.array([1,2,3,4])
print(a)
numpy.flip(a,0)
Out[26]:
In [11]:
%matplotlib notebook
a = pyplot.scatter(numpy.arange(fft1NP.size),fft1NP, s = 5)
pyplot.imshow(a)
pyplot.show()
In [ ]: