#!/usr/bin/env python # -*- coding: utf-8 -*- import matplotlib matplotlib.rcParams['backend'] = 'Agg' matplotlib.rc('font', size='10') matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!! from matplotlib.pyplot import * import os #import scipy.signal from scipy import signal from scipy.signal import fftconvolve from multiprocessing import Process, Pool, cpu_count import time from numpy import * from CMWT import * from SoundGenerator import * from pygolem_lite import Shot from pygolem_lite.modules import list2array,deconvolveExp,save_adv,saveconst from scipy.stats.mstats import mquantiles #from scipy.stats.mstats import mquantiles #RingCoilOrientation = (1,1,-1,-1, 1,-1,1,1,1,-1,-1,1,1,1,1,-1) #TODO důležité #RingCoilOrientation = (-1,1,1,1,-1,1,-1,-1,-1,1,-1,1,-1,-1,-1,-1) #TODO důležité RingCoilOrientation = (-1,-1,1,1,-1,1,-1,1,-1,1,-1,-1,1,-1,-1,-1) RingCoilOrientation = (-1,1,1,1,-1,1,-1,-1,-1,1,-1,1,-1,-1,-1,-1) #AEff = (68e-4, 140e-4, 138e-4, 140e-4, 68e-4, 134e-4,134e-4, 142e-4, 67e-4, 142e-4, 140e-4, 138e-4, 76e-4, 142e-4, 139e-4, 139e-4) #in m^2 AEff = [ 1.56253461e-02, 1.50618720e-02 , 1.45473064e-02 , 1.27421132e-02, 1.44558132e-02, 1.80693924e-02 , 1.56932626e-02 , 1.05514271e-02, 3.90988379e-03 , 5.96870657e-03 , 8.48762208e-03 , 8.47233307e-05, 1.70823107e-02 , 8.08770169e-03 , 1.44213854e-02 ,1.58269600e-02] AEff = [ 0.00209134 , 0.00446023 , 0.00428289 ,0.00374073 , 0.00320856 , 0.00182771, 0.00267794, 0.00284156 , 0.00145525 ,0.00222005 ,0.00206269 , 0.00108452, 0.00054154 , 0.00078395 , 0.00186135 , 0.00274475] #AEff = (68e-4, 140e-4, 138e-4, 140e-4, 130e-4, 134e-4,134e-4, 142e-4, 67e-4, 142e-4, 140e-4, 140e-4, 140e-4, 142e-4, 139e-4, 139e-4) #in m^2 MirnovCoil = (1,1,1,1)#TODO důležité 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 mad(x,axis = -1): #Median absolute deviation m = np.median(x,axis=axis) x = np.swapaxes(x, 0, axis) s = np.median(abs(x-m),axis=0)*1.4826 return s def stft(x, fs, framesz, hop): framesamp = int(2**np.ceil(np.log2(framesz*fs))) hopsamp = int(hop*fs) print hopsamp, framesamp, len(x) print xrange(0, len(x)-framesamp, hopsamp) w = np.hamming(framesamp) X = np.array([np.fft.rfft(w*x[i:i+framesamp]) for i in xrange(0, len(x)-framesamp, hopsamp)], copy=False) return X #def PrepareData(data,typ): #sign = MirnovCoil #WARNING závisí na aktuálním nastavení cívek!!TODO dát to do externího configu #n = size(data, 0) #envelope = zeros(shape(data)) #dt = data[0,1]-data[0,0] #n_smooth = int(0.002/dt)*2+1 #omezení časového rozlišení pro nízké frekvence 250Hz #for i in range(1,n): #data[i,:]*= sign[i-1] #baseline = fftconvolve((data[i,:]), ones(n_smooth)/n_smooth,mode = 'same' ) #data[i,:]-= baseline #envelope[i,:] = fftconvolve(abs(data[i,:]), ones(n_smooth/10)/n_smooth/10,mode = 'same' ) #plot(baseline, label = str(i)) #legend() ##show() #savefig('data.png') #close() #total_envelope = mean(envelope,0) #for i in range(1,n): #data[i,:]*= total_envelope/(envelope[i,:]+1e-6) #return data def plotData( data): dt = data[0,1]-data[0,0] n_smooth = 20 #omezení časového rozlišení pro nízké frekvence 250Hz n = size(data, 1) for i in range(1,size(data, 0)): #baseline = fftconvolve((data[i,:]), ones(n_smooth)/n_smooth,mode = 'same' ) dfft = fft.rfft(data[i,:]) dfft[n/16:] = 0 difft = fft.irfft(dfft) #plot(data[0,:], data[i,:]) plot(data[0,:], difft) xlim(0.016, 0.0171) ylim(-0.2, 0.2) #show() savefig('./graphs/raw_data.png') close() def loadconst(fname): with open(fname, 'r') as fhandle: return float(fhandle.readline()) #return the raw string #def LoadData(): #Data = Shot() #gd = Shot().get_data #das, m1 = gd('any', 'mirnov_1', return_channel = True) #das, m5 = gd('any', 'mirnov_5', return_channel = True) #das, m9 = gd('any', 'mirnov_9', return_channel = True) #das, m13 = gd('any', 'mirnov_13', return_channel = True) #Papouch = list2array( Data[das, [m1, m5, m9, m13]] ).T #plasma = Data['plasma'] #plasma_start = Data['plasma_start'] #plasma_end = Data['plasma_end'] #return Papouch,plasma_start, plasma_end,plasma def LoadData(): Data = Shot() #(data, axis=-1, type='linear', bp=0) gd = Shot().get_data das, data = gd('any', 'ring_1', return_channel = True) Papouch = list2array( Data[das, range(16)] ).T Papouch[1:,:]*= array(RingCoilOrientation)[:,None] Papouch[1:,:]/= array(AEff)[:,None] #print 'ploz' plasma = Data['plasma'] plasma_start = Data['plasma_start'] plasma_end = Data['plasma_end'] #plot() #plot markovicuv profil #X = signal.detrend(Papouch[1:,:], axis=1) #X = cumsum(X,axis=1)*1e-6 #t = Papouch[0,:] #for i in range(16): #plot(t,X[i]) #savefig('data%d.png'%i) #clf() #exit() #X = X[:,16000:16700] #t = t[16000:16700] #X = signal.detrend(X, axis=1, type='linear', bp=0) #print mad(X,axis=1)/mean(mad(X,axis=1))*array(AEff) #exit() #det = range(1,16) #from scipy.stats.mstats import mquantiles #lim = mquantiles(abs(X),0.99) #imshow(X, aspect='auto', extent=[t[0]*1e3,t[-1]*1e3,det[-1]+0.5, det[0]+0.5], #vmin=-lim, vmax=lim) #axis([t[0]*1e3,t[-1]*1e3 ,det[0]+0.5, det[-1]+0.5]) #colorbar(format='%.1e') #savefig('X.png') return Papouch,plasma_start, plasma_end,plasma #def LoadData(path, shot_num, skip_det_num): #path += str(shot_num)+'/' #try: #data = load(path+'data.npy') #if size(data,1) == 0: #raise "wrong data" #except: ##print 'reload' #data = list() #tvec = None #for i in range(0,16): #print 'load NIturbo_%2.2i' %(i+1) #if i in skip_det_num: #data.append(zeros(shape(data[0]))) #continue ##single_data = loadtxt(path+'NIturbo_%2.2i' %(i+1)+'.asc', usecols = (1,)) #single_data = loadtxt(path+'NIturbo_%2.2i' %(i+1), usecols = (1,)) #if tvec == None: ##tvec = loadtxt(path+'NIturbo_%2.2i' %(i+1)+'.asc', usecols = (0,)) #tvec = loadtxt(path+'NIturbo_%2.2i' %(i+1), usecols = (0,)) #data.append(single_data) #data = array(data) #if size(data,1) == 0: #raise "wrong data" #data = vstack((tvec, data)) #save(path+'data', single(data)) #try: #start = loadconst(path+'PlasmaStart')/1000 #end = loadconst(path+'PlasmaEnd' )/1000 #except: #start = nan #end = nan #print start,end #return data,start,end #def LoadData_npz(path, shot_num, skip_det_num): #FIXME provizorní verze!! #data = load('Nidatap.npz') ##data data['data'] #tvec = linspace(data['t_start'], data['t_end'], size(data['data'],0)) #data = data['data'] #data[:,skip_det_num] = 0 #data = vstack((tvec, data.T)) #start = data[0,8300] #end = data[0,19260] #return data,start,end ##def LoadData_mirnov(path, shot_num): ##path += str(shot_num)+'/' ##try: ##data = load(path+'data.npy') ##if size(data,1) == 0: ##raise "wrong data" ##if shot_num > 6000 and shot_num < 9500: ##data = data.T ##data[1:,1:] = diff(data[1:,:], axis = 1)/(data[0,1]-data[0,0]) ##except: ##data = list() ##tvec = None ##for i in range(0,4): ##print 'load PapouchSt_%2.2i' %(i+1) ##single_data = loadtxt(path+'PapouchSt_%2.2i' %(i+1), usecols = (1,)) ##if tvec == None: ##tvec = loadtxt(path+'PapouchSt_%2.2i' %(i+1), usecols = (0,)) ##data.append(single_data) ##data = array(data) ##if size(data,1) == 0: ##raise "wrong data" ##data = vstack((tvec, data)) ##save(path+'data', single(data)) ##try: ##start = loadconst(path+'PlasmaStart')/1000 ##end = loadconst(path+'PlasmaEnd' )/1000 ##except: ##start = nan ##end = nan ##print start,end ##return data,start,end def PlotSpec((freq, field, t, sufname ,vmin, vmax)): fig = figure() t = t*1000 ax = fig.add_axes([0.1, 0.1, 0.8, 0.85]) ax.set_yscale('log', nonposy='clip') img = ax.imshow(abs(field), extent=[t[0],t[-1] ,freq[-1], freq[0]], aspect='auto',vmin=vmin, vmax=vmax) minorLocator = MultipleLocator(1) ax.xaxis.set_minor_locator(minorLocator) ax.axis([t[0],t[-1] ,amin(freq), amax(freq)]) ax.set_xlabel('time [ms]') ax.set_ylabel('Frequency [Hz]') savefig('graphs/spectrogram'+str(sufname)+'.png',bbox_inches='tight') close() def PrepareData( t, X): #print X.shape #fig = figure('poloidal M') #t = t #ax = fig.add_axes([0.1, 0.1, 0.8, 0.85]) X = cumsum(X.T,axis=1)*1e-6 N = 101 taps =signal.firwin(N, 1e3/(0.5*1e6) , window='hamming') filtered_x = signal.lfilter(taps, 1.0, X[:,:]-X[:,0][:,None], axis=1) taps =signal.firwin(N, 1e5/(0.5*1e6) , window='hamming') X = X[:,:-N/2]-X[:,0][:,None]- filtered_x[:,(N+1)/2:] t = t[:-N/2] X = signal.lfilter(taps, 1.0, X, axis=1) X = X[:,(N+1)/2:] t = t[:-N/2] #print X.shape #exit() return t, X.T def PlotSignal(t, sig, tbeg, tend): ind = (t> tbeg)&(t < tend) t = t[ind]*1000 sig = sig[ind,:] sig = signal.detrend(sig, axis=0) sig/= std(sig, axis=0) lim = mquantiles(abs(sig),0.99) det = range(1,16) imshow(sig.T, aspect='auto', extent=[t[0],t[-1],det[-1]+0.5, det[0]+0.5], vmin=-lim, vmax=lim, interpolation='nearest') axis([t[0],t[-1] ,det[0]+0.5, det[-1]+0.5]) colorbar(format='%.1e') savefig('graphs/signal.png') clf() #exit() def PlotModes(field, t): fig = figure('poloidal M') t = t*1000 ax = fig.add_axes([0.1, 0.1, 0.8, 0.85]) print field.shape #exit() X = cumsum(field.T,axis=1)*1e-6 #t = Papouch[0,:] #X = signal.detrend(X, axis=1, type='linear', bp=0) print mad(X,axis=1)/mean(mad(X,axis=1)) #X = X[:,16000:16700] #t = t[16000:16700] #FC = 2e3/(0.5*1e6) #N = 3001 # number of filter taps #b = signal.firwin(N, cutoff=FC, window='hamming') #b/= mean(b) #plot(b) #savefig('b.png') #X -= signal.lfilter(b, 1, X, axis=0) #print mean(abs(X)), mean(abs(signal.lfilter(b, 1, X, axis=0) )) #FC = 5e4/(0.5*1e6) #b = signal.firwin(N, cutoff=FC, window='hamming') #clf() #plot(b) #savefig('graphs/f.png') #show() #X = signal.lfilter(b, 1, X, axis=1) #exit() N = 101 taps =signal.firwin(N, 1e3/(0.5*1e6) , window='hamming') #plot(taps) #savefig('X.png') #exit() filtered_x = signal.lfilter(taps, 1.0, X[:,:]-X[:,0][:,None], axis=1) taps =signal.firwin(N, 1e5/(0.5*1e6) , window='hamming') X = X[:,:-N/2]-X[:,0][:,None]- filtered_x[:,(N+1)/2:] t = t[:-N/2] X = signal.lfilter(taps, 1.0, X, axis=1) X = X[:,(N+1)/2:] t = t[:-N/2] #exit() #plot(taps) #savefig('X.png') #exit() #taps =signal.firwin(N, 1e3/(0.5*1e6) , window='hamming') #filtered_x = signal.lfilter(taps, 1.0, X[1,:]-X[1,0]) ##print X[1,:] #clf() #plot((X[:,:]-X[:,0][:,None]).T) #plot(filtered_x[:,N/2:].T) #plot(X.T) #savefig('X.png') #exit() #X = X[:,5300:5700] #t = t[5300:5700] #X = signal.detrend(X, axis=1, type='linear', bp=0) #print mad(X,axis=1)/mean(mad(X,axis=1)) #lim = mquantiles(abs(abs(X)),0.99) #det = range(1,16) #imshow(X, aspect='auto', extent=[t[0],t[-1],det[-1]+0.5, det[0]+0.5], #vmin=-lim, vmax=lim, interpolation='nearest') #axis([t[0],t[-1] ,det[0]+0.5, det[-1]+0.5]) #colorbar(format='%.1e') #savefig('X.png') #clf() #exit() ##fs = 1e6 ##framesz = 0.001 # with a frame size of 50 milliseconds ##hop = 0.0002 # and hop size of 20 milliseconds. ##print shape(X[10,:]) ##FX = stft(X[10,:]+sin(arange(len(X[10,:]))/1000)*10000, fs, framesz, hop) ##print abs(FX) ##print shape(FX) ##imshow(log(abs(FX)).T, aspect='auto') ##savefig('STFT.png') #exit() #imshow(X, aspect='auto') #savefig('graphs/X.png') #field = X #print shape(field) #exit() Fsig = abs(fft.rfft(X, axis = 0)) Fsig = signal.medfilt(Fsig, [ 1,101]) #imshow(abs(Fsig), aspect='auto', vmin = 0, vmax = lim) #savefig('X2.png') #exit() #FC = 1e3/(0.5*1e6) #N = 1001 # number of filter taps #b = signal.firwin(1001, cutoff=FC, window='hamming') #Fsig = signal.lfilter(b, 1, Fsig, axis=0) n_mod = size(Fsig,0) #modes = arange(n_modn_mod) Fsig = c_[Fsig[:,(n_mod+1)/2:],Fsig[:,:(n_mod+1)/2] ] #print Fsig.shape lim = mquantiles(abs(Fsig),0.99) img = ax.imshow(Fsig, extent=[t[0],t[-1] ,n_mod/2+0.5, -n_mod/2+0.5], aspect='auto', interpolation='nearest',vmin=0, vmax=lim) fig.colorbar(img) minorLocator = MultipleLocator(1) ax.xaxis.set_minor_locator(minorLocator) ax.yaxis.set_major_locator(minorLocator) ax.axis([t[0],t[-1] ,-n_mod/2+0.5, n_mod/2+0.5]) ax.set_xlabel('time [ms]') max_mode = argmax(mean(Fsig,axis=0))-n_mod/2 #print max_mode-n_mod/2 saveconst('mode_M', max_mode) saveconst('mode_M_abs', max(mean(Fsig,axis=0)) ) ax.set_ylabel('Poloidal mode number M') #savefig('graphs/poloidalM.png') savefig('graphs/poloidalM.png',bbox_inches='tight') close() #print 'saved' exit() def CalculateSpectrogram(): print "CalculateSpectrogram" t1 = time.time() signal,plasma_start,plasma_end,plasma = LoadData() #plasma = False #BUG!!!! #plasma_start = plasma_start*0.8 #plasma_end = min((plasma_end-plasma_start)*1.2+plasma_start,signal[0,-1]) #plot(signal.T) #savefig('graphs/data.png') #clf() if not plasma: plasma_start = 0 plasma_end = 0.04 dt = signal[0,1]-signal[0,0] n_start = argmin(abs(plasma_start - signal[0,:])) n_end = argmin(abs(plasma_end - signal[0,:])) t = signal[0,n_start:n_end] signal = signal[1:,n_start:n_end].T #plot(signal) #savefig('data.png') #clf() PlotModes(signal, t) exit() t, signal = PrepareData( t, signal) PlotSignal(t, signal, 16.3e-3, 16.7e-3) omega0 = 40#20 horiz_res = 2000 f_min = 1e3 #Hz f_max = 100e4 #Hz # !! searched MHD mdoes !! modes = [-4, -3,-2,-1,0,1,2,3, 4] N = size(signal,0) n_det = size(signal,1) t2 = time.time() signal_fft = fft.rfft(signal, axis = 0) for m in modes: print "wave order", m signal_m = copy(signal_fft) phase = arange(n_det)/double(n_det)*m exp_phase = matrix(exp(-2*pi*1j*phase)) signal_m*= exp_phase signal_m = sum(signal_m, axis = 1) signal_m = fft.irfft(signal_m) try: soundGenerator(signal_m,2000, 'mp3/sound'+str(m)) except Exception, e: print "sound gener. failed err:" , e.message print 'generated sound', time.time()-t2 #modes = [2,] spec_all = list() scale_all = list() print "started wavelets" #out = map(NTM_CWT, [ (signal, dt, 0.005,omega0,m ,horiz_res, f_min,f_max) for m in modes ]) p = Pool(cpu_count()) out = p.map(NTM_CWT, [ (signal, dt, 0.005,omega0,m ,horiz_res, f_min,f_max) for m in modes ]) p.close() p.join() for i in range(len(modes)): spec, scale = out[i] spec = single(abs(spec)) spec_all.append(spec) scale_all.append(scale) print 'calc time', time.time()-t1 t_plot = time.time() spec_all = array(spec_all) scale_all = array(scale_all) freq = (omega0 + sqrt(2.0 + omega0**2))/(4*pi * scale_all) contrast = 10 spec_all = log(1+contrast*abs(spec_all)) vmin = amin(spec_all) vmax = amax(spec_all) p = Pool(cpu_count()) p.map( PlotSpec, [(freq[i,...],spec_all[i,...],t, m,vmin, vmax) for i,m in enumerate(modes) ]) p.close() p.join() print 'plot. time', time.time()-t_plot def main(): for path in ['graphs', 'mp3']: if not os.path.exists(path): os.mkdir(path) if sys.argv[1] == "plots": CalculateSpectrogram() os.system('convert -resize 150x120\! graphs/spectrogram0.png icon.png') saveconst('status', 0) if __name__ == "__main__": main()