From 93bff7f72c325980c83c27ea01712d78723783fb Mon Sep 17 00:00:00 2001
From: aaaakshat <33050725+aaaakshat@users.noreply.github.com>
Date: Tue, 16 Nov 2021 17:35:46 -0500
Subject: [PATCH] Finish wiener filtering

---
 ch10.r | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 57 insertions(+)

diff --git a/ch10.r b/ch10.r
index 24a106f..3fbdde1 100644
--- a/ch10.r
+++ b/ch10.r
@@ -182,3 +182,60 @@ plot(y_corr, lwd = 4, type="l")
 grid()
 legend(10, 0.9, legend=c("R_y[k]"), col=c("black"), lty=1:1, lwd=4)
 
+## Predict sample
+
+# R code to predict the samples
+y = scan("./ch10_LPC_data_02.txt")
+K = 10
+N = length(y)
+
+y_corr = ccf(y, y, lag.max=2*N-1)$acf[,,1]
+R = toeplitz(y_corr[N + 1:K-1])
+lhs = y_corr[N + 1:K]
+h = solve(R, lhs)
+
+z = y[311:320]
+yhat = numeric(340)
+yhat[1:320] = y
+
+for (t in 1:20) {
+    predict = t(z) * h
+    z = c(z[2:10], predict)
+    yhat[320+t] = predict
+}
+
+plot(yhat, lwd=3, col="red", type="l")
+grid()
+legend(10, 0.9, legend=c("Prediction", "Input"), col=c("red", "gray"), lty=c(1, 2), lwd=2)
+lines(y, lwd=4, col="gray", lty=2)
+
+## Wiener filter 
+
+# R code for Wiener filtering
+library(pracma)
+w = 0.05 * rnorm(320)
+x = y + w
+
+Ry = ccf(y, y, lag.max=2*length(y)-1)$acf
+Rw = ccf(w, w, lag.max=2*length(w)-1)$acf
+Sy = fft(Ry)
+Sw = fft(Rw)
+H = Sy / (Sy + Sw)
+
+a = x[1:639]
+a[is.na(a)] = 0
+Yhat = H * fft(a[1:639])
+yhat = Re(ifft(Yhat))
+
+# Figure 1
+plot(x, lwd=5, col="gray", type="l")
+grid()
+legend(10, -0.7, legend=c("Noisy Input X[n]", "Wiener Filtered Yhat[n]", "Ground Truth Y[n]"), col=c("gray", "red", "black"), lty=c(1, 1, 2), lwd=2)
+lines(yhat[1:320], lwd=2, col="red")
+lines(y, lwd=2, lty=2)
+
+# Figure 2
+plot(Rw, lwd=4, col="blue", label="h[n]", type="l")
+legend(500, 0.85, legend=c("h[n]"), col=c("blue"), lty=c(1), lwd=c(4))
+grid()
+