ウィナー・ホッフ方程式をPythonで解いてみた
#ウィナーフィルタを最小二乗法で解く
]
] =
\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}
# -*- 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()