In [5]:
import glob
import csv
from operator import add
import scipy.stats as stat
#Dateiliste erstellen
files = []
for filename in glob.glob('/Users/robinweiss/Documents/RWeiss/final data/*.csv'):
files.append(filename)
#read files
seen = 0
trials = 0
catchseen = 0
catchtrials = 0
deltas = []
freqs = []
for filename in files:
data = csv.reader(open(filename, "rb"), delimiter=',')
data.next()
data.next()
data.next()
for row in data:
freqdiff = abs(int(row[0][0:-2])-int(row[1][0:-2]))
if float(row[3]) < float(row[4]):
if row[0] != row[1]:
trials += 1
if row[5] == 'y':
seen += 1
deltas.append((freqdiff,1))
freqs.append((int(row[0][0:-2]),int(row[1][0:-2]),1))
else:
deltas.append((freqdiff,0))
freqs.append((int(row[0][0:-2]),int(row[1][0:-2]),0))
if row[0] == row[1]:
catchtrials += 1
if row[5] == 'n':
catchseen += 1
else:
print row
d_delta = {}
for freq, count in deltas:
if freq not in d_delta:
d_delta[freq] = [0,0]
d_delta[freq][0] += 1
if count:
d_delta[freq][1] += 1
d_freqs = {}
for freq1, freq2, count in freqs:
if (freq1, freq2) not in d_freqs:
d_freqs[(freq1, freq2)] = [0,0]
d_freqs[(freq1, freq2)][0] += 1
if count:
d_freqs[(freq1, freq2)][1] += 1
d_all = d_freqs.copy()
d_all['overall\ncorrectness'] = (1, round(float(seen)/trials, 3))
d_all['false\npositives'] = (1, 1-round(float(catchseen)/catchtrials, 3))
delta_std = loadtxt("/Users/robinweiss/Documents/RWeiss/final data/delta_std.txt", delimiter=', ')
d_all_std = loadtxt("/Users/robinweiss/Documents/RWeiss/final data/d_all_std.txt", delimiter=', ')
contrast_std = loadtxt("/Users/robinweiss/Documents/RWeiss/final data/contrast_std.txt", delimiter=', ')
delta_std = [stat.sem(delta_std[:,i]) for i in range(delta_std.shape[1])]
d_all_std = [stat.sem(d_all_std[:,i]) for i in range(d_all_std.shape[1])]
contrast_std = [stat.sem(contrast_std[:,i]) for i in range(contrast_std.shape[1])]
#Simulation
def rectwindow(freq1, freq2, f_sample=10, width=1, sig=.1, filtering=0): #f_sample in kHz und width ist der Anteil des Fensters am Gesamtsignal (1 heisst, kein Fenster)
"Erstellt einen Rechteckfenster-Frequenzübergang freq1 auf freq2 mit Signallänge 2*sig (in s) und entsprechende FFT"
freq1, freq2, width = float(freq1), float(freq2), float(width)
f_sample = f_sample*1000.
f1 = tile(concatenate((ones(f_sample/(2*freq1)), zeros(f_sample/(2*freq1))), axis=1), sig*freq1) #erste Frequenz
f2 = tile(concatenate((ones(f_sample/(2*freq2)), zeros(f_sample/(2*freq2))), axis=1), sig*freq2) #zweite Frequenz
f = concatenate((f1, f2), axis=1)-0.5
window = concatenate((zeros(((width-1)/(2*width))*f.size), ones(f.size/width)), axis=1)
window = concatenate((window, zeros(((width-1)/(2*width))*f.size)), axis=1)
window = concatenate((window, zeros(f.size-window.size)), axis=1)
f = f*window
xaxis = linspace(0, 2*sig, f.size)
FTdata = fft.fft(f)/f.size
FT = FTdata[:(FTdata.size/2)]
if filtering:
u = arange(0,.12,.001)
visfilter = .005*(((u**3)/(.005**4))*exp(-u/.005)-((u**3)/(.01**4))*exp(-u/.01))
FTfilter = abs(fft.fft(visfilter, 1000)/21.6)[:300]
FTfilter = concatenate((FTfilter, zeros(len(FT)-len(FTfilter))), axis=1)
FT = FT*FTfilter
faxis = fft.fftfreq(f.size, (1./f_sample))
faxis = faxis[:(faxis.size/2)]
integral = sum(abs(FT)[:600])
#print integral
#plot(abs(FT)[0:100])
return FTdata, xaxis, faxis, integral;
#rectwindow(200,600,168.,1,3.)
#rectwindow(200,600,168.,1,3.,1)
integrals = []
for item in sorted(d_freqs):
hin = rectwindow(item[0],item[1],168,1.,3.)
rueck = rectwindow(item[1],item[0],168,1.,3.)
integrals.append((hin[3]+rueck[3]))
filterintegrals = []
for item in sorted(d_freqs):
hin = rectwindow(item[0],item[1],168,1.,3.,1)
rueck = rectwindow(item[1],item[0],168,1.,3.,1)
filterintegrals.append((hin[3]+rueck[3]))
In [6]:
figure(figsize=(7,5))
plot(list(integrals), [float(d_freqs[key][1])/d_freqs[key][0] for key in sorted(d_freqs)], 'o')
#plot([0,.12],[0,1], alpha=.5)
plt.xlabel('Energie', fontsize=15)
plt.ylabel('P', fontsize=15)
plt.yticks(fontsize=12)
plt.xticks(fontsize=12)
plt.ylim(-.1,1.1)
plt.title('x: Integral ueber das Spektrum 0-100Hz\ny: Wahrnehmungswahrscheinlichkeit ueber 6VPs gemittelt')
plt.tight_layout()
#plt.savefig('/home/weiss/ba/figures/scatter_integrals_experiment.pdf')
figure(figsize=(7,5))
plot(list(filterintegrals), [float(d_freqs[key][1])/d_freqs[key][0] for key in sorted(d_freqs)], 'o')
#plot([0,.12],[0,1], alpha=.5)
plt.xlabel('Energie', fontsize=15)
plt.ylabel('P', fontsize=15)
plt.yticks(fontsize=12)
plt.xticks(fontsize=12)
plt.ylim(-.1,1.1)
plt.title('x: Integral ueber das Spektrum 0-100Hz mit Filter gewichtet\ny: Wahrnehmungswahrscheinlichkeit ueber 6VPs gemittelt')
plt.tight_layout()
#plt.savefig('/home/weiss/ba/figures/scatter_filterintegrals_experiment.pdf')
In [7]:
#iFFT
def rectwindownn(freq1, freq2, f_sample=10, width=1, sig=.1): #f_sample in kHz und width ist der Anteil des Fensters am Gesamtsignal (1 heisst, kein Fenster)
"Erstellt einen Rechteckfenster-Frequenzübergang freq1 auf freq2 mit Signallänge 2*sig (in s) und entsprechende FFT"
freq1, freq2, width = float(freq1), float(freq2), float(width)
f_sample = f_sample*1000.
f1 = tile(concatenate((ones(f_sample/(2*freq1)), zeros(f_sample/(2*freq1))), axis=1), sig*freq1) #erste Frequenz
f2 = tile(concatenate((ones(f_sample/(2*freq2)), zeros(f_sample/(2*freq2))), axis=1), sig*freq2) #zweite Frequenz
f = concatenate((f1, f2), axis=1)-0.5
window = concatenate((zeros(((width-1)/(2*width))*f.size), ones(f.size/width)), axis=1)
window = concatenate((window, zeros(((width-1)/(2*width))*f.size)), axis=1)
window = concatenate((window, zeros(f.size-window.size)), axis=1)
f = f*window
xaxis = linspace(0, 2*sig, f.size)
FTdata = fft.fft(f)#/f.size
FT = FTdata[:(FTdata.size/2)]
faxis = fft.fftfreq(f.size, (1./f_sample))
faxis = faxis[:(faxis.size/2)]
return FTdata, xaxis, faxis;
u = arange(0,.12,.001)
visfilter = .005*(((u**3)/(.005**4))*exp(-u/.005)-((u**3)/(.01**4))*exp(-u/.01))
FTfilter = abs(fft.fft(visfilter, 1000)/21.6)[:300]
peakheights = []
for item in sorted(d_freqs):
FT, t, f = rectwindownn(item[0],item[1],168,1.,5.)
FTfilter = concatenate((FTfilter, zeros(len(FT)-len(FTfilter))), axis=1)
signal = fft.ifft(abs(FT)*FTfilter)
#print amax(array(abs(signal[200000:1300000])))
#print abs(signal[len(signal)/2])
peakheights.append(abs(signal[len(signal)/2]))
figure(figsize=(7,5))
plot(list(peakheights), [float(d_freqs[key][1])/d_freqs[key][0] for key in sorted(d_freqs)], 'o')
#plot([0,.12],[0,1], alpha=.5)
plt.xlabel('Hoehe der Peaks', fontsize=15)
plt.ylabel('P', fontsize=15)
plt.yticks(fontsize=12)
plt.xticks(fontsize=12)
plt.ylim(-.1,1.1)
plt.title('x: Hoehe des Peaks beim Filteroutput und Ruecktransformation\ny: Wahrnehmungswahrscheinlichkeit ueber 6VPs gemittelt')
plt.tight_layout()
#plt.savefig('/home/weiss/ba/figures/scatter_ifftpeaks_experiment.pdf')
In [22]:
# x-Werte gegeneinander
figure(figsize=(7,5))
plot(list(integrals), list(peakheights), 'o')
plot((0,0.12), (0,.00146), linewidth=.5)
plt.xlabel('Integrale (x-Werte)', fontsize=12)
plt.ylabel('Peaks (x-Werte)', fontsize=12)
plt.tight_layout()
#plt.savefig('/Users/robinweiss/Documents/RWeiss/figures/x-werte.pdf')
In [3]: