deconvolve {spectral}

Deconvolve Sampling Spectrum for Equidistant Sampling

Description

The function removes the probable alias peaks in the power spectral density. These projections originate from correlated gaps, missing values and interactions with noise. The function should be considered as *experimental* but with didactic background.

Usage

deconvolve(x, y, SNR.enable = T, SNR.level = 1)

Arguments

x

sampling instances

y

values

SNR.enable

binary value, include or exclude the noise

SNR.level

theshold in the sense of a multiple of mean() noise level

Details

In the special case of a non complete equidistant grid containing the data and missing values (NA), this function performs the deconvolution of Y = fft(y) from the sampling spectrum of the aquisition series x. The data is assumed to exist on a equidistant grid with missing values and gaps.

Given a one dimensional vector y of data this function reverses the spectral convolution of Y = S * X + N, if * describes the convolution operation and Y = F(y) denotes the discrete Fourier transform via the operator F(.). If, the sampling series x is considered to be purely deterministic, which should be the case for captured data, or the distortions (missing values, gaps) are *correlated* (see example), then there exists an analytic inversion of the convolution. Given the general definition of power spectral density |Y|^2 = |S * X + N|^2 the challenge is to prove |S * X + N|^2 ~ |S|^2 * |X|^2 + |N|^2. Here N describes a stochastic term of gaussian noise. This issue is solved in correlation space, where convolution becomes a multiplitation. The auto correlation function (acf) of y is given by Ry = F(|Y|^2). As a remark, IF we consider the special case of equispaced sampling, modeled by the Diraq distribution \delta(x), it is easy to show that the correlation function of a product is the product of individual correaltaion functions, F(|S*X|^2) = F(|S|^2) . F(|X|^2).

The aim is, to approximate S as the "true" spectrum. To the cost of the phase information, the result is the standardized power spectral density. The spectral noise term F(N) is approximated by a theshold in Fourier space. Here SNR.level sets the factor of mean(fft(y)) below which noise level is assumed. Above this value, the signal should be present. As a parameter to play with, SNR.enable enables or disables the noise term. This parameter was introduced to be consistent with present approaches, not considering the presence of noise.

Value

list of frequency f and spectral density function S

Examples

### Deconvolution ###
#
#
# we define a test function with gaps and noise
# we show:
# - the aliased Fourier spectrum and for comparison Lomb Spectrum
# - the corrected spectrum
#

## definition of the sampling series
x <- seq(0,pi/2,by = 5e-3)
n <- rnorm(length(x),sd = 0.1)

## definition of the test function
## with 2 frequencies
yf <- function(x)
{
  cos(2*pi*5.123*x) +
  cos(2*pi*17*x)
}

y <- yf(x)
y <- y - mean(y)

## define strongly correlated gaps
i <- NULL

i <- c(i,which(sin(2*pi / 0.3 * x) - 0.5 > 0))
i <- c(i,which(sin(2*pi / 0.04 * x + 1.123) - 0.5 > 0))
i <- sort(unique(i))


xs <- x
ys <- yf(xs) + n # add some noise
ys[i] <- NA

## for comparison we calculate a Lomb-Spectrum
LT <- spec.lomb(x = xs,y = ys
                ,f = seq(0,250,by = 0.02)
                ,mode = "generalized"
)

WS <- deconvolve(x = xs, y = ys,SNR.enable = 1,SNR.level = 1)
FT <- spec.fft(x = xs, y = ys,center = FALSE)
FTS <- spec.fft(x = xs, y = is.na(ys),center = FALSE)
LTS <- spec.lomb(x = xs, y = is.na(ys),f = seq(0,250,by = 0.02))

### results ###
#
# - signal spectrum (solid) dominant peaks at around f0 = {5, 17}
# - (minor) alias peaks (grey line, FFT dots) at f0 +/- fs
# - sampling spectrum (dashed) with fs = {3.3, 25} (dominant modes)
# - deconvolved spectrum (solid black) rejects the aliases and sampling
#
#

### time series

par(mfrow = c(1,1),mar = c(4,4,3,0.3))
curve(yf,0,max(x), col = "grey",n = 1000
      ,xlim = c(0,max(x)),ylim = c(-2,3)
      ,xlab = "Time", ylab = "y(t)"
      ,main = "Fragmented Time Series"
      )
points(xs,ys)
points(xs[is.na(ys)],yf(xs[is.na(ys)]),pch = 16,cex = 0.5)

legend("topright",c("y(t)","y(tn) + n(tn)","NA's")
       ,lty = c(1,NA,NA)
       ,lwd = c(1,NA,NA)
       ,pch = c(NA,1,16)
       ,col = c("darkgrey","black","black")
       ,bg = "white"
       ,cex = 0.8
)

## plot spectra
par(mfrow = c(1,1),mar = c(4,4,3,0.3))
with(FT,plot(fx,PSD,type="p",log = "x"
     # ,col="grey"
     ,xlim = c(1,100),ylim = c(1e-2,0.75)
     ,xlab = "f", ylab = "PSD"
     ,pch = 1
     ,lwd = 1
     ,main = "Spectra"
))
with(LT,lines(f,PSD,col = "grey",lwd = 4))
with(WS,lines(f,S, lwd = 2, col = "black"))
with(LTS,lines(f,PSD,lty = 2))
abline(h = c(1,0.5),lty = 3)
legend("topright",c("Fourier","Lomb","Decon.","Sampling")
       ,lty = c(NA,1,1,2)
       ,lwd = c(2,2,2,2)
       ,pch = c(1,NA,NA,NA)
       ,col = c("black","grey","black","black")
       ,bg = "white"
       ,cex = 0.8
       ,ncol = 2
)