Source code :: main
[Return]
[Download]#!/usr/bin/env python
# -*- coding: utf-8 -*-
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
from scipy.signal import *
from numpy import *
from matplotlib.pyplot import *
import numpy as np
from cwt import cwt
from matplotlib.ticker import MultipleLocator, FormatStrFormatter,LogLocator
import os
matplotlib.rcParams['xtick.direction'] = 'out'
matplotlib.rcParams['ytick.direction'] = 'out'
matplotlib.rcParams['xtick.major.size'] = 10
matplotlib.rcParams['xtick.minor.size'] = 7
matplotlib.rcParams['ytick.major.size'] = 10
matplotlib.rcParams['ytick.minor.size'] = 7
def Spectrogram(omega0,frequenciCutOff,signal,name,start, end):
print name
tvec = signal[:,0]
startAdv = max(start*0.8, tvec[0])
endAdv = min((end-start)*1.2+start, tvec[-1])
dt = tvec[1]-tvec[0]
N = len(tvec)
for i in range(3):
filtered = medfilt(signal[:,1], floor(N/100)+1)
res = abs(signal[:,1] - filtered)
ind = argsort(res)
n_out = floor(N/400)
ind_2 = in1d(range(N), ind[-n_out:]) & (res > 0.2*amax(filtered))
signal[ ind_2, 1] = filtered[ind_2] #remove "n_out" most distant points from moving median
spec, scale = cwt(signal[:,1], dt, dj=0.05, wf='morlet', p=omega0, extmethod='none', extlength='powerof2')
# approximate scales through frequencies
freq = (omega0 + np.sqrt(2.0 + omega0**2))/(4*np.pi * scale[1:])
contrast = 10
#spec = spec[:,(tvec > startAdv )*( tvec <endAdv)]
pole = log(1+contrast*np.abs(spec))
pole = pole[freq > frequenciCutOff, : ]
freq = freq[freq > frequenciCutOff]
fig = figure()
ax = fig.add_axes([0.1, 0.1, 0.8, 0.85])
img = ax.imshow(pole, extent=[startAdv*1000,endAdv*1000 ,freq[-1], freq[0]], aspect='auto')
ax.set_yscale('log', nonposy='clip')
minorLocator = MultipleLocator(1)
ax.xaxis.set_minor_locator(minorLocator)
ax.plot([start*1000, (start)*1000], [amin(freq), amax(freq)], 'w--')
ax.plot([end*1000, (end)*1000], [amin(freq), amax(freq)], 'w--')
ax.axis([startAdv*1000,endAdv*1000,amin(freq), amax(freq)])
ax.set_xlabel('time [ms]')
ax.set_ylabel('Frequency [Hz]')
savefig('spectrogram_'+name+'.png')
clf()
signal= signal[(tvec > startAdv )*( tvec<endAdv),1 ]
plot(1000*tvec[(tvec > startAdv )*( tvec<endAdv) ],signal,'k',linewidth=0.1)
#axvline(linewidth=4)
minimum = -max(abs(amin(signal)), abs(amax(signal)))
maximum = -minimum
axis([startAdv*1000,endAdv*1000,minimum, maximum])
plot([start*1000, start*1000], [minimum, maximum], 'r--')
plot([end*1000, end*1000], [minimum, maximum], 'r--')
xlabel('time [ms]')
ylabel('Intensity [a.u.]')
#show()
savefig('signal_'+name+'.png')
clf()
start = loadtxt('PlasmaStart')
end = loadtxt('PlasmaEnd')
omega0 = 15 #čím větší tím lepší frekvenční rozlišení a opačně
frequenciCutOff = 900
signal = loadtxt('NIdata.lvm')
for i in range(1,size(signal,1)):
try:
name = str(i)
Spectrogram(omega0,frequenciCutOff,signal[:,[0,i]],name,start,end)
except:
print 'Failed mirnov'
raise
os.system('convert -resize 150 spectrogram_1.png icon.png')
#try:
#start = 0
#end = 50
#omega0 = 30
#frequenciCutOff = 900
#signal = loadtxt('comedi_spectra.txt')
#signal[:,1] -=mean(signal[:,1])
#name = 'comedi'
#Spectrogram(omega0,frequenciCutOff,signal,name,start,end)
#except:
#print 'Failed comedi'
#raise
[Return]