From 0b9a7c7b6c3bc7aebc62101fe2e435f265736d50 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Wed, 3 Nov 2021 16:12:24 -0400 Subject: [PATCH 01/13] Add example 10.5 --- ch10.r | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 ch10.r diff --git a/ch10.r b/ch10.r new file mode 100644 index 0000000..bbdec77 --- /dev/null +++ b/ch10.r @@ -0,0 +1,17 @@ +############# +# Chapter 10.2 Mean and correlation functions + +## Example 10.5 + +# R code for Example 10.5 +x = matrix(0, 1000, 20) +t = seq(-2, 2, (2+2)/999) + +for (i in 1:20) { + x[,i] = runif(1) * cos(2 * pi * t) +} + +matplot(t, x, lwd=2, col="gray") +lines(t, 0.5*cos(2*pi*t), lwd=4, col="darkred") +grid() + From 0a36938c516fdd45c6d46563e3ea70c6381d45f7 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Wed, 3 Nov 2021 16:17:54 -0400 Subject: [PATCH 02/13] Add example 10.6 --- ch10.r | 15 ++++++++++++++- 1 file changed, 14 insertions(+), 1 deletion(-) diff --git a/ch10.r b/ch10.r index bbdec77..365ffea 100644 --- a/ch10.r +++ b/ch10.r @@ -12,6 +12,19 @@ for (i in 1:20) { } matplot(t, x, lwd=2, col="gray") -lines(t, 0.5*cos(2*pi*t), lwd=4, col="darkred") grid() +lines(t, 0.5*cos(2*pi*t), lwd=4, col="darkred") + +## Example 10.6 + +# R code for Example 10.6 +x = matrix(0, 1000, 20) +t = seq(-2, 2, (2+2)/999) +for (i in 1:20) { + x[,i] = cos(2*pi*t+2*pi*runif(1)) +} + +matplot(t, x, lwd=2, col="gray") +grid() +lines(t, 0*cos(2*pi*t), lwd=4, col="darkred") From f3d92c34c0c9acbe12d9ee6b30053ade41610305 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Wed, 10 Nov 2021 10:27:50 -0500 Subject: [PATCH 03/13] Add example 10.7 --- ch10.r | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/ch10.r b/ch10.r index 365ffea..16e0fa9 100644 --- a/ch10.r +++ b/ch10.r @@ -28,3 +28,35 @@ for (i in 1:20) { matplot(t, x, lwd=2, col="gray") grid() lines(t, 0*cos(2*pi*t), lwd=4, col="darkred") + +## Example 10.8 + +# R code for Example 10.7 +x = matrix(0, 21, 20) +t = 0:20 + +for (i in 1:20) { + x[,i] = runif(1) ^ t +} + +stem <- function(x,y,pch=1,...){ + if (missing(y)){ + y = x + x = 1:length(x) + } + for (i in 1:length(x)){ + lines(c(x[i],x[i]), c(0,y[i]), ...) + } + lines(c(x[1]-2,x[length(x)]+2), c(0,0), ...) +} + +matplot(t, x, pch=1, col="gray", lwd=2) +grid() + +for (i in 1:ncol(x)) { + stem(t, x[,i], col="gray", lwd=2) +} +stem(t, 1/(t+1), col="darkred", lwd=2, pch=1, type="o") + + + From fd7215c8e39b1cb983821c0a96bb860c9e1cccf3 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Wed, 10 Nov 2021 16:14:55 -0500 Subject: [PATCH 04/13] Finish examples 10.11 --- ch10.r | 19 +++++++++++++++++++ 1 file changed, 19 insertions(+) diff --git a/ch10.r b/ch10.r index 16e0fa9..6763fb3 100644 --- a/ch10.r +++ b/ch10.r @@ -58,5 +58,24 @@ for (i in 1:ncol(x)) { } stem(t, 1/(t+1), col="darkred", lwd=2, pch=1, type="o") +## Example 10.11 +# R code for Example 10.11: Plot the time function +t = seq(-2, 2, len = 1000) +x = matrix(0, 1000, 20) + +for (i in 1:20) { + x[,i] = runif(1) * cos(2*pi*t) +} + +matplot(t, x, lwd=2, col="gray", type="l", lty="solid", xaxp=c(-2,2, 8)) +grid() +lines(t, 0.5*cos(2*pi*t), lwd=5, col="darkred") +points(numeric(20), x[501,], lwd=2, col="darkorange") +points(numeric(20) + 0.5, x[626,], lwd=2, col="blue", pch=4) + +# R code for Example 10.11: Plot the autocorrelation function +t = seq(-1, 1, len=1000) +R = (1/3)*outer(cos(2*pi*t), cos(2*pi*t)) +image(t, t, R, col=topo.colors(255), xlab="t_1", ylab="t_2") From 1ee487b26bac78c935f668b38ecf25eb265690f7 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Thu, 11 Nov 2021 13:32:24 -0500 Subject: [PATCH 05/13] Finish example 10.14 --- ch10.r | 49 +++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) diff --git a/ch10.r b/ch10.r index 6763fb3..d83b5db 100644 --- a/ch10.r +++ b/ch10.r @@ -79,3 +79,52 @@ t = seq(-1, 1, len=1000) R = (1/3)*outer(cos(2*pi*t), cos(2*pi*t)) image(t, t, R, col=topo.colors(255), xlab="t_1", ylab="t_2") +## Example 10.12 + +# R code for Example 10.11: Plot the time function +t = seq(-2, 2, len = 1000) +x = matrix(0, 1000, 20) + +for (i in 1:20) { + x[,i] = cos((2*pi*t) + (2*pi*runif(1))) +} + +matplot(t, x, lwd=2, col="gray", type="l", lty="solid", xaxp=c(-2,2, 8)) +grid() +lines(t, 0*cos(2*pi*t), lwd=5, col="darkred") +points(numeric(20), x[501,], lwd=2, col="darkorange") +points(numeric(20) + 0.5, x[626,], lwd=2, col="blue", pch=4) + +# R code for Example 10.12: Plot the autocorrelation function +t = seq(-1, 1, len=1000) +R = toeplitz(0.5*(cos(2*pi*t))) +image(t, t, R, col=topo.colors(255), ylim=c(1,-1), xlab="t_1", ylab="t_2") + +############# +# Chapter 10.3 Wide sense stationary processes + +## Example 10.14 + +# R code to demonstrate autocorrelation + +# Figure 1 +Xa = rnorm(1000) +Xa2 = rnorm(1000) + +plot(Xa, type="l", lwd=2, col="blue") +grid() +lines(Xa2, lwd=2) + +# Figure 2 +N = 1000 +T = 1000 +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 +} + +plot(xc[1,], type="l", lwd=2, col="darkblue") +lines(xc[2,], lwd=2) + From 05591a8ca0e550221703ab80bbcc2eb4a6d8db7e Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Thu, 11 Nov 2021 16:08:47 -0500 Subject: [PATCH 06/13] Complete chapter 10.5 --- ch10.r | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/ch10.r b/ch10.r index d83b5db..bf15c6d 100644 --- a/ch10.r +++ b/ch10.r @@ -128,3 +128,32 @@ for (i in 1:N) { plot(xc[1,], type="l", lwd=2, col="darkblue") lines(xc[2,], lwd=2) +############# +# Chapter 10.5 Wide sense stationary processes + +## Example 10.15 + +# R code for Example 10.15 +t = seq(-10, 10, by=0.001) +L = length(t) +X = rnorm(L) +h = 10 * sapply((1 - abs(t)), max, 0) / 1000 +Y = convolve(X, h, type=c("circular")) + +# Figure 1 +plot(t, X, lwd=1, col="gray", type="l") +grid() +legend(6, 3.8, legend=c("X(t)", "μ_x(t)", "Y(t)", "μ_y(t)"), col=c("gray", "black", "darkorange", "yellow"), lty=c(1, 1, 1, 3), lwd=c(1, 4, 3, 4), bg="white") +abline(h=0, lwd=4, lty=1, col="yellow") +lines(t, Y, lwd=3, col="darkorange") +abline(h=0, lwd=4, lty=3) + +# Figure 2 +h2 = convolve(h, h, type="open") +Rx = numeric(40001) +Rx[20001] = 0.2 +plot(seq(-20, 20, by=0.001), Rx, lwd=2, col="gray", type="l", xlim=c(-2,2), ylim=c(-0.05, 0.2)) +grid() +legend(-2, 0.19, legend=c("R_x(t)", "R_y(t)"), col=c("gray", "darkorange"), lty=1:1, lwd=2) +lines(seq(-20, 20, by=0.001), h2, lwd=2, col="darkorange") + From 1944f22d86d4c8525f1aa53c34b3cb65d2cf5e9f Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Thu, 11 Nov 2021 16:11:20 -0500 Subject: [PATCH 07/13] Add chapter 10.6 data --- ch10_LPC_data.txt | 320 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 ch10_LPC_data.txt diff --git a/ch10_LPC_data.txt b/ch10_LPC_data.txt new file mode 100644 index 0000000..8d71513 --- /dev/null +++ b/ch10_LPC_data.txt @@ -0,0 +1,320 @@ + 0 + 0 + 0 + 0 + 0.0000 + 0 + 0 + -0.0087 + 0 + -0.0087 + 0 + -0.0087 + -0.0087 + -0.0087 + -0.0087 + -0.0087 + 0 + 0 + -0.0087 + -0.0087 + 0.0000 + 0 + -0.0087 + -0.0087 + -0.0174 + -0.0087 + -0.0087 + -0.0087 + -0.0174 + -0.0174 + -0.0174 + -0.0261 + -0.0087 + -0.0087 + -0.0174 + -0.0174 + -0.0087 + -0.0174 + -0.0174 + 0 + -0.0087 + -0.0087 + -0.0174 + -0.0261 + -0.0261 + -0.0087 + 0 + -0.0087 + -0.0174 + -0.0174 + -0.0087 + -0.0087 + -0.0174 + -0.0174 + -0.0087 + -0.0174 + -0.0174 + -0.0087 + 0 + 0 + -0.0087 + 0 + -0.0087 + -0.0174 + -0.0261 + -0.0261 + -0.0174 + -0.0087 + 0 + -0.0087 + -0.0174 + -0.0087 + 0 + -0.0087 + -0.0087 + 0 + 0.0087 + 0.0261 + 0.0174 + 0.0174 + 0.0174 + 0.0261 + 0.0174 + 0.0261 + 0.0174 + 0.0087 + -0.0087 + -0.0348 + -0.0435 + -0.0348 + -0.0348 + -0.0435 + -0.0348 + -0.0261 + -0.0348 + -0.0261 + -0.0261 + -0.0348 + -0.0261 + -0.0348 + -0.0348 + -0.0174 + 0.0087 + 0.0174 + 0.0348 + 0.0435 + 0.0435 + 0.0435 + 0.0522 + 0.0435 + 0.0348 + 0.0348 + 0.0174 + 0.0261 + 0.0174 + 0.0087 + 0.0174 + 0.0087 + 0 + -0.0174 + -0.0261 + -0.0348 + -0.0435 + -0.0348 + -0.0435 + -0.0348 + -0.0348 + -0.0348 + -0.0348 + -0.0435 + -0.0348 + -0.0435 + -0.0522 + -0.0435 + -0.0348 + -0.0348 + -0.0261 + -0.0174 + -0.0261 + -0.0174 + -0.0261 + -0.0261 + -0.0174 + -0.0261 + -0.0261 + -0.0174 + -0.0087 + -0.0087 + 0 + 0.0087 + 0.0174 + 0.0174 + 0.0087 + 0.0261 + 0.0261 + 0.0174 + 0.0174 + 0.0261 + 0.0261 + 0.0174 + 0.0261 + 0.0174 + 0.0174 + 0.0174 + 0.0087 + 0.0087 + 0 + 0 + -0.0087 + -0.0087 + -0.0174 + -0.0174 + -0.0261 + -0.0261 + -0.0261 + -0.0174 + -0.0174 + -0.0261 + -0.0174 + -0.0174 + -0.0087 + 0 + -0.0087 + -0.0174 + -0.0087 + -0.0174 + -0.0087 + 0.0087 + 0 + 0.0000 + 0 + 0.0087 + 0.0174 + 0.0174 + 0.0261 + 0.0174 + 0.0174 + 0.0087 + 0 + 0.0087 + 0.0087 + 0.0261 + 0.0522 + 0.0609 + 0.0435 + 0.0174 + -0.0087 + -0.0522 + -0.0783 + -0.0870 + -0.0696 + -0.0609 + -0.0522 + -0.0348 + -0.0435 + -0.0696 + -0.0783 + -0.0696 + -0.0783 + -0.0783 + -0.0609 + -0.0261 + 0 + 0.0174 + 0.0435 + 0.0609 + 0.0522 + 0.0435 + 0.0261 + 0.0174 + 0.0261 + 0.0435 + 0.0522 + 0.0609 + 0.0522 + 0.0435 + 0.0348 + 0.0087 + -0.0174 + -0.0261 + -0.0348 + -0.0522 + -0.0522 + -0.0435 + -0.0522 + -0.0522 + -0.0609 + -0.0696 + -0.0783 + -0.0783 + -0.0696 + -0.0609 + -0.0435 + -0.0174 + -0.0087 + -0.0174 + -0.0174 + -0.0261 + -0.0261 + -0.0087 + -0.0087 + -0.0174 + -0.0174 + -0.0087 + 0.0087 + 0.0174 + 0.0087 + 0 + -0.0087 + -0.0087 + 0 + 0 + 0.0087 + 0.0087 + 0.0174 + 0.0174 + 0.0087 + 0.0174 + 0.0087 + 0.0087 + 0.0174 + 0.0087 + 0.0174 + 0.0174 + 0.0261 + 0.0261 + 0.0174 + 0.0087 + -0.0087 + -0.0174 + -0.0261 + -0.0261 + -0.0348 + -0.0261 + -0.0174 + -0.0261 + -0.0261 + -0.0261 + -0.0261 + -0.0261 + -0.0174 + 0 + 0 + 0.0087 + 0 + 0.0087 + 0.0348 + 0.0348 + 0.0348 + 0.0348 + 0.0435 + 0.0957 + 0.1304 + 0.1217 + 0.0783 + 0.0348 + -0.0348 + -0.1130 + -0.1478 + -0.1043 \ No newline at end of file From 6a737a531e6637ba91f4425c96fd7a365fb8458c Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Thu, 11 Nov 2021 16:25:51 -0500 Subject: [PATCH 08/13] Solve Yule Walker equation --- ch10.r | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) diff --git a/ch10.r b/ch10.r index bf15c6d..24a106f 100644 --- a/ch10.r +++ b/ch10.r @@ -157,3 +157,28 @@ grid() legend(-2, 0.19, legend=c("R_x(t)", "R_y(t)"), col=c("gray", "darkorange"), lty=1:1, lwd=2) lines(seq(-20, 20, by=0.001), h2, lwd=2, col="darkorange") +############# +# Chapter 10.6 Optimal linear filter + +## Solve the Yule Walker equation + +# R code to solve the Yule Walker Equation + +y = scan("./ch10_LPC_data.txt") +K = 10 +N = 320 +y_corr = ccf(y, y, lag.max=2*length(y)-1)$acf +R = toeplitz(y_corr[N + 1:K-1]) +lhs = y_corr[N + 1:K] +h = solve(R, lhs) + +# Figure 1 +plot(y, lwd = 4, col="blue", type="l") +grid() +legend(10, 0.1, legend=c("Y[n]"), col=c("blue"), lty=1:1, lwd=4) + +# Figure 2 +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) + From 39b1cd51eeb5f217a72d21459c4fb5063f8dd398 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Tue, 16 Nov 2021 17:35:32 -0500 Subject: [PATCH 09/13] Add LPC data part 2 --- ch10_LPC_data_02.txt | 320 +++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 ch10_LPC_data_02.txt diff --git a/ch10_LPC_data_02.txt b/ch10_LPC_data_02.txt new file mode 100644 index 0000000..a57ddf4 --- /dev/null +++ b/ch10_LPC_data_02.txt @@ -0,0 +1,320 @@ +-0.2222 +-0.1481 +-0.0370 +-0.1111 +-0.3333 +-0.5556 +-0.5185 +-0.4444 +-0.3704 +-0.2593 +-0.0370 +0.1481 +0.3333 +0.4444 +0.4815 +0.4074 +0.2593 +0.0370 +0.0741 +0.1852 +0.2593 +0.4074 +0.5556 +0.3704 +0.0741 +-0.1481 +-0.2963 +-0.3704 +-0.4074 +-0.2593 +-0.1852 +-0.2222 +-0.2222 +-0.1852 +-0.2222 +-0.2963 +-0.3333 +-0.3704 +-0.3704 +-0.2963 +-0.1852 +0 +0.0741 +0.0741 +0 +-0.2222 +-0.4074 +-0.4074 +-0.3333 +-0.1852 +-0.0741 +-0.0370 +-0.0741 +-0.1852 +-0.2593 +-0.2593 +-0.1852 +-0.1481 +-0.0370 +0 +0.0741 +0.0741 +0.0370 +0.0741 +0.0370 +0.0741 +0.2222 +0.2593 +0.3333 +0.3704 +0.4074 +0.3704 +0.2593 +0.1481 +0.1111 +0.0741 +0 +0.0741 +0.0741 +0.0370 +0 +-0.0741 +-0.1111 +-0.1481 +-0.1852 +-0.1111 +-0.0741 +-0.0370 +0.0741 +0.1111 +-0.0370 +-0.1111 +-0.1111 +-0.1852 +-0.2222 +-0.1852 +-0.1111 +-0.0741 +0.0370 +-0.0370 +-0.0370 +-0.0370 +-0.1481 +-0.1852 +-0.1481 +-0.1111 +-0.0370 +0.0741 +0.1481 +0.1111 +0.1111 +0.0741 +0.0370 +-0.0370 +0 +0.1481 +0.1481 +0.1852 +0.1111 +0.0370 +-0.0370 +-0.1481 +-0.0370 +0 +0.1852 +0.4074 +0.4074 +0.2593 +0 +-0.2963 +-0.5185 +-0.5926 +-0.4074 +-0.1111 +0.0370 +-0.0370 +-0.2222 +-0.5185 +-0.7037 +-0.7037 +-0.4815 +-0.2593 +-0.0370 +0.0741 +0.0741 +0.1111 +0.0741 +0.0741 +0.1481 +0.2222 +0.1481 +0.2593 +0.4074 +0.5185 +0.4444 +0.4074 +0.2222 +-0.0741 +-0.2222 +-0.2593 +-0.1111 +0.0370 +0.1111 +0.1852 +0.0370 +-0.2222 +-0.4815 +-0.5185 +-0.4074 +-0.1852 +-0.0370 +0.0370 +0 +-0.0741 +-0.1481 +-0.2593 +-0.2593 +-0.2222 +-0.1111 +-0.0741 +-0.0741 +0 +-0.0370 +-0.1481 +-0.2593 +-0.3333 +-0.2593 +-0.2593 +-0.1111 +0.0741 +0.0741 +-0.0370 +-0.0741 +-0.1111 +-0.0741 +-0.0370 +0.1111 +0.2963 +0.2963 +0.2593 +0.1852 +0 +-0.1111 +-0.0370 +0.0741 +0.1111 +0.1111 +0.0741 +0 +-0.0370 +-0.1481 +-0.2222 +-0.1852 +-0.1852 +-0.1111 +-0.0370 +0.0741 +0.0370 +-0.0741 +-0.1111 +-0.0741 +-0.0741 +0.0370 +0.1852 +0.1852 +0.0741 +0.0370 +0 +0 +0.0741 +0.1111 +0.2222 +0.1852 +0.1111 +0.1111 +0 +0.0370 +0.0741 +0.2593 +0.6296 +0.7778 +0.5185 +0.1111 +-0.2963 +-0.7037 +-0.7407 +-0.3704 +-0.0741 +0.0741 +0 +-0.3333 +-0.7407 +-1.0000 +-0.9259 +-0.5556 +-0.1111 +0.1111 +0.2222 +0.0370 +-0.1111 +-0.1852 +-0.0741 +0.2222 +0.4074 +0.5556 +0.5556 +0.4444 +0.3704 +0.2222 +0.2222 +0.1481 +0.0741 +-0.0370 +0 +0.0741 +0.0370 +-0.0741 +-0.2593 +-0.4815 +-0.5926 +-0.5556 +-0.3333 +-0.0741 +0.0741 +0 +-0.2593 +-0.4444 +-0.5926 +-0.4815 +-0.1852 +0.1111 +0.1852 +0.0370 +-0.1852 +-0.5185 +-0.6667 +-0.5926 +-0.3704 +-0.0741 +0.0370 +-0.1111 +-0.2963 +-0.4444 +-0.3704 +-0.2593 +-0.1111 +0.1481 +0.2593 +0.3704 +0.2963 +0.2963 +0.3333 +0.2963 +0.2963 +0.3704 +0.4444 +0.4815 +0.4815 +0.3704 +0.2593 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 10/13] 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() + From 3374e1cb014132ccfd55b2deff0ac9563cb0c85b Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Tue, 16 Nov 2021 20:36:44 -0500 Subject: [PATCH 11/13] Add wiener data --- ch10_wiener_deblur_data.txt | 320 ++++++++++++++++++++++++++++++++++++ 1 file changed, 320 insertions(+) create mode 100644 ch10_wiener_deblur_data.txt diff --git a/ch10_wiener_deblur_data.txt b/ch10_wiener_deblur_data.txt new file mode 100644 index 0000000..a9ecd78 --- /dev/null +++ b/ch10_wiener_deblur_data.txt @@ -0,0 +1,320 @@ + -0.0652 + -0.0554 + -0.0422 + -0.0272 + -0.0123 + 0.0006 + 0.0102 + 0.0160 + 0.0183 + 0.0183 + 0.0182 + 0.0203 + 0.0273 + 0.0410 + 0.0628 + 0.0924 + 0.1281 + 0.1666 + 0.2038 + 0.2351 + 0.2563 + 0.2648 + 0.2601 + 0.2435 + 0.2185 + 0.1897 + 0.1617 + 0.1386 + 0.1230 + 0.1158 + 0.1166 + 0.1234 + 0.1340 + 0.1458 + 0.1565 + 0.1644 + 0.1685 + 0.1686 + 0.1649 + 0.1583 + 0.1503 + 0.1425 + 0.1363 + 0.1329 + 0.1325 + 0.1347 + 0.1378 + 0.1398 + 0.1383 + 0.1317 + 0.1193 + 0.1022 + 0.0834 + 0.0670 + 0.0578 + 0.0599 + 0.0761 + 0.1065 + 0.1485 + 0.1965 + 0.2432 + 0.2802 + 0.2995 + 0.2950 + 0.2632 + 0.2042 + 0.1211 + 0.0199 + -0.0917 + -0.2055 + -0.3138 + -0.4105 + -0.4911 + -0.5536 + -0.5975 + -0.6239 + -0.6347 + -0.6325 + -0.6199 + -0.5991 + -0.5719 + -0.5395 + -0.5023 + -0.4602 + -0.4122 + -0.3572 + -0.2942 + -0.2224 + -0.1419 + -0.0539 + 0.0389 + 0.1328 + 0.2233 + 0.3052 + 0.3742 + 0.4267 + 0.4608 + 0.4766 + 0.4758 + 0.4619 + 0.4391 + 0.4124 + 0.3861 + 0.3637 + 0.3475 + 0.3382 + 0.3351 + 0.3362 + 0.3387 + 0.3398 + 0.3365 + 0.3267 + 0.3092 + 0.2843 + 0.2533 + 0.2189 + 0.1848 + 0.1548 + 0.1324 + 0.1200 + 0.1186 + 0.1271 + 0.1429 + 0.1624 + 0.1816 + 0.1972 + 0.2068 + 0.2099 + 0.2074 + 0.2013 + 0.1940 + 0.1878 + 0.1846 + 0.1851 + 0.1892 + 0.1960 + 0.2045 + 0.2132 + 0.2209 + 0.2260 + 0.2272 + 0.2228 + 0.2110 + 0.1903 + 0.1598 + 0.1199 + 0.0721 + 0.0197 + -0.0329 + -0.0811 + -0.1206 + -0.1485 + -0.1641 + -0.1685 + -0.1647 + -0.1569 + -0.1493 + -0.1454 + -0.1476 + -0.1564 + -0.1711 + -0.1898 + -0.2100 + -0.2295 + -0.2465 + -0.2599 + -0.2692 + -0.2743 + -0.2755 + -0.2730 + -0.2670 + -0.2575 + -0.2448 + -0.2293 + -0.2116 + -0.1928 + -0.1741 + -0.1564 + -0.1403 + -0.1259 + -0.1126 + -0.0997 + -0.0860 + -0.0706 + -0.0532 + -0.0338 + -0.0133 + 0.0071 + 0.0260 + 0.0419 + 0.0536 + 0.0601 + 0.0612 + 0.0567 + 0.0470 + 0.0325 + 0.0137 + -0.0090 + -0.0351 + -0.0643 + -0.0961 + -0.1296 + -0.1639 + -0.1972 + -0.2278 + -0.2537 + -0.2727 + -0.2832 + -0.2840 + -0.2741 + -0.2534 + -0.2222 + -0.1812 + -0.1317 + -0.0756 + -0.0152 + 0.0463 + 0.1055 + 0.1591 + 0.2042 + 0.2391 + 0.2632 + 0.2775 + 0.2840 + 0.2853 + 0.2838 + 0.2810 + 0.2776 + 0.2729 + 0.2655 + 0.2534 + 0.2353 + 0.2103 + 0.1785 + 0.1412 + 0.1001 + 0.0575 + 0.0153 + -0.0247 + -0.0618 + -0.0955 + -0.1258 + -0.1526 + -0.1755 + -0.1939 + -0.2065 + -0.2119 + -0.2090 + -0.1971 + -0.1763 + -0.1478 + -0.1138 + -0.0772 + -0.0413 + -0.0098 + 0.0144 + 0.0287 + 0.0314 + 0.0221 + 0.0009 + -0.0308 + -0.0710 + -0.1172 + -0.1664 + -0.2155 + -0.2612 + -0.3006 + -0.3308 + -0.3494 + -0.3548 + -0.3460 + -0.3234 + -0.2883 + -0.2436 + -0.1931 + -0.1413 + -0.0928 + -0.0516 + -0.0203 + -0.0003 + 0.0087 + 0.0077 + -0.0013 + -0.0165 + -0.0360 + -0.0585 + -0.0824 + -0.1061 + -0.1274 + -0.1439 + -0.1532 + -0.1531 + -0.1429 + -0.1228 + -0.0948 + -0.0619 + -0.0280 + 0.0035 + 0.0296 + 0.0488 + 0.0612 + 0.0681 + 0.0717 + 0.0748 + 0.0795 + 0.0872 + 0.0983 + 0.1119 + 0.1260 + 0.1384 + 0.1467 + 0.1492 + 0.1451 + 0.1346 + 0.1187 + 0.0993 + 0.0784 + 0.0583 + 0.0403 + 0.0256 \ No newline at end of file From 5d86e7f9b319e8b4a581bbe953846200b8cc04c7 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Tue, 16 Nov 2021 20:36:49 -0500 Subject: [PATCH 12/13] Finish wiener deblurring --- ch10.r | 55 ++++++++++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 48 insertions(+), 7 deletions(-) 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) + +############# + From fcbd90dac9448fce034a44065628d58c41edc9a3 Mon Sep 17 00:00:00 2001 From: aaaakshat <33050725+aaaakshat@users.noreply.github.com> Date: Tue, 16 Nov 2021 20:37:22 -0500 Subject: [PATCH 13/13] Rename ch07.R -> ch07.r --- ch07.R => ch07.r | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename ch07.R => ch07.r (100%) diff --git a/ch07.R b/ch07.r similarity index 100% rename from ch07.R rename to ch07.r