-
Notifications
You must be signed in to change notification settings - Fork 0
/
utilfunc.py
400 lines (336 loc) · 11.2 KB
/
utilfunc.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
# coding: utf-8
# Funciones útiles
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sig
def plot_sig(x, fs=1, title='Señal temporal'):
"""
Grafica una señal en el dominio temporal.
Parameters
----------
x : 1darray
array de una dimensión que se desea graficar.
fs : int, opcional
frecuencia de muestreo de x. Por defecto es 1.
title : str, opcional
el título del plot. Por defecto es 'Señal temporal'.
Returns
-------
Figure :
Gráfico 2D de x.
"""
t = np.linspace(0, x.size/fs, x.size)
plt.figure(figsize=(16,7))
plt.plot(t, x)
plt.title(title)
plt.xlabel('Tiempo[s]')
plt.ylabel('Amplitud')
plt.grid()
plt.show()
def stft(x,fs,largo,solap):
'''
Calcula la STFT de una señal.
Subdivide la señal y la multiplica por ventanas de tipo Hann,
aplica la FFT a cada ventana y devuelve información tiempo-frecuencia
de la señal original.
Parameters
----------
x: 1darray
señal temporal
fs: int
frecuencia de sampleo
largo: float
duración de cada ventana en segundos
solap: float
superposicion entre ventanas en segundos
Returns
-------
X: list
STFT de la señal.
'''
N = int(largo*fs) # Largo de ventana en muestras
P = int(solap*fs) # Paso en muestras
n = np.arange(N)
hann = 0.5 - 0.5*np.cos(2*np.pi*n/N) # Ventana Hann
X = []
for i in range(0,x.size-N, P):
y = x[i:i+N]*hann
Y = np.fft.rfft(y)
X.append(Y)
return X
def plot_stft(STFT, fs, T, title='STFT'):
'''
Grafica una STFT.
Devuelve un plot donde el eje horizontal es el tiempo y
el eje vertical es la frecuencia. La amplitud se representa
con pseudocolor cmap = 'magma'.
Parameters
----------
STFT: list
Lista con la STFT a graficar
fs: int
Frecuencia de muestreo de la señal
T: float
Duración en segundos de la señal
title: str, opcional
Título del gráfico. Por defecto "'STFT'"
'''
STFT_MAG = np.asarray(np.abs(STFT))
f = np.linspace(0, fs/2, STFT_MAG.shape[1])
t = np.linspace(0, T, STFT_MAG.shape[0])
plt.figure(1,figsize=(20,10))
plt.pcolormesh(t, f, STFT_MAG.T, cmap='magma')
plt.ylabel('Frecuencia (Hz)')
plt.xlabel('Tiempo (s)')
plt.title(title)
plt.show()
def plot_stft2(STFT, f, t, title='STFT'):
'''
Grafica la STFT. Utilizar para STFT calculada con la función de scipy.
Devuelve un plot donde el eje horizontal es el tiempo y
el eje vertical es la frecuencia. La amplitud se representa
con pseudocolor cmap = 'magma'.
Parameters
----------
STFT: 2Darray
Array con la STFT a graficar
f: array
Eje frecuencial
t: array
Eje temporal
title: str
Título del gráfico. Por defecto 'STFT'
'''
plt.figure(1,figsize=(20,10))
plt.pcolormesh(t, f, np.abs(STFT), cmap='magma')
plt.ylabel('Frecuencia (Hz)')
plt.xlabel('Tiempo (s)')
plt.title(title)
plt.show()
def ssf(x, fs, largo=0.01, solap=None, silencio=0.3, b=1):
'''
Reducción de ruido por substracción espectral.
Aplica reducción de ruido por el método de substracción espectral a una señal.
Requiere que la señal no tenga información útil en el inicio. Procesa la señal
con reducción de ruido residual.
Parameters
----------
x: 1darray
Señal a filtrar
fs: int
Frecuencia de muestreo
largo: float, opcional
Longitud en segundos de la ventana para realizar la STFT. Por defecto largo = 0.01 seg
solap: float, opcional
El solapamiento en segundos entre ventana y ventana en la STFT. Por defecto es la mitad del largo. Debe ser menor
al largo.
silencio: float, opcional
Tiempo en segundos al inicio de la señal donde se estimará el ruido. Se requiere que no haya información útil.
b: float, opcional
Parámetro de sobre-sustracción. Usualmente 1<b<2. Por defecto b=1.
Returns
-------
tiempo: array
Vector de tiempo de la señal
salida: array
Señal filtrada
'''
if solap != None:
solap = int(solap*fs) # Pone el paso en muestras
f, t, X = sig.stft(x, fs, nperseg=int(largo*fs), noverlap=solap) # Calcula STFT
X = X.T # Trasponer la matriz
ventanas = []
fases = []
for i in range(len(X)):
ventanas.append(np.abs(X[i]))
fases.append(np.angle(X[i]))
ruido = np.array(ventanas[0])
muestras = int(silencio * fs) # Cantidad de muestras sin señal útil
M = int(muestras/len(ventanas[0]))
for i in range(1,M):
ruido += ventanas[i]
ruido /= M # Estimador de ruido
ruido *= b # Oversubstract
for i in range(len(ventanas)):
for k in range(len(ventanas[i])):
if ventanas[i][k] > ruido[k]:
ventanas[i][k] -= ruido[k]
else: ventanas[i][k] = 0
# Reducción de ruido residual:
maxres = ventanas[0]
for i in range(1,M):
for k in range(len(ventanas[i])):
if ventanas[i][k] > maxres[k]:
maxres[k] = ventanas[i][k] # Selecciona las magnitudes máximas en el silencio inicial luego de la substracción
for i in range(1,len(ventanas)-1):
for k in range(len(ventanas[i])):
if ventanas[i][k] < maxres[k]:
ventanas[i][k] = min(ventanas[i+1][k], ventanas[i][k], ventanas[i-1][k])
# Elige el valor mínimo en ventanas adyacentes para esa frecuencia
out = []
for i in range(len(ventanas)):
out.append(ventanas[i] * np.exp(1j*fases[i])) # Agregamos la magnitud y la fase
out = np.asarray(out)
tiempo, salida = sig.istft(out.T,fs, noverlap=solap) # Antitransformamos
return tiempo, salida
def SNR(x, sigma=None, inicio=0, fin=100):
'''
Calcula la relación señal a ruido (SNR).
Calcula la SNR de una señal con media 0según una porción
de señal donde se estima la desviación estandar. Se define como:
SNR = RD{|x|}/sigma donde RD es el rango dinámico.
Parameters
----------
x : 1darray
Señal a calcular SNR.
sigma : float
Desviación estandar del ruido. Por defecto es None
inicio : int, opcional
index desde el cual se calcula la desviación estándar. Por defecto es 0.
fin : int, opcional
index hasta el cual se calcula la desviación estándar. Por defecto es 100.
Returns
-------
SNR: float
Relación señal a ruido de x.
'''
if sigma==None:
sigma = np.std(x[inicio:fin])
return np.around((np.max(abs(x)) - np.min(abs(x))) / sigma, 3)
def polosyceros(b,a, lim=1.1):
'''
Grafica polos y ceros.
Devuelve un gráfico de polos y ceros y sus valores a partir de los
coeficientes de la función de transferencia. Los coeficientes deben
estar ordenados en potencias decrecientes de z.
Parameters
----------
b : array-like, list
Coeficientes del numerador de la función de transferencia H(z) ordenados
en potencias decrecientes de z.
a : array-like, list
Coeficientes del denominador de la función de transferencia H(z) ordenados
en potencias decrecientes de z.
lim : float, opcional
límites del gráfico. Por defecto es 1.1. El gráfico tendrá las mismas
proporciones en ambos ejes.
Returns
-------
Figure:
Gráfico 2D
Print:
Polos y ceros.
'''
z, p, k = sig.tf2zpk(b, a)
print('Ceros:', np.around(z,5))
print('Polos:', np.around(p,5))
tita = np.linspace(0,2*np.pi, 100)
x = np.cos(tita)
y = np.sin(tita)
plt.figure(figsize=(7,7))
plt.plot(x, y, '--k')
for i in z:
plt.scatter(i.real,i.imag, c='b', marker='o')
for i in p:
plt.scatter(i.real, i.imag, c='r', marker='x')
plt.grid(linestyle='--')
plt.xlabel('Re{z}')
plt.ylabel('Im{z}')
plt.title('Diagrama de Polos y Ceros')
plt.xlim([-lim,lim])
plt.ylim([-lim,lim])
plt.show()
def suavizado(f, amp, octava):
"""
Realiza un filtrado por bandas.
Parameters
----------
f : numpy array
eje de frecuencias
amp : numpy array
amplitudes en dB
octava : int
fracción de octava en la cual se realiza el filtrado. p.ej. octava = 3
resulta en un filtrado de 1/3 de octava. Si es 0, no hay filtrado
Returns
-------
ampsmooth : numpy array
señal suavizada
"""
ampsmooth = amp
if octava != 0: # Posibilidad de no filtrar con octava = 0
for n in range(f.size):
finf = f[n] * 2 ** (-1/(2*octava)) # Calcula frecuencia superior
fsup = f[n] * 2 ** (1/(2*octava)) # Calcula frecuencia inferior
if finf <= f[0]:
idxinf = 0
else:
idxinf = np.argmin(np.abs(f[:n+1] - finf)) # índice de la frecuencia inferior
if fsup >= f[-1]:
idxsup = f.size - 1
else:
idxsup = np.argmin(np.abs(f[:] - fsup)) # índice de la frecuencia superior
temp = 10 ** (0.1 * amp[idxinf:idxsup+1]) # Suma las presiones en la banda
ampsmooth[n] = 10 * np.log10(sum(temp) / len(amp[idxinf:idxsup+1])) # Promedio ponderado
return ampsmooth
def plot_impz(b, a = 1, l=100):
"""
Grafica la respuesta al impulso de una función de transferencia.
Parameters
----------
b : numpy array o int
numerador de la función de transferencia en potencias decrecuentes
negativas de z
a : numpy array o int, optional
denominador de la función transferencia. The default is 1.
l : int, optional
cantidad de muestras del impulso. The default is 100.
Returns
-------
None.
"""
if type(a)== int: #FIR
l = len(b)
else: # IIR
l = 500
impulse = np.repeat(0.,l)
impulse[0] =1.
x = np.arange(0,l)
response = sig.lfilter(b, a, impulse)
plt.plot(x, response); plt.grid()
plt.ylabel('Amplitud')
plt.xlabel('n (muestras)')
plt.title('Respuesta al Impulso')
def plot_magyfas(w, H, l=15, a=5):
"""
Grafica magnitud y fase de una fft.
Parameters
----------
w : numpy array
eje de frecuencias.
H : numpy array
transformada de Fourier de una señal.
l : int, optional
ancho de la figura. The default is 15.
a : int, optional
alto de la figura. The default is 5.
Returns
-------
None.
"""
H_mag = np.abs(H)
H_fase = np.angle(H)
plt.figure(1, figsize=(l,a))
plt.plot(w, H_mag, 'r')
plt.title('Magnitud de H(w)')
plt.xlabel('Frecuencia (rad/s)')
plt.ylabel('Magnitud')
plt.grid()
plt.show()
plt.figure(1, figsize=(l,a))
plt.plot(w, H_fase, 'g')
plt.title('Fase de H(w)')
plt.xlabel('Frecuencia (rad/s)')
plt.ylabel('Fase (rad/s)')
plt.grid()
plt.show()