Source code :: main

#!/usr/bin/python2
# -*- coding: utf-8 -*-
#fASTER VERSION THAN MAIN_2, MORE COMPLICATED


####     Microwaves 2.0
#This algorithm calculate transformation of the sin signal from microwaves density measurement to the phase/amplitude space. 
# First step of the calculation is estimate of the base frequency and calculation of the complex exponential
#with the same frequency.In the second step is signal multiplied by this exponential
#and resulting low frequency signal is smoothed over Gaussian window. Finally complex phase and amplitude are calculated.  

# Authors: Tomas Odstrcil, Ondrej Grover

from time import time
t = time()
#import matplotlib 
#matplotlib.rcParams['backend'] = 'Agg'
#matplotlib.rc('font',  size='10')
#matplotlib.rc('text', usetex=True)  # FIXME !! nicer but slower !!!
#import matplotlib.pyplot as plt

from numpy import *
from scipy.fftpack import fft, ifft,fftfreq
from scipy.signal import fftconvolve, medfilt 
from pygolem_lite import save_adv,load_adv,saveconst, Shot
from pygolem_lite.modules import multiplot, get_data, paralel_multiplot
import os
import sys
from matplotlib.pylab import *

print 'include time ',time()-t 


def BlockConv(signalin, impulse,mode='full'):
    S = len(impulse) # impulse length
    L = len(signalin) # signal length
    n_ext = ((L/S)+1)*S-L
    
    signalin = hstack((signalin, zeros(n_ext)))
    L = len(signalin) # signal length

    # find fftsize as the next power of 2
    # beyond S*2-1

    N = 2**int(log2(S)+2)


    # signal blocks for FFT
    imp = zeros(N,dtype= impulse.dtype)
    sig = zeros(N,dtype=signalin.dtype)

    # output signal
    sigout = zeros((L/2)*2+S+1, dtype=signalin.dtype)

    # copy impulse for FFT
    imp[:S] = impulse
    spec_imp = fft(imp)

    for i in xrange( L/S):

	p = S*i
	# get block from input
	sig[:S] = signalin[p:p+S]
	# perform convolution and overlap-add
	sigout[p:p+2*S] += ifft(spec_imp*fft(sig))[:2*S]

    if mode == 'full':
	sigout = sigout[:-n_ext-1]
    if mode == 'same':
	sigout = sigout[(S-1)/2:-(S+1)/2-n_ext-1]
    return sigout




def Demodulation(y,win,dt):
    t = time()
   
    y-= mean(y, axis = 0)
    n = size(y,0)
    N = 2 ** int(ceil(log2(n)))
    fourier = fft(y[:,0], N) #calulcate the fourier transfrom for the sine data


    #substract the varing offset of the signals
    reduct = 5*999
    n_ = (n/reduct)*reduct
    c = reshape(y[:n_,:],(n_/reduct,reduct,2))
    offset = mean(c,axis=1)
    c -= offset[:,newaxis,:]
    c = swapaxes(c,0,1)
    c = c.reshape(-1,2)


    #find the carrier frequency
    max_frequency_index =  argmax(abs(fourier[:N/2]))
    
    f = fftfreq(N,dt)
    s = slice(max_frequency_index-100,max_frequency_index+100)
    amplitude = abs(fourier[s])
    f_carrier = sum(f[s]*amplitude)/sum(amplitude)
 

    
    #find a unharmonics factor
    amplitude = linalg.norm(amplitude)
    norm1 = linalg.norm(fourier[:N/2])
    
    k = sqrt(norm1**2-amplitude**2)/norm1
    
    
    fourier[:] = 0 #cancel out all other frequencies
    fourier[max_frequency_index] = 1

    
    cmpl_exp = ifft(fourier)[:n]

    gauss = exp(-arange(-3*win,3*win)**2/win**2)  
    gauss/= sum(gauss)
  
    signal = list()
    for i in range(size(y,1)):
	signal.append(BlockConv(y[:,i]*cmpl_exp,gauss,mode='same' ))
	#it can be faster by application of the block convolution 
	
    signal = array(signal, copy = False).T    

    amplitude = abs(signal)
    phase = angle(signal)
    phase = unwrap(phase, axis = 0)

	
    print 'calc. time', time()-t
    return amplitude,phase,f_carrier,k,norm1/(n/2)

