Weighted mean spectrum

Calculate the mean spectrum from a list of individual spectra including weighting of the individual spectra and, if needed, interpolation to the highest resolution frequency grid. The default weighting produces the simple arithmetic mean. Interpolation is only performed when the frequency axes of the individual spectra have different lengths and/or differ in frequency discretization.

Usage

MeanSpectrum(specList, iRemoveLowest = 1, weights = rep(1, length(specList)))

Arguments

specList

list of spectra, i.e. objects of class "spec" (see SpecMTM for details), where each spectrum has to be a list of the vectors freq, spec and dof of a common length.

iRemoveLowest

integer; number of lowest frequencies to remove from each individual spectral estimate (e.g. to remove detrending bias) prior to the interpolation and averaging.

weights

numeric vector of weights; its length must match the number of elements in specList. The default setting of identical weights yields the simple arithmetic average of the input spectra; change this parameter accordingly to apply non-arithmetic averaging (see example).

Value

object of class "spec" with the weighted mean spectrum, amended by the element nRecord which gives the number of records contributing to each mean spectral estimate.

Author

Thomas Laepple and Thomas Münch

Examples

# Simple arithmetic average
f1 <- 1 : 5
f2 <- f1
s1 <- rep(1, length(f1))
s2 <- rep(3, length(f2))
dof1 <- rep(1, length(f1))
dof2 <- rep(1, length(f2))

spectra <- list(list(freq = f1, spec = s1, dof = dof1),
                list(freq = f2, spec = s2, dof = dof2))

MeanSpectrum(spectra, iRemoveLowest = 0)
#> $freq
#> [1] 1 2 3 4 5
#> 
#> $spec
#> [1] 2 2 2 2 2
#> 
#> $dof
#> [1] 2 2 2 2 2
#> 
#> $nRecord
#> [1] 2 2 2 2 2
#> 
#> $lim.1
#> [1] 7.377759 7.377759 7.377759 7.377759 7.377759
#> 
#> $lim.2
#> [1] 0.05063562 0.05063562 0.05063562 0.05063562 0.05063562
#> 
#> $pval
#> [1] 0.05
#> 
#> attr(,"class")
#> [1] "spec"

# Weighted mean with interpolation
f1 <- 1 : 5
f2 <- 3 : 8
s1 <- rep(1, length(f1))
s2 <- rep(3, length(f2))
dof1 <- rep(1, length(f1))
dof2 <- rep(1, length(f2))

spectra <- list(list(freq = f1, spec = s1, dof = dof1),
                list(freq = f2, spec = s2, dof = dof2))

MeanSpectrum(spectra, iRemoveLowest = 0, weights = c(1, 2))
#> $freq
#> [1] 1 2 3 4 5 6 7 8
#> 
#> $spec
#> [1] 1.000000 1.000000 2.333333 2.333333 2.333333 3.000000 3.000000 3.000000
#> 
#> $dof
#> [1] 1 1 2 2 2 1 1 1
#> 
#> $nRecord
#> [1] 1 1 2 2 2 1 1 1
#> 
#> $lim.1
#> [1]  5.023886  5.023886  8.607385  8.607385  8.607385 15.071659 15.071659
#> [8] 15.071659
#> 
#> $lim.2
#> [1] 0.0009820691 0.0009820691 0.0590748853 0.0590748853 0.0590748853
#> [6] 0.0029462074 0.0029462074 0.0029462074
#> 
#> $pval
#> [1] 0.05
#> 
#> attr(,"class")
#> [1] "spec"
# with some detrending bias removal
MeanSpectrum(spectra, iRemoveLowest = 1, weights = c(1, 2))
#> $freq
#> [1] 2 3 4 5 6 7 8
#> 
#> $spec
#> [1] 1.000000 1.000000 2.333333 2.333333 3.000000 3.000000 3.000000
#> 
#> $dof
#> [1] 1 1 2 2 1 1 1
#> 
#> $nRecord
#> [1] 1 1 2 2 1 1 1
#> 
#> $lim.1
#> [1]  5.023886  5.023886  8.607385  8.607385 15.071659 15.071659 15.071659
#> 
#> $lim.2
#> [1] 0.0009820691 0.0009820691 0.0590748853 0.0590748853 0.0029462074
#> [6] 0.0029462074 0.0029462074
#> 
#> $pval
#> [1] 0.05
#> 
#> attr(,"class")
#> [1] "spec"