############################################################################ # # Name: signal_model.py # Python: 3.3.2 on win32 # Autor: Adrian Haas # ############################################################################ # Signal modeling ############################################################################ # # Creating a signal model. Source: [1] Monson. H. Hayes, STATISTICAL DIGITAL # PROCESSING AND MODELING, The Autocorrelation Method p. 178 ff. The # functions were translated from matlab to python code. Input video length # is set to 5s. # ############################################################################# import numpy from subprocess import Popen, PIPE from pylab import* import struct from scipy import signal as sgn import matplotlib.animation as animation ############################################################################# #functions ############################################################################# def periodogram(sig): #see [1] page 394 N=len(sig) try: re=10*log10((abs(numpy.fft.fft(sig))**2)/(N+1)) except: re=numpy.zeros(N) return re def covar(x,p): #autocorrelation matrix x=signal, p Matrix size # this is an alternative to [1], [1] uses convm p 573 akf=numpy.correlate(x,x,"full") matrix=numpy.fromfunction(lambda i,j: akf[-j+i+int((len(akf)-1)/2)], (p,p),dtype=int) matrixconj=numpy.conjugate(matrix) matrix=numpy.tril(matrix,0)+numpy.triu(matrixconj,1) return matrix def acm(x,p): #autocorrelation method #see [1] page 181 N=len(x) if p>=N: print("ERROR,Model order too large") else: X=covar(x,p+1) Xq=X[:p,:p] x1=X[0,1:p+1] a, residual, rank, singularValues = linalg.lstsq(Xq, -x1) a=numpy.hstack((numpy.array([1]),a)) Xpre=numpy.dot(X,a) err=abs(sum((a*X[0,:]))) return a, err def mem(x,p): #This function create the plot graphs a,e=acm(x,p) fourierlen=len(x) padfaktor=fourierlen-len(a) apad0=numpy.zeros(padfaktor) apad=numpy.hstack((a,apad0)) try: Px=20*(log10(e)-log10(abs(numpy.fft.fft(apad)))) Pnorm=abs(Px)-20*log10(len(x)*numpy.pi) #norm not clear? except: Pnorm=numpy.zeros(len(x)) return Pnorm,a,e ############################################################################# #read audio ############################################################################# movie="model_b.wav" command = [ 'ffmpeg', '-i', movie, '-f', 's16le', '-acodec', 'pcm_s16le', '-ar', '44100', '-ac', '2', '-'] pipe = Popen(command, stdout=PIPE, bufsize=10**8) raw_audio = pipe.stdout.read() pipe.stdout.close() pipe.terminate() #transter to numpy array, framelength 2 means 16bit, convert b to int signal=numpy.frombuffer(raw_audio, dtype='int16') #stero evens are channel 1, odds channel 2 channel1=signal[0::2] channel2=signal[1::2] ############################################################################# #processing ############################################################################# #TP 3kHz filter3000=numpy.array( [-0.0, -3.8709847156706606e-05, -8.586965379592158e-05, 1.0856232397063891e-05, 0.00040674304289794707, 0.0011677015599131234, 0.0021897836926228226, 0.0031647804483967133, 0.003616466873188358, 0.0030130182349014122, 0.0009375383936684646, -0.002722953306815577, -0.007622751401134849, -0.012889280177131902, -0.01717879798817474, -0.018865987705366256, -0.016343983547811174, -0.008380051790808125, 0.005548717954883988, 0.0250194931870736, 0.04859236966197352, 0.07393670584287788, 0.0981463255740456, 0.118192315914501, 0.1314291413460933, 0.1360544217687075, 0.13142914134609338, 0.11819231591450105, 0.09814632557404564, 0.07393670584287804, 0.04859236966197369, 0.02501949318707366, 0.0055487179548841074, -0.008380051790808021, -0.016343983547811174, -0.018865987705366193, -0.017178797988174684, -0.012889280177131866, -0.007622751401134831, -0.002722953306815559, 0.000937538393668481, 0.003013018234901431, 0.003616466873188368, 0.003164780448396719, 0.0021897836926228326, 0.001167701559913129, 0.00040674304289794994, 1.0856232397065834e-05, -8.586965379592114e-05, -3.8709847156707114e-05, -0.0]) #filter a1=sgn.lfilter(filter3000,1,channel1) #downsample 42.2kHz to ~4kz a1=a1[::10] #Periodogram a1p=periodogram(a1) #ACM-Estimation mema1,a1n,a1z=mem(a1,8) ############################################################################# #Animation ############################################################################# #Plot init parameters a1p0=numpy.max(a1p) mema1_min=0.9*numpy.min(mema1) mema1_max=1.1*numpy.max(mema1) fsize=20 time=len(a1) line=numpy.zeros(len(a1)) mema1,a1n,a1z=mem(line,8) #drawing support/animation functions def piticks(length): plt.xticks([0, length/2,length-1], ['$0$', r'$\frac{\pi}{2}$', r'$\pi$'],fontsize=fsize) def halfspec(sigplot): leng=len(sigplot)/2 plt.plot(sigplot[:leng]) plt.xlim(0,leng) piticks(leng) def ani(i): j=i/25/5*time line[:j]=a1[:j] a1p=periodogram(line) mema1,a1n,a1z=mem(line,8) s_w.set_ydata(line) s_w2.set_ydata(a1p) s_w3.set_ydata(mema1) s_w4.set_ydata(a1n) return s_w,s_w2,s_w3,s_w4, #setup plot fig=plt.figure("Animation") fig.set_size_inches(19.2,10.8)#pixel size ax=fig.add_subplot(311) ax.set_title("Zeitfenster B [s]",fontsize=fsize) s_w,=ax.plot(line, linewidth=1) border=1.1*numpy.max(abs(a1)) plt.ylim(-border,border) plt.xlim(0,len(a1)) length=len(a1) plt.xticks([0, length/4, 3*length/4,length-1], ['1','2','4','5'],fontsize=fsize) plt.tick_params( which='both', top='off', left='off', labelleft='off', right='off') ax=fig.add_subplot(312) ax.set_title("Periodogramm",fontsize=fsize) a1p=periodogram(line) s_w2,=ax.plot(a1p,linewidth=1) plt.ylim(0,1.1*a1p0) plt.xlim(0,len(a1)/2) piticks(len(a1)/2) plt.tick_params( which='both', top='off', left='off', labelleft='off', right='off') ax=fig.add_subplot(325) ax.set_title("Autoregressives Modell",fontsize=fsize) s_w4,=ax.plot(a1n,"ro",markersize=12,linewidth=1) plt.ylim(-1.1,1.1) plt.tick_params(axis='both', which='major', labelsize=fsize) plt.tick_params( which='both', top='off', right='off') ax=fig.add_subplot(326) ax.set_title("Schätzung: Autokorrelationsmethode",fontsize=fsize) s_w3,=ax.plot(mema1,linewidth=1) plt.ylim(mema1_min,mema1_max) plt.xlim(0,len(a1)/2) piticks(len(a1)/2) plt.tick_params( which='both', top='off', left='off', labelleft='off', right='off') movie_ani = animation.FuncAnimation(fig, ani, frames=25*5,interval=1, blit=False) ############################################################################# #write video ############################################################################# Writer = animation.writers['ffmpeg'] writer = Writer(fps=25, metadata=dict(artist='Haas Adrian'), bitrate=7200,extra_args=['-vcodec', 'libx264', '-pix_fmt', 'yuv420p']) movie_ani.save('video.mp4', writer=writer) ############################################################################# # add audio to video ############################################################################# #multiplex channels signalout=numpy.zeros(len(channel1)*2,dtype='int16') signalout[0::2]=channel1 signalout[1::2]=channel2 #change int16 to bytes, use the struct module <:little endian, h: integer(2) order='<'+str(len(signalout))+'h' audio_bytes=struct.pack(order, *signalout ) #adding sound to video command = [ 'ffmpeg', '-f', 's16le', '-acodec', 'pcm_s16le', '-ar', '44100', '-ac', '2', '-i', '-', '-y', '-i', 'video.mp4', '-vcodec', 'copy', 'video_sound.mp4'] pipe = Popen(command, stdout=PIPE,stdin=PIPE,stderr=PIPE) pipe.communicate(audio_bytes) #stdin.write instead communicate not work ?! pipe.stdin.close() pipe.terminate() #show()