def LoadData():
    Data = Shot()
    Bt_trigger = Data['Tb']

    gd = Shot().get_data
    tvec, density1  = gd('any', 'density')
    tvec, density2  = gd('any', 'density_2')
    start  = Data['plasma_start']
    end  = Data['plasma_end']

       
    return tvec, start, end, density1,density2,Bt_trigger


    
    
def graphs():
	    
    tvec, phase_pila = load_adv('results/phase_saw')
    tvec, phase_sinus = load_adv('results/phase_sinus')
    tvec, phase = load_adv('results/phase_substracted')
    tvec, phase_corr = load_adv('results/phase_corrected')

    tvec, amplitude = load_adv('results/amplitude_sinus')   
    tvec, n_e = load_adv('results/electron_density')

    #print mean(n_e)
   
    
    data = [[get_data([tvec,-phase_pila+mean(phase_pila)], 'phase 1', 'phase [rad]', xlim=[0,40], fmt="--"), 
	    get_data([tvec,-phase_sinus+mean(phase_sinus)], 'phase 2', 'phase [rad]', xlim=[0,40], fmt="--" ), 
	    get_data([tvec,phase], 'substracted phase', 'phase [rad]', xlim=[0,40], fmt="k" ), 
	    get_data([tvec,phase_corr], 'corrected phase', 'phase [rad]', xlim=[0,40], fmt="k:" )], 
	    get_data([tvec,amplitude], 'amplitude', 'amplitude [a.u.]', xlim=[0,40],ylim=[0,None]  )]
    multiplot(data, ''  , 'graphs/demodulation', (10,6) )


    data = get_data('electron_density', 'Average electron density', '$<n_e>$ [$10^{19}\,m^{-3}$]', data_rescale=1e-19)
    multiplot(data, ''  , 'graphs/electron_density', (9,3) )
    paralel_multiplot(data, '', 'icon', (4,3), 40)


