Search
Vygenerujte signál $\cos(2\pi ft)$ a spočítejte a zobrazte jeho spektrum
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')
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 # 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()