#!/usr/bin/python3

print('*** Preprost zvočni spektralni analizator - S53MV 28.09.2020 ***')
print('*** popravek odziva za EliteBook 8540w - Andrej 11.2.2012 ***')
print('*** meritev faznega suma (interpolacija in povprecenje v frekvenci-sai) - Andrej 15.2.2012 ***')

fs=192000    #Frekvenca vzročenja zvočnega signala (max 192000Hz)
pts=2048    #Število točk spektra 0...fs/2 (potenca 2 za učinkovit FFT)
avg=128      #Število povprečenj spektra

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 matplotlib.animation as animation

sd.default.device = 0
fskala=np.linspace(0,fs/2,2*pts+1)[:-1]   #Izračunaj frekvenčno skalo [Hz]

### podatki za kalibracijo sistema
dSB = -62.86+15.83 #dB
L = 1500 #metrov
tau = L*1.4676/3e8

### popravek odziva
odziv = np.loadtxt("5mhz-polinom-2048tock.csv", delimiter=',')
a = odziv[:,1]
popravek = np.max(a[4:-1])-a[4:-1]


def spekter(fs,pts,avg,fskala): #Zajem signala z zvočno kartico in povprečenje spektra [dB]
    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
    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(2*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],n=4*pts)[:-1]))
        k+=pts
    s=10*np.log10(nz+0.5*s2/pts/pts/avg)            #Povprečje [dB]
    print('  Max:',"%.1f"%(np.amax(s[4:-1])),'dB @',"%.1f"%(fskala[np.argmax(s[4:-1])]),'Hz')
    return s[4:-1]+popravek

### izris z animacijo
f=fskala[4:-1]
rbw = np.diff(f)[0]
print(rbw)
initsum = spekter(fs,pts,avg,fskala)
fig,ax = plt.subplots()
raw,pn, = ax.plot(f,initsum,f,initsum)
ax.set_xlim(100, 100000)
ax.set_ylim(-130,-30)
ax.set_ylabel("Fazni šum [dBc/Hz]")
ax.set_xlabel("Frekvenčni odmik [Hz]")
plt.yticks(np.arange(-130,-29,10))
plt.xscale('log')
ax.grid(True)

def meritev(nic):
    s = spekter(fs,pts,avg,fskala)
    Pcal = np.amax(s)
    K = Pcal - 3 - dSB
    fsum = s - K - 3 - 6 - 20*np.log10(np.sin(f*np.pi*tau)) -10*np.log10(rbw)
    raw.set_data(f,s)
    pn.set_data(f,fsum)
    ax.figure.canvas.draw()
    return raw,pn

ani = animation.FuncAnimation(fig,meritev,interval=50)
plt.show()

#shrani csv file
s = spekter(fs,pts,avg,fskala)
Pcal = np.amax(s)
K = Pcal - 3 - dSB
fsum = s - K - 3 - 6 - 20*np.log10(np.sin(f*np.pi*tau)) -10*np.log10(rbw)
data = np.column_stack((f,fsum,s))
np.savetxt('rezultat.csv',data,delimiter=',')
print('*** Konec ***')