====== 13. Aplikace III ======

===== Signál v časové a frekvenční oblasti =====

  * [[https://cs.wikipedia.org/wiki/Joseph_Fourier|Joseph Fourier]]
  * [[https://cs.wikipedia.org/wiki/Fourierova_transformace|Fourierova transformace]]


==== Úloha 1. - generace harmonického signálu a výpočet spektra ====

Vygenerujte signál $\cos(2\pi ft)$ a spočítejte a zobrazte jeho spektrum


<code python>

</code>


==== Spontánní otoakustické emise (SOAE) ====

Akustické signály evokované ve vnitřním uchu bez specifické vnější akustické stimulace. To znamená, že tyto signály je možné nahrát mikrofonem vloženým do zvukovodu, aniž bychom ucho budili nějakým konkrétním signálem.

Podle jedné recenzované studie jsou SOAE přítomny u zhruba 70% normálně slyšících lidí. Kromě lidí jsou přítomny také u řady dalších obratlovců. Jejich frekvence je stabilní. Stejné hodnoty bychom nahráli u daného subjektu při opakovaném měření po jistém čase. Ale u malých dětí bývá největší výskyt emisí zhruba nad 3 kHz, kdežto u dospělých je to na frekvencích přibližně nad 1 kHz. Z toho lze usoudit, že se jejich frekvence s věkem může měnit.


Úloha 2. - V adresáři SOAE jsou wav soubory se signály nahranými sondou vloženou do zvukovodu. Signály obsahují SOAE. Vykreslete spektrum těchto signálů. Použijte funkci rfft(), která počitá pouze spektrum pro kladné frekvence. Použijte průměrovnání ve frekvenci pro potlačení šumu a zvýraznění emisí.


<code python>
import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import read


def rfft(x):
    """calculates half of signal spectrum (assumed that signal is real)"""
    if x.ndim==1:  # check whether x is 1D
        N = int(np.ceil(len(x)/2)) # half of the spectrum
        xc = np.fft.fft(x)
        X = 2*xc[0:N]/len(x)  # take first half of the spectrum and normalize
    else:
        print('Input must be a 1D array')
        X = None
    return X

fs, signal = read('SOAE/s015_soae_prave.wav')






</code>


==== Interakce SOAE s rozmítaným sinem (swept-sine) ====


<code python>
import scipy.io
import numpy as np
import matplotlib.pyplot as plt
from scipy import signal

Sxx = 0
for i in range(21):
    filename = 's037/LSWrec_s037_221020_140227_F2b_1000_F2e_4000_L1_35_L2_50F2F1_120_Left_{:02d}.mat'.format(i+1)
    
    mat = scipy.io.loadmat(filename)

    fs = int(mat['fs']) # sampling frequency
    rec_signal = mat['recordedSignal']  # recorded signal by the probe microphone (contains emission)
    
    fx, tx, Sx = signal.spectrogram(rec_signal[:,0],fs,('kaiser',40),2499,600)
    Sxx += Sx


plt.pcolormesh(tx, fx, 20*np.log10(Sxx), shading='gouraud')
plt.xlabel('Time [sec]')
plt.ylabel('Frequency [Hz]')
plt.ylim((0,4e3))
plt.clim((-240,-200))


plt.show()


</code>


==== EEG signal ====

  * Příklad na EEG je převzat z Backyard Brains https://backyardbrains.com/experiments/quantifyyoureeg#prettyPhoto
  * EEG signál byl nahrán pro


<code python>
# eeg signal
# nacte signal z wav souboru, provede podvzorkovani (je to zejmena low frequency signal)
# zobrazi casovy prubeh a spektrogram

import numpy as np
import matplotlib.pyplot as plt
from scipy.io.wavfile import read

# load wav file

file = 'EEG_Alpha_SampleData/TimBrain_VisualCortex_BYB_Recording.wav'
fs, data = read(file)

plt.figure('eeg',figsize=(10,10))
plt.plot(data)

# resample 
length_data=np.shape(data)
length_new=length_data[0]*0.05
ld_int=int(length_new)
from scipy import signal
data_new=signal.resample(data,ld_int)

plt.figure('Spectrogram',figsize=(10,10))
d, f, t, im = plt.specgram(data_new, NFFT= 256, Fs=500, noverlap=250)
plt.ylim(0,90)
plt.colorbar(label= "Power/Frequency")
plt.ylabel('Frequency [Hz]')
plt.xlabel('Time [s]')
plt.show()
</code>


<code python>
import csv
import pandas as pd

# ulozi vektor frekvenci do csv souboru
matrixf=np.array(f).T
np.savetxt('Frequencies.csv', matrixf)
df = pd.read_csv("Frequencies.csv", header=None, index_col=None)
df.columns = ["Frequencies"]
df.to_csv("Frequencies.csv", index=False)

# Pro vyber dat (frekvenci) se vytvori vektor s pozicemi (pro 7 Hz az 12 Hz).         

position_vector=[]
length_f=np.shape(f)
l_row_f=length_f[0]
for i in range(0, l_row_f):
    if f[i]>=7 and f[i]<=12:
        position_vector.append(i)

# spocita prumer hodnot spektra pro frekvence mezi 7 a 12 Hz
        
length_d=np.shape(d)
l_col_d=length_d[1]
AlphaRange=[]
for i in range(0,l_col_d):
    AlphaRange.append(np.mean(d[position_vector[0]:max(position_vector)+1,i]))


# data se pak vyhladi pomoci funkce:    
def smoothTriangle(data, degree):
    triangle=np.concatenate((np.arange(degree + 1), np.arange(degree)[::-1])) # up then down
    smoothed=[]

    for i in range(degree, len(data) - degree * 2):
        point=data[i:i + len(triangle)] * triangle
        smoothed.append(np.sum(point)/np.sum(triangle))
    # Handle boundaries
    smoothed=[smoothed[0]]*int(degree + degree/2) + smoothed
    while len(smoothed) < len(data):
        smoothed.append(smoothed[-1])
    return smoothed   


# vykresleni dat
plt.figure('AlphaRange',figsize=(10,10))
y=smoothTriangle(AlphaRange, 100)
plt.plot(t, y)
plt.xlabel('Time [s]',fontsize=12)
plt.xlim(0,max(t))
plt.show()
</code>