Finish wiener deblurring

ch10.r

@@ -122,7 +122,7 @@ X = matrix(rnorm(N*T), N, T)
 xc = matrix(0, N, 2*T-1)
 
 for (i in 1:N) {
-    xc[i,] = ccf(X[i,], X[i,], lag.max=2*T-1)$acf/T
+    xc[i,] = ccf(X[i,], X[i,], lag.max=2*T-1, pl=FALSE)$acf/T
 }
 
 plot(xc[1,], type="l", lwd=2, col="darkblue")
@@ -167,7 +167,7 @@ lines(seq(-20, 20, by=0.001), h2, lwd=2, col="darkorange")
 y = scan("./ch10_LPC_data.txt")
 K = 10
 N = 320
-y_corr = ccf(y, y, lag.max=2*length(y)-1)$acf
+y_corr = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
 R = toeplitz(y_corr[N + 1:K-1])
 lhs = y_corr[N + 1:K]
 h = solve(R, lhs)
@@ -207,35 +207,76 @@ for (t in 1:20) {
 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)
+lines(y, lwd=4, col="gray", lty=3)
 
 ## Wiener filter 
 
 # R code for Wiener filtering
 library(pracma)
+y = scan("./ch10_LPC_data_02.txt")
 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
+Ry = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
+Rw = ccf(w, w, lag.max=2*length(w)-1, pl=FALSE)$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))
+yhat = Re(ifft(as.vector(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)
+lines(y, lwd=2, lty=3)
 
 # 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()
 
+## Wiener deblurring
+
+# R code to solve the Wiener deconvolution problem
+library(pracma)
+
+conv_same = function(x, y) {
+    s = length(y) / 2
+    e = length(x) + s - 1
+    return (convolve(x, y, type="open")[s:e])
+}
+
+y = scan("./ch10_wiener_deblur_data.txt")
+g = ones(32, 1)/32
+w = 0.02 * rnorm(320)
+s = length(y)/2
+e = length(g) + s - 1
+x = conv_same(y, g) + w
+
+Ry = ccf(y, y, lag.max=2*length(y)-1, pl=FALSE)$acf
+Rw = ccf(w, w, lag.max=2*length(w)-1, pl=FALSE)$acf
+Sy= fft(Ry)
+Sw = fft(Rw)
+a = g[1:639]
+a[is.na(a)] = 0
+G = fft(a[1:639])
+
+H = (Conj(G) * Sy) / (abs(G) ^ 2 * Sy + Sw)
+b = x[1:639]
+b[is.na(b)] = 0
+Yhat = H * fft(b[1:639])
+yhat = Re(ifft(as.vector(Yhat)))
+
+plot(x, lwd=4, col="gray", type="l", ylim=c(-0.7, 0.7))
+grid()
+legend(150, -0.45, 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(16:(320+15), yhat[1:320], col="red", lwd=2)
+lines(1:320, y, lwd=2, lty=3)
+
+#############
+