Finish wiener filtering

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()
+