13. Aplikace III

• Joseph Fourier
• Fourierova transformace

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()

• 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()