def main():
    
    for path in ['graphs', 'results' ]:
	if not os.path.exists(path):
	    os.mkdir(path)
	    
    if sys.argv[1] ==  "analysis":
	
	win = 30e-6 #[s]

	t = time()
	#if the phase difference is negativ, change order of dens1,dens2
	try:
	    tvec,start, end, density2,density1,Bt_trigger = LoadData()
	except:
	    return 
	
	dt = (tvec[-1]-tvec[0])/len(tvec)
	density1 = density1[tvec>0]
	density2 = density2[tvec>0]
	tvec = tvec[tvec>0]

	if amax(density1)  >  amax(density2):
	    density1,density2 = density2,density1
	
	
	print 'load time ', time()-t
	signals = vstack((density1,density2)).T
	amplitude,phase,f_carrier,k,norm_ampl = Demodulation(signals,win/dt,dt)  
	#signals/= amplitude
	#amplitude,phase,f_carrier = Demodulation(signals,win/dt,dt)  

	
	downsample = int(win/dt/2)    
	phase_pila = phase[::downsample,0]
	phase_sinus = phase[::downsample,1]
	tvec = tvec[::downsample]
	
	phase = phase_pila-phase_sinus
	phase -= median(phase[tvec<Bt_trigger])
	switched = sign(mean(phase[(tvec > start) & (tvec < end)])) == -1
	if switched:
	    phase *= -1   # rotate the density of cabels were switched 
	#amplitude = amplitude[::downsample,0]
	#else:
	amplitude = amplitude[::downsample,0]
	amplitude *= norm_ampl/median(amplitude)
	apl0 = median(amplitude[tvec < start])
	
	
	
	##############  detekce skoku ################
	t0 = time()
	
	dwin = 20;
	N = len(phase)
	phase_diff = zeros(N)
	for i in arange(N):
	    p_tmp = phase[max(i-dwin,0):min(i+dwin,N-1)]
	    phase_diff[i] = amax(p_tmp) - amin(p_tmp)
	    
	ind = medfilt((amplitude < 0.8*apl0) & (phase_diff > 2), 3)  # remove standalone points
	ind = where(ind)[0]
	ind_skip = where(diff(ind)> 1)[0]
	ind_skip = unique(concatenate([[ind[0]], ind[ind_skip], ind[ind_skip+1] , [ind[-1]]]))  # find indexes with skips

	phase_new = phase.copy()
	if mod(len(ind_skip), 2) == 0 and len(ind_skip) < 20: # fix only the simple issues 
	    Nskip = len(ind_skip)
	    for i in arange(0,Nskip,2):
		i0 = ind_skip[i]
		i1 = ind_skip[i+1]
		#print i0, i1
		phase_new[i1:] += phase_new[i0] - phase_new[i1]
		phase_new[i0:i1] = nan

	print "time detekce skoku", time() - t0

		
	#i0 = ind[0]
	#while i0 < amax(ind):
	    
	##for i0 in ind:
	    #i1 = i0
	    #for i2 in ind[ind > i1]:
		#if i2 - i1 > 1:
		    #break
		#else:
		    #i1 = i2
	    #print  i1, i2
	    #i0 = i1
	    
	
	plot(phase_diff)
	plot(amplitude / amax(amplitude) * amax(phase_diff))
	#plot(ind*amax(phase_diff))
	savefig('diff.png')
	close()
	
	plot(phase)
	plot(phase_new)
	#plot(amplitude / amax(amplitude) * amax(phase_diff))
	#plot(ind*amax(phase_diff))
	savefig('phase.png')
	close()
	
	save_adv('results/phase_saw', tvec, phase_pila)
	save_adv('results/phase_sinus', tvec, phase_sinus)
	save_adv('results/phase_substracted', tvec, phase)
	save_adv('results/amplitude_sinus', tvec, amplitude)    
	save_adv('results/phase_corrected', tvec, phase_new)    

	p = exp(1-mean((norm_ampl/amplitude)**2.5))
	
	from scipy.constants import c,m_e,epsilon_0,e
	ind_plasma = (tvec > start) & (tvec < end)
	a = 0.01   #[m]
	f_0 = 75e9 #[Hz]
	lambda_0 = c/f_0
	n_e = 4*pi*m_e*epsilon_0*c**2/(e**2*lambda_0)*phase
	save_adv('results/electron_density_line', tvec, n_e)
	saveconst('results/electron_density_mean', mean(n_e[ind_plasma]))

	n_e /= 2*a
	save_adv('results/electron_density', tvec, n_e)

	phase_skip  = abs(phase[0] - phase[-1])/2*pi 
	negativity = 1-sum(phase[ind_plasma])/sum(abs(phase[ind_plasma]))
	cum_var = mean(abs(diff(phase[ind_plasma]))) / mean(abs(phase[ind_plasma]))
	
	saveconst('results/carrier_freq', abs(f_carrier))
	saveconst('results/harmonics_distortion', k)
	saveconst('results/norm_ampl', norm_ampl)
	saveconst('results/probability', p)
	saveconst('results/reliability', phase_skip + negativity+cum_var*10 )

	#norm_ampl
	print 'carrier_freq ', f_carrier
	print 'harmonics_distortion ',k
	print 'norm_ampl ',norm_ampl
	print 'probability ',p
	print "cum_var", cum_var
	print "negativity", negativity
	print "phase_skip", phase_skip

	
	


    if sys.argv[1] ==  "plots":
	graphs()
	saveconst('status', 0)



if __name__ == "__main__":
    main()