In [1]:

    

from matplotlib.mlab import psd

def rectwindow(freq1, freq2, f_sample=10, width=3, N=100):
    "Erstellt einen rect windowed Frequenzübergang Freq1 auf Freq2 mit N Perioden Freq1 Rechteckfunktion+selbe Anz Data Freq2 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=1), N)
    f2 = tile(concatenate((ones(f_sample/(2*freq2)), zeros(f_sample/(2*freq2))), axis=1), round(N*20*freq2/(20*freq1))+1)
    f2 = f2[:(f2.size-(f2.size-f1.size))]
    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
    window_neg = window*(-1)
    f = f*window_neg
    #xaxis = linspace(0, (N/freq1)+(N/freq2), f.size)
    #plot(xaxis, f)
    #plt.xlabel('t in s')
    #plt.ylim(-1,1)
    #figure()
    #FT = fft.fft(f)/f.size
    #FT = FT[:(FT.size/2)]
    #faxis = fft.fftfreq(f.size, (1./f_sample))
    #faxis = faxis[:(faxis.size/2)]
    #plot(faxis, abs(FT))
    #plt.xlim(0, (freq1+freq2))
    #plt.xlabel('Freq in Hz')
    #plt.ylabel('Amplitude')
    #figure()
    a, freqs = psd(f, 12000, f_sample, noverlap=6000)
    #plt.xlim(0, (freq1+freq2))
    return 20*log10(a), freqs;

def rectwindow3d(freq1, freq2, f_sample=10, data=100, step=0.1):
    FT, faxis = rectwindow(freq1, freq2, f_sample)
    print "Samplingfrequenz:", f_sample, "kHz"
    S = zeros((data-10, FT.size))
    l = arange(1, data*step, step)
    for n, item in enumerate(l):
        sn, faxis = rectwindow(freq1, freq2, f_sample, item)
        S[n,:] = sn.T
    return S, faxis, l;

S, faxis, l = rectwindow3d(300,500,15,200)

figure()
im = imshow(S[:,:1000], extent=[0, faxis[1000], 9, 1], aspect=('auto'))
plt.colorbar()


Smax = amax(S, axis=0)
S_norm = zeros(S.shape)
for i, item in enumerate(Smax):
    S_norm[:,i] = -item/abs(S[:,i])


figure()
im = imshow(S_norm[:,:1000], extent=[0, faxis[1000], 9, 0], aspect=('auto'), cmap='hot')
plt.colorbar()


    

    




Samplingfrequenz: 15 kHz






    
Out[1]:





<matplotlib.colorbar.Colorbar instance at 0x7f07683f0c68>






    













    

In [5]:

    

from mpl_toolkits.mplot3d import Axes3D

x, y = meshgrid(faxis[:1000], l)
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
ax.plot_surface(x, y, S[:,:1000], cmap=plt.get_cmap("autumn"))
plt.show()


    

In [ ]: