filterTSeriesSSA(series, borders.wl, M = rep(floor(length(series)/3), times = n.steps), n.comp = rep(40, times = n.steps), harmonics = rep(0, times = n.steps), tolerance.harmonics = 0.05, var.thresh.ratio = 0.005, grouping = c("grouping.auto", "groupSSANearestNeighbour")[1], groupingMethod = "wcor", repeat.extr = rep(1, times = length(borders.wl)), recstr.type = "subtraction", pad.series = c(0, 0), SSA.methods = c("nutrlan", "propack", "eigen", "svd"), center.series = TRUE, call.freq = quote(calcFrequency(series.t)), n.steps = switch(class(borders.wl), list = length(borders.wl), dim(borders.wl)[2]), plot.spectra = FALSE, second.axis = TRUE, open.plot = TRUE, print.stat = TRUE, xlim = c(), debugging = FALSE, ...)Definition of the period borders (borders.wl):
borders.wl contains the borders of the different periodicity bands to
extract. Units are sampling frequency of the series. borders.wl has to be a
list with one element per desired decomposition step (see also 'stepwise
extraction'). In the case of a single band to be extracted, this vector has
to consist of two values (e.g. c(
Window length (M): For optimal detection of spectral components the window length or embedding dimension M should be an integer multiple of the period of the oscillation that is to be extracted. See also ?ssa (here the same parameter is called L)
Padding (pad.series): For padding, the series should start and end exactly at the start and end of a major oscillation (e.g. a yearly cycle). Additionally the length to use for padding should be a integer multiple of this period length. The padding is solved internally by adding the indicated part of the series at the start and at the end of the series. This padded series is only used internally and only the part of the series with original data is returned in the results.
High frequency part (recstr.type): Two different ways to compute the high frequency part (e.g. the part with a frequency between 0 and the lowest border.wl value) are implemented. The standard way (if recstr.type == 'subtraction') is to sum all eigentriples with lower frequencies and subtract this sum from the original series to compute this as the high frequency residual. Otherwise (if (recstr.type != 'subtraction')) only the eigentriples with such high frequencies that were actually extracted are used to build this spectral band.
Stepwise extraction: In general the whole algorithm can be run stepwise. This means that first a certain spectral band is computed (with all possible harmonies etc. ...) and subtracted from the original series. Subsequently this process can be repeated with the residual as often as wanted. This allows for the adaptation of, for example, M, harmonics, ... to the particular oscillation to be extracted. Additionally it often this often leads to a clearer signal separation. To implement several steps, each of M, n.comp and harmonics needs to be a vector of the length of the amount of steps. For each step the corresponding element of this vector is used. borders.wl needs to contain one list entry per step which needs to be a vector containing all borders of the bands to extract in this step (see 'Definition of the period borders').
Repeated extraction (repeat.extr): Especially for the trend it may be advisable to repeat the filtering step for this particular band. In this case the result of the first filtering (e.g. the sum of all eigentriples within this band) is filtered again n times with the same period borders. Finally only the final filtered components are summed and treated as the actual spectral band. In many cases this is helpful to exclude high frequency parts from the trend or low frequency components. It is only possible to repeat the extraction for steps where single bands are extracted.
Visualize results (plot.spectra): In the case that diagnostic plots should be plotted (plot.spectra == TRUE) one (or more) pseudospectra are plotted. Each point in these represents one group of eigentriples determined. For each step a separate plot is produced. Colored regions represent the specified spectral bands. grey lines in the background represent a Fourier Spectrum of the original series.
ssa,calcFrequency#create series consisting of two oscillations and noise
series.ex <- sin(2 * pi * 1:1000 / 100) + 0.7 * sin(2 * pi * 1:1000 / 10) +
rnorm(n = 1000, sd = 0.4)
#prepare graphics
layout(matrix(c(1, 2, 3, 4, 5, 6, 7, 8), ncol = 2))
par(tcl = 0.2, mgp = c(2, 0, 0), mar = c(0, 4, 0, 0), oma = c(2, 0, 2, 0),
ps = 10, cex = 1)
plot.new()
#perform decomposition
data.decomposed <- filterTSeriesSSA(series = series.ex,
borders.wl = list(a = c(8, 12), b = c(80, 120)
, c = c(0, 10, 100, Inf)),
M = c(30, 200, 100),
n.comp = c(10, 20, 20),
harmonics = c(1, 0, 0),
plot.spectra = TRUE, open.plot = FALSE)
#plot series and spectral parts
plot(series.ex)
plot(data.decomposed$dec.series[1, ], ylab = '')
plot(data.decomposed$dec.series[2, ], ylab = '')
plot(colSums(data.decomposed$dec.series[-c(1:2), ]), ylab = '')
mtext(side = 2, outer = TRUE, at = -(1 / 8) + ((4:1) * (1 / 4)),
c('orig.series', '1.step', '2.step', '3.step'), las = 3, cex = 1.5, line = -1)
mtext(side = 3, outer = TRUE, at = -(1 / 4) + ((1:2) * (1 / 2)),
c('pseudospectra', 'identified components'), las = 1, cex = 1.5, line = 1)
Run the code above in your browser using DataLab