1. ↓ output × 
  2. ############################################################################
#
#   Name:           stochastic_process_07_messung.py
#   Python:         3.3.2 on win32
#   Autor:          Adrian Haas
#   
############################################################################
#  
############################################################################
#
# The script is part of the files:
#
#       stochastic_process_01_titel.py
#       stochastic_process_02_theory.py
#       stochastic_process_03_simulation.py
#       stochastic_process_04_simulation.py
#       stochastic_process_05_messung.py
#       stochastic_process_06_messung.py
#       stochastic_process_07_messung.py
#
# The output video of this files shows a theory, simulation and measurement
# part for a stochastic process. Measurements were made with rtl_sdr. 
# 
############################################################################

import numpy
from subprocess import Popen, PIPE
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from scipy import optimize

from videoplot import*
import cumvar

############################################################################
#Define functions
############################################################################

#PREPARATION

#load frames from rtl_sdr record and take the first frame
channel_293MHz=numpy.load("channel_293MHz.npy")
nr=len(channel_293MHz) #nr of frames
l=len(channel_293MHz[0]) #len of frames
rp=channel_293MHz[0]

#define ranges
t=numpy.arange(0,1,1/l)
historange=numpy.arange(-50,50,0.1)


#CUMULATIVE FUNCTIONS

mean_cum_div=numpy.arange(1,len(rp)+1,1)
mean_cum=numpy.cumsum(rp)/mean_cum_div

var_cum=cumvar.cumvar(mean_cum,rp)


#MODEL FUNCTION 1

mean=mean_cum[l-1]
var=var_cum[l-1]

#gaussian function var
f_gauss=lambda mu,var,tr: (1/(2*var*numpy.pi)**0.5*
                           numpy.exp(-(tr-mu)**2/(2*var)))

#model 1 output
model1=f_gauss(mean,var,historange)


#MODEL FUNCTION 2

#optimize for the standard deviation instead variance,
#std^2=var, otherwise the optimizer fails

#gaussian function std
f_gauss_std=lambda mu,std,tr: (1/(2*std**2*numpy.pi)**0.5
                               *numpy.exp(-(tr-mu)**2/(2*std**2)))

#two gaussian function std
def f_gauss_sum(t,p):
  mu1, std1, mu2, std2, alpha = p
  return alpha*f_gauss_std(mu1,std1,t)+(1-alpha)*f_gauss_std(mu2,std2,t)

#histogram function
histo=numpy.histogram(rp, bins=historange, density=True)
x=historange[1:] #!output of histogram is 1 sample less

#error function: (two gaussian) - histogram
errfunc = lambda p, x, y: f_gauss_sum(x, p) - y

# Initial guess for the parameters
p0 = [0, 2, 1, 0.4,0.9]

#optimizer
p1, success = optimize.leastsq(errfunc, p0[:], args=(x, histo[0]))

#model 2 output
model2=f_gauss_sum(x,p1)


############################################################################
#Animation function
############################################################################

def ani(i):
  #using videoplot module 
  ax2.set(i,0,fig)

  ax3.fade_in(i,25,50,fig)
  ax4.fade_in(i,25,50,fig)

  ax2_animater.reset(i,50)
  ax2_animater.animate(i,50,125,ax2)
  ax3_animater.animate(i,50,125,ax3)
  ax4_animater.animate(i,50,125,ax4)

  ax5.fade_in(i,125,150,fig)
  ax6.fade_in(i,125,150,fig)

  ax2_animater.reset(i,150)
  ax3_animater.reset(i,150)
  ax4_animater.reset(i,150)

  ax2_animater.animate(i,150,225,ax2)
  ax3_animater.animate(i,150,225,ax3)
  ax4_animater.animate(i,150,225,ax4)
  ax5_animater.animate(i,150,225,ax5)

  ax7.set(i,250,fig)
  ax5_animater.set_model1(i,250,model1[1:],ax5)

  ax8.set(i,325,fig)
  ax9.set(i,325,fig)
  ax5_animater.set_model2(i,325,model2,ax5)

  ax2.fade_out(i,400,425,fig)       
  ax3.fade_out(i,400,425,fig) 
  ax4.fade_out(i,400,425,fig) 
  ax5.fade_out(i,400,425,fig)
  ax6.fade_out(i,400,425,fig) 
  ax7.fade_out(i,400,425,fig) 
  ax8.fade_out(i,400,425,fig) 
  ax9.fade_out(i,400,425,fig)

  return


############################################################################
#Init plot 
############################################################################

#define figure without frame and define size
fig = plt.figure(frameon=False)

