Nothing
# tfestimate.R
# Copyright (C) 2020 Geert van Boxtel <gjmvanboxtel@gmail.com>
# Original Octave code:
# Copyright (C) 2006 Peter V. Lanspeary <pvl@mecheng.adelaide.edu.au>
#
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU General Public License
# as published by the Free Software Foundation; either version 3
# of the License, or (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# Version history
# 20201107 GvB setup for gsignal v0.1.0
#------------------------------------------------------------------------------
#' Transfer Function Estimate
#'
#' Finds a transfer function estimate for signals.
#'
#' \code{tfestimate} uses Welch's averaged periodogram method.
#'
#' @param x input data, specified as a numeric vector or matrix. In case of a
#' vector it represents a single signal; in case of a matrix each column is a
#' signal.
#' @param window If \code{window} is a vector, each segment has the same length
#' as \code{window} and is multiplied by \code{window} before (optional)
#' zero-padding and calculation of its periodogram. If \code{window} is a
#' scalar, each segment has a length of \code{window} and a Hamming window is
#' used. Default: \code{nextpow2(sqrt(length(x)))} (the square root of the
#' length of \code{x} rounded up to the next power of two). The window length
#' must be larger than 3.
#' @param overlap segment overlap, specified as a numeric value expressed as a
#' multiple of window or segment length. 0 <= overlap < 1. Default: 0.5.
#' @param nfft Length of FFT, specified as an integer scalar. The default is the
#' length of the \code{window} vector or has the same value as the scalar
#' \code{window} argument. If \code{nfft} is larger than the segment length,
#' (seg_len), the data segment is padded \code{nfft - seg_len} zeros. The
#' default is no padding. Nfft values smaller than the length of the data
#' segment (or window) are ignored. Note that the use of padding to increase
#' the frequency resolution of the spectral estimate is controversial.
#' @param fs sampling frequency (Hertz), specified as a positive scalar.
#' Default: 1.
#' @param detrend character string specifying detrending option; one of:
#' \describe{
#' \item{\code{"long-mean"}}{remove the mean from the data before
#' splitting into segments (default)}
#' \item{\code{"short-mean"}}{remove the mean value of each segment}
#' \item{\code{"long-linear"}}{remove linear trend from the data before
#' splitting into segments}
#' \item{\code{"short-linear"}}{remove linear trend from each segment}
#' \item{\code{"none"}}{no detrending}
#' }
#'
#' @return A list containing the following elements:
#' \describe{
#' \item{\code{freq}}{vector of frequencies at which the spectral variables
#' are estimated. If \code{x} is numeric, power from negative frequencies is
#' added to the positive side of the spectrum, but not at zero or Nyquist
#' (fs/2) frequencies. This keeps power equal in time and spectral domains.
#' If \code{x} is complex, then the whole frequency range is returned.}
#' \item{\code{trans}}{NULL for univariate series. For multivariate series,
#' a matrix containing the transfer function estimates between different
#' series. Column \eqn{i + (j - 1) * (j - 2)/2 } of \code{coh} contains the
#' cross-spectral estimates between columns \eqn{i} and \eqn{j} of \eqn{x},
#' where \eqn{i < j}.}
#' }
#'
#' @examples
#' fs <- 1000
#' f <- 250
#' t <- seq(0, 1 - 1/fs, 1/fs)
#' s1 <- sin(2 * pi * f * t) + runif(length(t))
#' s2 <- sin(2 * pi * f * t - pi / 3) + runif(length(t))
#' rv <- tfestimate(cbind(s1, s2), fs = fs)
#' plot(rv$freq, 10*log10(abs(rv$trans)), type="l", xlab = "Frequency",
#' ylab = "Tranfer Function Estimate (dB)", main = colnames((rv$trans)))
#'
#' h <- fir1(30, 0.2, window = rectwin(31))
#' x <- rnorm(16384)
#' y <- filter(h, x)
#' tfe <- tfestimate(cbind(x, y), 1024, fs = 500)
#' plot(tfe$freq, 10*log10(abs(tfe$trans)), type="l", xlab = "Frequency",
#' ylab = "Tranfer Function Estimate (dB)", main = colnames((tfe$trans)))
#'
#' @note The function \code{tfestimate} (and its deprecated alias \code{tfe})
#' is a wrapper for the function \code{pwelch}, which is more complete and
#' more flexible.
#'
#' @author Peter V. Lanspeary, \email{pvl@@mecheng.adelaide.edu.au}.\cr
#' Conversion to R by Geert van Boxtel, \email{G.J.M.vanBoxtel@@gmail.com}.
#'
#' @seealso \code{\link{pwelch}}
#'
#' @rdname tfestimate
#' @export
tfestimate <- function(x, window = nextpow2(sqrt(NROW(x))), overlap = 0.5,
nfft = ifelse(isScalar(window), window, length(window)),
fs = 1,
detrend = c("long-mean", "short-mean",
"long-linear", "short-linear", "none")) {
pw <- pwelch(x, window, overlap, nfft, fs, detrend)
rv <- list(freq = pw$freq, trans = pw$trans)
rv
}
#' @rdname tfestimate
#' @export
tfe <- tfestimate
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.