forked from baranovmv/RoomResponse
-
Notifications
You must be signed in to change notification settings - Fork 1
/
linear_model.py
39 lines (34 loc) · 1.15 KB
/
linear_model.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
#!env python3
import numpy as np
import wavio
import scipy.signal
import scipy.fftpack
import math
import sys
from scipy.fftpack import fft, ifft
# Generate the best linear approximation (BLA) of the system
# Note: The Fourier Transform IS the best least squares fit!
# Credit to Pascal Brunet:
# Nonlinear System Modeling and Identification of Loudspeakers
# https://repository.library.northeastern.edu/files/neu:336724/fulltext.pdf
def BLA(u, y, M, P, N, w):
k = np.round(w/2/np.pi*N) + 1;
K = len(k)
Grz = np.zeros((K,M))
Var_ng = np.zeros((K,M))
realiz_u = np.reshape(u, ((*P+1* N), M))
realiz_y = np.reshape(y, ((*P + 1 * N), M))
for m in range(M):
frames_u = np.reshape(realiz_u[N+1:, m], N, P)
U = fft(frames_u)
avgU = np.mean(U, axis=1)
frames_y = np.reshape(realiz_y[N+1:,m], N, P)
Y = fft(frames_y)
avgY = np.mean(Y, axis=1)
varY = np.var(Y, axis=1)
Grz[:, m] = avgY[k] / avgU[k]
Var_ng[:,m] = 1 / P*varY[k] / (np.abs(avgU[k]) ** 2)
Gbla = np.mean(Grz,2)
VarN = 1/M*np.mean(Var_ng, axis=1)
VarBla = 1/M * np.var(Grz, axis=1)
return Gbla, VarN, VarBla