deconvolve {spectral} | R Documentation |
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
)