fig.set_size_inches(19.2,10.8) #1 inch 100 pixel

# background plot
ax=plt.Axes(fig, [0., 0., 1., 1.])
ax.set_axis_bgcolor('black')

#define all artists

rp2=rp/numpy.max(rp)
ax2_animater=c_loadsignal(rp2,t)
ax2=signalbox(fig, [0.1, 0.55, 0.3, 0.4],rp2,t)
ax2.set_title(r"Normed frame 293 MHz",color="white",fontsize=20)
ax2.set_ylim(-1.1,1.1)
ax2.set_axis_bgcolor("black")
ax2.spines['top'].set_color("white")
ax2.spines['right'].set_color("white")
ax2.spines['bottom'].set_color("white")
ax2.spines['left'].set_color("white")
ax2.sig[0].set_color("yellow")

ax3_animater=c_loadsignal(mean_cum,t)
ax3=signalbox(fig, [0.1, 0.3, 0.3, 0.2],ax3_animater.sigout,t)
ax3.set_title(r"Cumulative samples of mean $\mu$",color="white",fontsize=20)
ax3.set_ylim(numpy.min(mean_cum),numpy.max(mean_cum))
ax3.set_axis_bgcolor("black")
ax3.spines['top'].set_color("white")
ax3.spines['right'].set_color("white")
ax3.spines['bottom'].set_color("white")
ax3.spines['left'].set_color("white")
ax3.sig[0].set_color("yellow")

ax4_animater=c_loadsignal(var_cum,t)
ax4=signalbox(fig, [0.1, 0.05, 0.3, 0.2],ax4_animater.sigout,t)
ax4.set_title(r"Cumulative samples of variance $\sigma^2$",
              color="white",fontsize=20)
ax4.set_ylim(numpy.min(var_cum),numpy.max(var_cum)*1.1)
ax4.set_axis_bgcolor("black")
ax4.spines['top'].set_color("white")
ax4.spines['right'].set_color("white")
ax4.spines['bottom'].set_color("white")
ax4.spines['left'].set_color("white")
ax4.sig[0].set_color("yellow")

ax5_animater=c_loadhisto_2(rp,t,historange)
ax5_animater.set()
multi=numpy.vstack((ax5_animater.histo[0],ax5_animater.model1[1:]-1,
                    ax5_animater.model2[1:]-1))
ax5=signalbox(fig, [0.5, 0.4, 0.4, 0.5],multi.T,historange[1:])
ax5.set_title("Probability density function",color="white",fontsize=20)
ax5.set_ylim(0,numpy.max(histo[0])*1.05)
ax5.set_axis_bgcolor("black")
ax5.spines['top'].set_color("white")
ax5.spines['right'].set_color("white")
ax5.spines['bottom'].set_color("white")
ax5.spines['left'].set_color("white")
ax5.sig[0].set_color("yellow")
ax5.sig[1].set_color("cyan")
ax5.sig[1].set_linewidth(3)
ax5.sig[2].set_color("red")
ax5.sig[2].set_linewidth(3)

ax6=textbox(fig, [0.5, 0.3, 0.02, 0.02],"Histogram",20,"yellow")

ax7=textbox(fig, [0.5, 0.25, 0.02, 0.02],
            r"Model 1 with $\mu$:"+str(round(mean,3))+", $\sigma^2$:"
            +str(round(var,3)),20,"cyan")

ax8=textbox(fig, [0.5, 0.2, 0.02, 0.02],
            r"Model 2 with $\alpha$:"+str(round(p1[4],3))+"$\mu_{1}$:"
            +str(round(p1[0],3))+", $\sigma_{1}^2$:"+str(round(p1[1]**2,3))
            +"$\mu_{2}$:"+str(round(p1[2],3))+", $\sigma_{2}^2$:"
            +str(round(p1[3]**2,3)),20,"red")

ax9=textbox(fig, [0.5, 0.15, 0.02, 0.02]
            ,"          determined with least square optimization",20,"red")


fig.add_axes(ax) #add after definition of all artists!


############################################################################
#animate and save video
############################################################################


movie_ani = animation.FuncAnimation(fig, ani, frames=426,interval=1,
                                    blit=False, repeat=False)

Writer = animation.writers['ffmpeg']
writer = Writer(fps=25, metadata=dict(artist='Haas Adrian'),
                bitrate=45675,extra_args=['-vcodec', 'libx264',
                                          '-pix_fmt', 'yuv420p'])

movie_ani.save('video_stochastic_process_07.mp4', writer=writer)

  3. ↑ script ×