In [5]:
def rectwindow(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.
#print "Samplingfrequenz:", f_sample/1000, "kHz"
f1 = tile(concatenate((ones(f_sample/(2*freq1)), zeros(f_sample/(2*freq1))), axis=0), sig*freq1) #erste Frequenz
f2 = tile(concatenate((ones(f_sample/(2*freq2)), zeros(f_sample/(2*freq2))), axis=0), sig*freq2) #zweite Frequenz
f = concatenate((f1, f2), axis=0)-0.5
window = concatenate((zeros(((width-1)/(2*width))*f.size), ones(f.size/width)), axis=0)
window = concatenate((window, zeros(((width-1)/(2*width))*f.size)), axis=0)
window = concatenate((window, zeros(f.size-window.size)), axis=0)
f = f*window
xaxis = linspace(0, 2*sig, f.size)
#xaxis = arange(0, 2.*sig, 1./f.size)
#figure(figsize=(10,3))
#plot(xaxis, f)
#plt.xlabel('t [s]', fontsize=15)
#plt.ylabel('Stimulus', fontsize=15)
#plt.xticks([0, 0.05, 0.1], [0,1,2], fontsize=15)
#plt.yticks([-0.5, 0.5], ['aus', 'an'], fontsize=15)
#plt.ylim(-.7,.7)
#plt.xlim(0,0.1)
#plt.tight_layout()
#plt.savefig("/Users/robinweiss/Documents/RWeiss/figures/signalsim.pdf")
#figure(figsize=(10,2))
FTdata = fft.fft(f)#/f.size
FT = FTdata[:(FTdata.size/2)]
faxis = fft.fftfreq(f.size, (1./f_sample))
faxis = faxis[:(faxis.size/2)]
#plot(faxis[0:600], abs(FT)[0:600], linewidth=2)
#plt.xlabel('Frequenz [Hz]', fontsize=15)
#plt.xlim(0,40)
#plt.ylim(0,.0005)
#plt.ylabel('Amplitude', fontsize=15)
#plt.yticks([0, 0.0003], fontsize=15)
#plt.xticks(fontsize=15)
#plt.tight_layout()
#plt.savefig('/Users/robinweiss/Documents/RWeiss/figures/zacken_200_400.pdf')
#print abs(FT[:101])
integral = sum(abs(FT)[:600])
#print integral
#figure()
#plt.psd(f, f.size/2, f_sample)
#plt.xlim(0, 800)
return FTdata, xaxis, faxis, integral;
'''
rectwindow(128, 128, 16.384, 10, 5.) #2er-Potenzen mit 1Hz pro bin
rectwindow(128, 512, 16.384, 10, 5.)
rectwindow(200, 400, 24, 1., .05)
#rectwindow(300, 400, 24, 10., 5.)
integrals200 = []
integrals400 = []
integrals600 = []
integrals800 = []
freqs = [200,300,400,500,600,700,800]
for i in freqs:
x = rectwindow(200,i,168,1.,3.)
integrals200.append(x[3])
for i in freqs:
x = rectwindow(400,i,168,1.,3.)
integrals400.append(x[3])
for i in freqs:
x = rectwindow(600,i,168,1.,3.)
integrals600.append(x[3])
for i in freqs:
x = rectwindow(800,i,168,1.,3.)
integrals800.append(x[3])
print integrals200
x = arange(7)
figure(figsize=(10,6))
plot(x, integrals200, marker='o', markersize=8, linewidth=2, label='200 Hz')
plot(x, integrals400, marker='o', markersize=8, linewidth=2, label='400 Hz')
plot(x, integrals600, marker='o', markersize=8, linewidth=2, label='600 Hz')
plot(x, integrals800, marker='o', markersize=8, linewidth=2, label='800 Hz')
plt.xticks(arange(7), freqs, fontsize=15)
plt.xlim(-.2, 6.2)
plt.yticks([0,.02,.04], fontsize=15)
plt.xlabel('Frequenz [Hz]', fontsize=15)
plt.ylabel('Amplitude', fontsize=15)
plt.legend(loc='best', fontsize=15)
plt.tight_layout()
#plt.savefig("/Users/robinweiss/Documents/RWeiss/figures/integrals.pdf")
'''
# Was passiert bei der iFFT? Mit lowpass und bandpass
# oben in der Funktion Normierung auskommentieren!! wegen Rücktransformation: Asymmetrie siehe fft doc
FT1, t, f, i= rectwindow(200,400,24.,1.,5.)
FT2, t, f, i = rectwindow(200,600,24.,1.,5.)
FT3, t, f, i = rectwindow(200,800,24.,1.,5.)
x = arange(0,.12,.001)
y = .005*(((x**3)/(.005**4))*exp(-x/.005)-((x**3)/(.01**4))*exp(-x/.01))
Y = abs(fft.fft(y, 1000)/21.6)[:300]
Y = concatenate((Y, zeros(FT1.size-300)))
#signal = fft.ifft(abs(FT1)*Y)
figure(figsize=(10,6))
plot(t, fft.ifft(abs(FT1)*Y), linewidth=2, label='200Hz, 400Hz')
plot(t, fft.ifft(abs(FT2)*Y), linewidth=2, label='200Hz, 600Hz')
plot(t, fft.ifft(abs(FT3)*Y), linewidth=2, label='200Hz, 800Hz')
plt.xticks(range(10), fontsize=15)
plt.yticks(fontsize=15)
plt.xlim(4,6)
plt.ylim(-.0015, .0003)
plt.grid()
plt.legend(loc='lower right', fontsize=15)
plt.xlabel('t [s]', fontsize=15)
plt.ylabel('Amplitude', fontsize=15)
plt.tight_layout()
#plt.savefig("/Users/robinweiss/Documents/RWeiss/figures/ifft.pdf")
'''
#lowpass = zeros((1,504000))
#lowpass[0,0:680] = c[0:680]
#lowpass[0,502000:] = c[502000:]
sig = fft.ifft(lowpass*c.shape)
figure()
plot(x, real(sig[0,:]))
plt.title('Lowpass-iFFT Flicker 200->700Hz')
plt.xlabel('t in s')
plt.ylabel('Amplitude')
plt.xlim(5.7,6)
'''
'''
figure()
bandpass = concatenate((signal.tukey(800, sym=False), zeros(501081)), axis=1)
bandpass = concatenate((bandpass, signal.tukey(800, sym=False)), axis=1)
bandpass = abs(c)*bandpass
sig = fft.ifft(bandpass*c.shape)
figure()
plot(x, real(sig))
plt.title('bandpass-iFFT Flicker 200->700Hz')
plt.xlabel('t in s')
plt.xlim(5,7)
#plt.ylim(-0.005, 0.005)
figure()
autocorr = correlate(sig[0,240000:260000], sig[0,240000:260000], mode='full')
u = arange(0, autocorr.size)
plot(u[20000:24000], autocorr[20000:24000])
plt.title('Autocorrelation')
figure()
plot(f[:1500], abs(a[:1500]), label=('200Hz'))
plt.title('200Hz FFT')
figure()
#plot(f[:500], abs(b[:500]), label=('700Hz'))
plot(f[:500], real(c[:500]), label=('real'))
plot(f[:500], imag(c[:500]), label=('imag'))
plot(f[:500], abs(c[:500]), label=('abs'))
pylab.legend(loc='upper right')
plt.title('FFT-Spektren im niedrigen Frequenzbereich')
plt.xlabel('Frequenz in Hz')
plt.ylabel('Amplitude')
'''
Out[5]:
In [ ]: