Fix wiener filter for real value for the noise (#17)

* Fix wiener filter for real value for the noise

* Update documentation for Wiener deconvolution of noisy time series

* Add example with blurring of a time series to documentation code

* Add example of blurred, noisy time series to documentation

* Add constant noise example to documentation

* Fix runtests

* Add missing test cases to wiener

doc/src/index.md

@@ -105,11 +105,11 @@ series.
 We first construct the noisy signal:
 
 ``` {.sourceCode .julia}
-using LombScargle, Deconvolution, Plots
-t = linspace(0, 10, 1000) # observation times
-x = sinpi(t) .* cos.(5t) .- 1.5cospi.(t) .* sin.(2t) # the original signal
+using LombScargle, Deconvolution, Plots, Statistics
+t = range(0, stop=10, length=1000) # observation times
+x = sinpi.(t) .* cos.(5t) - 1.5cospi.(t) .* sin.(2t) # the original signal
 n = rand(length(x)) # noise to be added
-y = x + 3(n - mean(n)) # observed noisy signal
+y = x + 3(n .- mean(n)) # observed noisy signal
 ```
 
 In order to perform the Wiener deconvolution, we need a signal that has
@@ -146,12 +146,13 @@ wiener with a simple signal that is the sum of these three models:
 signal = m1 + m2 + m3 # signal for `wiener`
 noise = rand(length(y)) # noise for `wiener`
 polished = wiener(y, signal, noise)
-# Compare...
 plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
-plot!(t, y, label="Observed signal") # ...original and observed signal
+plot!(t, y, label="Observed signal")
+savefig("time-series-observed.png")
 plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
-plot!(t, polished, label="Recovered with Wiener") # ...original and recovered signal
-plot!(t, signal, label="Lomb–Scargle model") #...and best fitting Lomb–Scargle model
+plot!(t, polished, label="Recovered with Wiener")
+plot!(t, signal, label="Lomb–Scargle model")
+savefig("time-series-recovered.png")
 ```
 
 ![image](wiener/time-series-observed.png)
@@ -166,6 +167,44 @@ With real-world data the Lomb–Scargle periodogram may not work as good
 as in this toy-example, but we showed a possible strategy to create a
 suitable signal to use with wiener function.
 
+#### Blurred noisy time series
+Additionally to noise, also a blurring kernel can applied to the image.
+First we define the blurring kernel as a Gaussian kernel.
+Using `ifftshift` we move the center to index position 1 which is required
+for the Wiener deconvolution algorithm.
+
+``` {.sourceCode .julia}
+# Gaussian blurring kernel
+kernel = exp.( - 10 .* (t .- 5).^2)
+kernel ./= sum(kernel) # normalize kernel to sum of 1
+kernel = ifftshift(kernel) # move center to index pos 1
+```
+
+The blurring can be applied simply in Fourier space. After the blurring we
+add the noise.
+
+``` {.sourceCode .julia}
+y_blurred = real(ifft(fft(kernel) .* fft(x))) + noise
+```
+
+The deblurred image can be obtained with the following call.
+
+``` {.sourceCode .julia}
+deblurred = wiener(y_blurred, signal, noise, kernel)
+```
+
+Additionally to the `noise` array we also pass `kernel` to `wiener`.
+
+``` {.sourceCode .julia}
+plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
+plot!(t, deblurred, label="Deblurred with Wiener")
+savefig("time-series-deblurred.png")
+```
+
+![image](wiener/time-series-observed-blurred.png)
+
+![image](wiener/time-series-deblurred.png)
+
 #### Blurred image
 
 Here is an example of use of wiener function to perform the Wiener
@@ -202,20 +241,34 @@ noise2 = rand(size(img)) # Create another additive noise
 # Polish the image with Deconvolution deconvolution
 polished2 = wiener(blurred_img, img2, noise2, blurring)
 
+# Wiener also works using a real number instead of a noise array
+polished3 = wiener(blurred_img, img2, 1000, blurring)
+polished4 = wiener(blurred_img, img2, 10000, blurring)
+
+
 # Compare...
 view(img) # ...the original image
 view(blurred_img) # ...the blurred image
 view(polished) # ...the polished image
 view(polished2) # ...the second polished image
+view(polished3) # ...the third polished image
+view(polished4) # ...the fourth polished image
 ```
 
-| Original image                                     | Blurred image                                              |
-| :------------------------------------------------- | :--------------------------------------------------------- |
-| ![](wiener/original.jpg)                           | ![](wiener/blurred.jpg)                                    |
+| Original image                                            | Blurred image                                              |
+| :-------------------------------------------------        | :--------------------------------------------------------- |
+| ![](wiener/original.jpg)                                  | ![](wiener/blurred.jpg)                                    |
+
+| Image restored with exact power spectrum and noise        | Image restored with imperfect reference noise and spectrum |
+| :-------------------------------------------------        | :--------------------------------------------------------- |
+| ![](wiener/polished.jpg)                                  | ![](wiener/polished2.jpg)                                  |
+
+| Image restored with imperfect spectrum and constant noise | Image restored with imperfect spectrum and constant noise  |
+| :-------------------------------------------------------- | :--------------------------------------------------------- |
+| ![](wiener/polished3.jpg)                                 | ![](wiener/polished4.jpg)                                  |
+
 
-| Image restored with exact power spectrum and noise | Image restored with imperfect reference noise and spectrum |
-| :------------------------------------------------- | :--------------------------------------------------------- |
-| ![](wiener/polished.jpg)                           | ![](wiener/polished2.jpg)                                  |
+Without knowing the noise array exactly the contrast drops significantly. Some postprocessing of the contrast can enhance the quality further.
 
 ### Richardson-Lucy deconvolution
 

doc/src/wiener/cameraman.jl

@@ -45,8 +45,13 @@ img2 = channelview(testimage("livingroom")) # Load another image
 noise2 = rand(Float64, size(img)) # Create another additive noise
 # Polish the image with Deconvolution deconvolution
 polished2 = wiener(blurred_img, img2, noise2, blurring)
+# Wiener also works using a real number instead of a noise array
+polished3 = 0.5 .* wiener(blurred_img, img2, 1000, blurring)
+polished4 = 0.5 .* wiener(blurred_img, img2, 10000, blurring)
 
 save("original.jpg", Images.clamp01nan.(img))
 save("blurred.jpg", Images.clamp01nan.(blurred_img))
 save("polished.jpg", Images.clamp01nan.(polished))
 save("polished2.jpg", Images.clamp01nan.(polished2))
+save("polished3.jpg", Images.clamp01nan.(polished3))
+save("polished4.jpg", Images.clamp01nan.(polished4))

doc/src/wiener/time-series.jl

@@ -19,11 +19,22 @@ using Random
 
 Random.seed!(7)
 
-using LombScargle, Deconvolution, Plots, Statistics
+using LombScargle, Deconvolution, Plots, Statistics, FFTW
 t = range(0, stop=10, length=1000) # observation times
 x = sinpi.(t) .* cos.(5t) - 1.5cospi.(t) .* sin.(2t) # the original signal
+
+# Gaussian blurring kernel
+kernel = exp.( - 10 .* (t .- 5).^2)  
+kernel ./= sum(kernel) # normalize kernel to sum of 1
+kernel = ifftshift(kernel) # move center to index pos 1
+
 n = rand(length(x)) # noise to be added
-y = x + 3(n .- mean(n)) # observed noisy signal
+noise = 3 .* (n .- mean(n))
+y = x + noise # observed noisy signal
+# blurred and noise signal
+y_blurred = real(ifft(fft(kernel) .* fft(x))) + noise
+
+
 # Lomb-Scargle periodogram
 p = lombscargle(t, y, maximum_frequency=2, samples_per_peak=10)
 plot(freqpower(p)...)
@@ -33,12 +44,24 @@ m2 = LombScargle.model(t, y, findmaxfreq(p, [0.5, 1])[1]) # second model
 m3 = LombScargle.model(t, y, findmaxfreq(p, [1, 1.5])[1]) # third model
 
 signal = m1 + m2 + m3
-noise = rand(length(y)) # noise for `wiener`
 polished = wiener(y, signal, noise)
+deblurred = wiener(y_blurred, signal, noise, kernel)
+
+# Plots 
 plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
 plot!(t, y, label="Observed signal")
 savefig("time-series-observed.png")
+
+plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
+plot!(t, y_blurred, label="Blurred signal")
+savefig("time-series-observed-blurred.png")
+
 plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
 plot!(t, polished, label="Recovered with Wiener")
 plot!(t, signal, label="Lomb–Scargle model")
 savefig("time-series-recovered.png")
+
+plot(t, x, size=(900, 600), label="Original signal", linewidth=2)
+plot!(t, deblurred, label="Deblurred with Wiener")
+plot!(t, signal, label="Lomb–Scargle model")
+savefig("time-series-deblurred.png")

src/wiener.jl

@@ -52,8 +52,14 @@ function wiener(input::AbstractArray, signal::AbstractArray,
 end
 
 # Noise as a real (it's converted to an array)
-wiener(input::AbstractArray, signal::AbstractArray, noise::Real) =
-    wiener(input, signal, fill!(similar(input), one(eltype(input)))*noise)
+function wiener(input::AbstractArray, signal::AbstractArray, noise::Real)
+    @assert size(input) == size(signal)
+    input_ft = fft(input)
+    signal_power_spectrum = abs2.(fft(signal))
+    noise_power_spectrum = noise .* ones(eltype(input), size(input))
+    return _wiener_no_blur(input_ft, signal_power_spectrum,
+                           noise_power_spectrum)
+end
 
 ## With blurring
 function wiener(input::AbstractArray, signal::AbstractArray,
@@ -68,9 +74,17 @@ function wiener(input::AbstractArray, signal::AbstractArray,
 end
 
 # Noise as a real (it's converted to an array)
-wiener(input::AbstractArray, signal::AbstractArray,
-       noise::Real, blurring::AbstractArray) =
-           wiener(input, signal, fill!(similar(input), one(eltype(input)))*noise, blurring)
+function wiener(input::AbstractArray, signal::AbstractArray,
+                noise::Real, blurring::AbstractArray)
+    @assert size(input) == size(signal) == size(blurring)
+    input_ft = fft(input)
+    signal_power_spectrum = abs2.(fft(signal))
+    noise_power_spectrum = noise .* ones(eltype(input), size(input))
+    blurring_ft = fft(blurring)
+    return _wiener_with_blur(input_ft, signal_power_spectrum,
+                             noise_power_spectrum, blurring_ft)
+end
+
 
 """
     wiener(input, signal, noise[, blurring])

test/runtests.jl

@@ -8,27 +8,41 @@ Random.seed!(42) # Fixed random seed
 x = range(0, stop = 10, length = 15)
 s = sinpi.(x) .- 1.5cos.(x)
 n = rand(length(s))
-
-@testset "No blurring" begin
+blurring_ft  = exp.(-0.001 * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
+blurred_s_ft = fftshift(blurring_ft) .* fft(s) .+ fft(n)
+blurred_s    = real(ifft(blurred_s_ft))
+blurring     = ifft(fftshift(blurring_ft))
+@testset "No blurring with noise scalar parameter" begin
     @test wiener(s + n, s, 0.5) ≈
-        [-1.5300909616077885,-0.4607553796475224,-1.73313390900302,0.8555916473930066,
-         2.2692599396286277,0.7989152977213652,1.1431959877178226,-0.8118586562136341,
-         -1.6519310564150556,-0.4882369786466372,-1.0387206070662094,
-         -0.8324097927723406,1.8208048770862535,0.7427114498907227,1.0587917379150844]
+        [-1.1111473161507488, -0.17041312679112874, -1.3469825775608049, 1.2447585407259756, 
+          2.578860561637874, 1.09191362235022, 1.5339504088753226, -0.34582554585559355, 
+          -1.316232389748769, -0.21271451029702804, -0.6157543266519075, -0.3338395657941805, 
+          2.1195369797597614, 1.060270051101612, 1.4446941593151053]
 end
 
-@testset "Blurring" begin
-    blurring_ft  = exp.(-0.001 * (range(-div(length(x), 2), length = length(x)) .^ 2) .^ (5 // 6))
-    blurred_s_ft = fftshift(blurring_ft) .* fft(s) .+ fft(n)
-    blurred_s    = real(ifft(blurred_s_ft))
-    blurring     = ifft(fftshift(blurring_ft))
+@testset "Blurring with noise scalar parameter" begin
     @test wiener(blurred_s, s, 0.5, blurring) ≈
-        [-1.5288865581538547,-0.4592713125163481,-1.7360819833552192,0.8525544925861452,
-         2.272966823544607,0.8021732396525801,1.14091217864045,-0.8166101187683248,
-         -1.6455701546601862,-0.4926325667795202,-1.0323127293402468,
-         -0.8380228416032021,1.8214493186523952,0.7447075584686396,1.0566048110361044]
+        [-1.1086114591808454, -0.17138676636801406, -1.3482585663111306, 1.2427663740234156, 
+         2.5813314219359276, 1.0940885915582341, 1.532411498782748, -0.34834962964626176, 
+         -1.311080379275804, -0.21869338550577228, -0.6084828281998079, -0.3366932366271657, 
+         2.117735346098039, 1.0624233782615877, 1.4429704301669912]
+end
+
+@testset "No blurring with noise array" begin
+    @test wiener(blurred_s, s, ifftshift(n)) ≈ [-1.4080797956919306, -0.49286583995619976, 
+        -1.6205631520484833, 0.7973315264441291, 2.041487785980698, 0.7068033074341512, 
+        1.303796832927927, -0.5691518774600562, -1.7799047917607576, -0.7841432287697161, 
+        -0.9149839099122377, -0.4555342804458623, 1.761794603633166, 0.5975996848388563, 0.9542148100114622]
+end
+
+@testset "Blurring with noise array" begin
+    @test wiener(blurred_s, s, ifftshift(n), blurring) ≈ [-1.4250670444336333, -0.48221282618655026, 
+        -1.636981677852502, 0.8020690789792746, 2.053748123047759, 0.6974830808841243, 1.3177198257314031, 
+        -0.5714565338646487, -1.7922287374907637, -0.7809316185676288, -0.9185450961060387, 
+        -0.46118637814202684, 1.7785714975739182, 0.5887392863050276, 0.967947374977039]
 end
 
+
 ##### Richardson-Lucy deconvolution
 
 @testset "Richardson-Lucy deconvolution tests" begin