#!/usr/bin/python3

print('*** AZ/EL vrtiljaka in zvočni spektralni analizator - S53MV 07.05.2022 ***')
cime='/dev/ttyUSB0'	#ime zaporednega vmesnika vrtiljaka KR5600
fs=96000		#Frekvenca vzročenja zvočnega signala (max 192000Hz)
pts=1024		#Število točk spektra 0...fs/2 (potenca 2 za učinkovit FFT)
avg=256			#Število povprečenj spektra
el0=16			#Konstante pretvorbe podatkov vrtiljaka (cela števila)
el90=130
azm180=17
az180=247

import sounddevice as sd        #Poiskati na spletu in naložiti s pip3!
import numpy as np              #Uporaba učinkovitih funkcij numpy
import matplotlib.pyplot as plt #Risanje rezultata z matplotlib
import serial
sd.default.device=0             #Hardverska izbira vhoda
com=serial.Serial(cime,9600,timeout=0)	#odpri comport,9600bps,8bit,noparity

Hzkm=1627   #Faktor pretvorbe frekvence v dolžino [km]
fskala=np.linspace(0,fs/2/Hzkm,pts)     #Izračunaj dolžinsko skalo [km]

def spekter(fs,pts,avg,fskala): #Položaj vrtiljaka, zajem signala z zvočno kartico in povprečenje spektra [dB]
	com.write(bytes([80]))	#Vprašaj za parametre vrtiljaka KR5600
	s=sd.rec(pts*(1+avg),samplerate=fs,channels=1)  #Vzorčenje ADC zvočne kartice, 2D polje!
	sd.wait()		#Počakaj do konca vzorčenja
	azel='Vrtiljak '	#Prečitaj azimut in elevacijo vrtiljaka KR5600
	v=com.read(99)
	if len(v)!=20:		#Napaka pri prenosu podatkov?
		azel+='ni podatka!!!'
	else:
		azb=int(v[3])	#Byte azimut
		elb=int(v[4])	#Byte elevacija
		azel+=str(azb)+'/'+str(elb)
		az=(azb-azm180)*360.0/(az180-azm180)-180	#Preračunaj azimut
		if az<0:
			az+=360
		el=(elb-el0)*90.0/(el90-el0)			#Preračunaj elevacijo
		azel+=str(' AZ=%05.1f'%az)+str('/EL=%04.1f'%el)
	print(azel, end=' ')	#Izpis azimuta in elevacije
	nz=1.e-12		#Dodatek za neničelni argument logaritma
	s2=np.square(s)		#Izpis jakosti/max [dB]
	print('Jakost:'+"%.1f"%(10*np.log10(nz+np.mean(s2)))+'dB Vrh:'+"%.1f"%(10*np.log10(nz+np.amax(s2))),'dB',end=' ')
	o=1+np.cos(np.linspace(-np.pi,np.pi,2*pts))	#Okno dvignjeni kosinus
	k=0
	s2=np.zeros(pts)
	while k<pts*avg:				#Seštevanje moči FFT
		s2=s2+np.square(np.abs(np.fft.rfft(o*s.reshape(-1)[k:k+2*pts])[:-1]))
		k+=pts
	s=10*np.log10(nz+0.5*s2/pts/pts/avg)		#Povprečje [dB]
	print('Max:'+"%.1f"%(np.amax(s))+'dB @'+"%.1f"%(fskala[np.argmax(s)])+'km')
	return s

wfbp=None
while wfbp!=True:		#Prvo risanje spektra z normalizacijo, ponovi s klikom miške
	fig,ax=plt.subplots()
	line,=ax.plot(fskala,spekter(fs,pts,avg,fskala))
	wfbp=plt.waitforbuttonpress(0.01)
	while wfbp==None:	#Zanka risanja spektra, izhod=pritisk tipke ali miške
		line.set_ydata(spekter(fs,pts,avg,fskala))
		fig.canvas.draw()
		wfbp=plt.waitforbuttonpress(0.01)
	plt.close()
com.close()		#Zapri com
print('*** Konec ***')	#pritisk tipke na tipkovnici