waterfall {spectral} | R Documentation |
Estimate the local frequencies
Description
A waterfall
-diagramm displays the local frequency in dependence of
or spatial vector. One can then locate an event in time or space.
Usage
waterfall(
y = stop("y value is missing"),
x = NULL,
nf = 3,
type = "b",
width = 7
)
Arguments
y |
numeric real valued data vector |
x |
numeric real valued spatial vector. (time or space) |
nf |
steepness of the bandpass filter, degree of the polynomial. |
type |
type of weightening function: "poly", "sinc", "bi-cubic","gauss", can be abbreviated |
width |
normalized maximum "inverse" width of the bandpass |
Details
Each frequency is evaluated by calculating the amplitude demodulation, which
is equivalent to the envelope function of the band pass filtered signal.
The frequency of interest defines automatically the center frequency fc
of the
applied band pass with the bandwidth BW
:
BW = fc / width, BW < width -> BW = width, BW > width -> BW = fc / width
The frequency is normalized so the minimal frequency is 1
.
With increasing frequency the bandwidth becomes wider, which lead to a variable
resolution in space and frequency. This is comparable to the wavelet
(or Gabor) transform,
which scales the wavelet (window) according to the frequency.
However, the necessary bandwidth is changed by frequency to take the
uncertainty principle into account. Slow oscillating events are measured precisely
in frequency and fast changing processes can be determined more exact in space.
This means for a signal with steady
increasing frequency the waterfall
function will produce a diagonally
stripe. See the examples below.
Value
a special fft
-object is returned. It has mode "waterfall" and
x
and fx
present, so it is only plotable.
Missing values
Given a regualar grid x_i = \delta x \cdot i
there might be missing values
marked with NA
, which are treated by the function as 0's.
This "zero-padding" leads to a loss of signal energy being
roughly proportional to the number of missing values.
The correction factor is then (1 - Nna/N)
as long as Nna / N < 0.2
.
As long as the locations of missing values are randomly
distributed the implemented procedure workes quite robust. If, in any case,
the distribution becomes correlated the proposed correction is faulty and
projects the wrong energies.
The amplitudes and PSD values are compensated to show up an estimate of the "correct" value. Therefore this method is experimental
Examples
#### noisy signal with amplitude modulation ####
x <- seq(0,3, length.out = 1000)
# original data
# extended example from envelope function
y <- 1*(abs(x-1.5))*sin(10*2*pi*x) + ifelse(x > 1.5,sin(15*(1+0.25*(x - 1.5))*2*pi*x),0)
ye <- base::Re(envelope(y))
par(mfrow=c(2,1),mar=c(1,3.5,3,3),mgp=c(2.5,1,0))
# plot results
plot(x,y,type="l",lwd=1,col="darkgrey",lty=2,ylab="y",main="Original Data",xaxt="n",xlab="")
lines(x,ye)
legend("bottomright",c("modulated","envelope"),col=c("grey","black"),lty=c(2,1))
par(mar=c(3.5,3.5,2,0))
wf <- waterfall(y,x,nf = 3)
# rasterImage2(x = wf$x, y = wf$fx, z = wf$A
# ,ylim = c(0,60))
plot(wf,ylim=c(0,40),main="Waterfall")
#### uncertainty principle ####
#
# take a look at the side effects
# at [0,30] and [1,0]
#
# With a large steepness e.g. n = 50 you will gain
# artefacts.
#
# if frequency is not stationary
# PSD becomes > 1 depending on the type of band filter.
#
###############################
x <- seq(0,1, length.out=1500)
y <- sin(100*x*x)
FT <- spec.fft(x = x, y = y)
wf <- waterfall(y,x)
par(mfrow=c(2,1),mar=c(1,3.5,3,3),mgp=c(2.5,1,0))
# plot results
plot(x,y,type="l",lwd=1,col="darkgrey",lty=2,ylab="y",main="Original Data",xaxt="n",xlab="")
par(mar=c(3.5,3.5,2,0))
plot(wf
,ylim=c(0,40),main="Waterfall"
)
abline(h = 25, lty = 3, lwd = 3, col = "grey")
range(wf$PSD,na.rm = TRUE)
range(wf$A)
###### effect of missing values #####
#
# 10% random missing values cause a
# distortion and a miss scaling of
# the PSD value, which becomes >1 now.
# This depends on the type of band pass
# filter selected.
#
#####################################
x <- seq(0,5, length.out=500)
y <- sin(2*pi * 15 * x + 2*1*cos(2*pi*0.5*x))
# delete 10% of the data
y[sample(length(y),size = 50)] <- NA
wf <- waterfall(y,x,type = "b")
par(mfrow=c(2,1),mar=c(1,3.5,3,3),mgp=c(2.5,1,0))
# plot results
plot(x,y,type="l",lwd=1,col="darkgrey",lty=2,ylab="y",main="Original Data",xaxt="n",xlab="")
par(mar=c(3.5,3.5,2,0))
plot(wf
,ylim=c(10,20),main="Waterfall"
)
abline(h = 25, lty = 3, lwd = 3, col = "grey")
# check the PSD range
range(wf$PSD)
range(wf$A)