====== 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
==== 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í.
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')
==== Interakce SOAE s rozmítaným sinem (swept-sine) ====
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()
==== 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
# 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()
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()