diff --git a/ch10.r b/ch10.r index 3fbdde1..1923c77 100644 --- a/ch10.r +++ b/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,17 +207,18 @@ 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) @@ -225,17 +226,57 @@ 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) + +############# +