In [1]:

    

from matplotlib.mlab import psd

def hannwindow(freq1, freq2, f_sample=10, width=3, sig=.5):
    "Erstellt einen hamming 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), sig*freq1)
    f2 = tile(concatenate((ones(f_sample/(2*freq2)), zeros(f_sample/(2*freq2))), axis=1), sig*freq2)
    f = concatenate((f1, f2), axis=1)-0.5
    xaxis = linspace(0, 2*sig, f.size)
    window = concatenate((zeros(((width-1)/(2*width))*f.size), hanning(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
    #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, 29929, f_sample)
    #plt.xlim(0, (freq1+freq2))
    #y = concatenate((S, abs(FT)))
    return 10*log10(a), freqs;


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


figure()
im = imshow(S[:,:1000], extent=[0, faxis[1000], 20, 1], aspect=('auto'))
plt.xlabel('Frequenz in Hz')
plt.ylabel('Fensterbreite in Anteilen am Gesamtsignal')
plt.title('Frequenzuebergang 300Hz -> 700Hz, sampling rate 42kHz')
cbar = plt.colorbar()
cbar.set_label('PSD in dB/Hz')
plt.savefig('hann_3d.png')


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], 20, 1], aspect=('auto'), cmap='hot')
plt.xlabel('Frequenz in Hz')
plt.ylabel('Fensterbreite in Anteilen am Gesamtsignal')
cbar = plt.colorbar()
cbar.set_label('Amplitude in % des hoechsten Werts je Frequenz')
plt.savefig('hann_diff.png')


    

    




(190, 14965)
(14965, 1)
Samplingfrequenz: 42 kHz






    













    

In [2]:

    

from mpl_toolkits.mplot3d import Axes3D

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


    

In [7]:

    




    

In [3]: