I am trying to create a Wiener filter and check for filter values from 2 to 20. I get 18 graphs as required but I keep getting an error even though I get an output. I know my code is not the cleanest and I only started Python a few months ago as Lockdown project. Any and all help will be appreciated. The code is written according to formulae provided in the question.
import numpy as np
import scipy as sp
from scipy import linalg
from scipy import signal
import matplotlib.pyplot as plt
import scipy.io as sio
N = 109068
data = sio.loadmat('audiodata.mat')
y = data['audio']
y = np.reshape(y,(109068,1))
v2 = data['noise']
v2 = np.reshape(v2,(109068,1))
for t in range(1,18):
for n in range (t+1,20):
rv2 = np.zeros(t+1)
ryv2 = np.zeros(t+1)
rv2[t] += 1/(n-t)*v2[n]*v2[n-t]
ryv2[t] += 1/(n-t)*y[n]*v2[n-t]
Rv2 = sp.linalg.toeplitz(rv2)
Rv_2 = sp.linalg.inv(Rv2)
w_opt = np.matmul(Rv_2,ryv2)
v1 = sp.signal.lfilter(w_opt,1,v2)
y_filt = y-v1
ts = np.arange(0., y.shape[0]/22050, 1/22050)
fig, (ax1, ax2, ax3) = plt.subplots(3,sharex=True,figsize=(15,10))
ax1.plot(ts, y,'b')
ax1.plot(ts, y_filt,'g')
ax2.plot(ts, v1)
ax2.plot(ts, v2)
ax3.plot(ts, y_filt)
plt.xlabel('Time[s]')
plt.ylabel('PCM')
plt.show
The error is get is this:
---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
<ipython-input-244-7a0d629580e8> in <module>
19 ryv2[t] += 1/(n-t)*y[n]*v2[n-t]
20 Rv2 = sp.linalg.toeplitz(rv2)
---> 21 Rv_2 = sp.linalg.inv(Rv2)
22 w_opt = np.matmul(Rv_2,ryv2)
23 v1 = sp.signal.lfilter(w_opt,1,v2)
~\Anaconda3\lib\site-packages\scipy\linalg\basic.py in inv(a, overwrite_a, check_finite)
972 inv_a, info = getri(lu, piv, lwork=lwork, overwrite_lu=1)
973 if info > 0:
--> 974 raise LinAlgError("singular matrix")
975 if info < 0:
976 raise ValueError('illegal value in %d-th argument of internal '
LinAlgError: singular matrix
Any thoughts ?
I got the mistake. I just needed another for loop above t to correctly iterate for t itself.