ウィナー・ホッフ方程式をPythonで解いてみた

 #ウィナーフィルタを最小二乗法で解く

{\boldsymbol d}=[d_1, \cdots, d_{L_s}]

{\boldsymbol U}=[ {\boldsymbol u}_1, \cdots,  {\boldsymbol u} _ {L _ s}  ] =

\begin{bmatrix} u _ 1 & u _ 2 & \cdots & u _ K & \cdots & u _ {L _ s - 1} & u _ {L _ s} \\ 0 & u _ 1 & \cdots & u _ {K - 1} & \cdots & u _ {L _ s - 2} & u _ {L _ s - 1} \\ \vdots & & \ddots & \vdots & & \vdots & \vdots \\ 0 & 0 & \cdots & u _ 1 & \cdots & u _ {L _ s - K} & u _ {L _ s - K + 1} \\ \end{bmatrix}

\hat{{\boldsymbol R}} _ u = {\boldsymbol U}{\boldsymbol U}^H

\hat{{\boldsymbol r}} _ {ud} = {\boldsymbol U}{\boldsymbol d}^H

\hat{{\boldsymbol \omega}} _ {LS} =  \hat{{\boldsymbol R}_u} ^ {-1}  \hat{{\boldsymbol r}} _ {ud}

# -*- coding: utf-8 -*-
"""
Created on Sat Jun 22 22:14:03 2019

@author: user
"""

import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sg


z  = np.arange(0, 512)
h0 = 1 / (1 + np.exp(0.0001 * (z - 128) ** 2))
#h0=1/(1+np.exp(0.03*z))

Wref  = h0 / np.sum(h0) * 2
plt.plot(Wref)
plt.title("$\omega_{ref}$")
plt.show()

x  = np.linspace(0, 1, 4096)
#y1 = np.sin(2 * np.pi * x) + 0.2 * np.sin(np.pi * 20 * (x + 0.1))
y1 = np.sin(2 * np.pi * x) + 0.2*np.sin(np.pi * 20 * x + 0.1 * np.pi)
y2 = sg.lfilter(Wref, [1], y1)

plt.plot(x, y1)
plt.plot(x, y2)
plt.legend(["$y1$", "$y2=y1 * \omega_{ref}$"])
plt.show()


Ls  = 1024
K   = 512
U   = np.array([[ y1[l - k] if l >= k else 0 for l in range(0, Ls)] for k in range(0, K)])
Ru  = U @ U.T
d   = y2[:Ls]
rud = U @ d.T
Wls  = np.linalg.solve(Ru, rud)

plt.plot(Wref)
plt.plot(Wls)
plt.legend(["$\omega_{ref}$", "$\omega_{LS}=(U U^T)^{-1} U d^T$"])
plt.show()

plt.plot(Wls - Wref);
plt.legend(["$\omega_{LS}-\omega_{ref}$"])
plt.show()

plt.plot(x, y2);
plt.plot(x, sg.lfilter(Wls, [1], y1))
plt.legend(["$y2$", "$y1 * \omega_{LS}$"])
plt.show()

plt.plot(x, y2 - sg.lfilter(Wls, [1], y1))
plt.legend(["$y2 - y1 * \omega_{LS}$"])
plt.show()