R/rexp_MultSines.R

In arnejohannesholmin/echoIBM: Simulation of multibeam sonar data

#*********************************************
#*********************************************
#' Returns an array of correlated vectors of exponentially distributed variables, where each vector has autocorrelation resulting from superimposed sine waves of equal length, frequency and amplitude, but with uniformly distributed starting point. The values are generated by dividing the time into intervals, and extracting the central peak of the superimposed wave for each sampling interval. Correlation between vectors is acomplished by also including sine waves from other vectors in the superimposed wave, to a degree given by 'olpn'. For low numbers of sine waves the square of the probability distribution described by Barakat 1974 emerges, but strong deviations from the exponential distribution requires very low values of 'L' and 'w' and few elements in 'olpn', as the Rayleigh distribution kicks in as low as around 4 individual waves. WARNING: The expectation of the exponential variables is dependent on the parameters used, and for the default setting it is equal to 7.350175.
#'
#' @param J  is the length of the beams.
#' @param I  is the number of beams.
#' @param L  is the number of fish in each voxel.
#' @param N  is the number of sample points in each voxel (N=4*P, the sampling points are expected to end up close to peaks along the sine waves).
#' @param P  is the is the frequency of the sine waves, meaning the number of periods per voxel.
#' @param w  is the length of the sine waves in units of the time intervals constituting the voxels.
#' @param olpn  is a vector, matrix or array of correlations, representing constant correlation for all voxels, variable correlation between beams, and variable correlation for all voxels.
#' @param seed  is the seed of the function.
#' @param shape  is is used to transform to a Weibull variable of the desired shape and scale=1.
#'
#' @return
#'
#' @examples
#' \dontrun{}
#'
#' @importFrom TSD dim_all zeros
#' @importFrom stats runif
#'
#' @export
#' @rdname rexp_MultSines
#'
rexp_MultSines<-function(J=1324, I=25, L=3, N=40, P=10, w=3, olpn=c(0.5,1,0.5), shape=1, mean=1, seed=NULL){
	
	############### LOG: ###############
	# Start: 2012-01-24 - Clean version.
	# Update: 2012-02-17 - Added the option 'buffer'.
	# Update: 2012-02-26 - The c++ function rexp_MultSines.cpp was changed to identify the central peak instead of averaging between a number of peaks near the center of voxels. This saved time, and the output out does needs not be devided by the number of peaks. Also deleted the parameter 'buffer'.
	# Update: 2013-09-23 - Implemented expected value equal to 1 inside the function.
	# Last: 2013-12-04 - Added 'shape'.
		

	##################################################
	##################################################
	##### Preparation #####
	# Define the output array here, for use in the call of the C++ function:
	out <- zeros(J*I)
	# Check whether 'olpn' has the correct length, which is (2 * ci + 1) * I:
	if(length(dim(olpn))==0){
		olpn <- matrix(olpn, nrow=length(olpn), ncol=I)
	}
	ci <- (NROW(olpn)-1)/2
	# Define the minimum distance between maxima in the superimposed wave, used for extracting only the proper peaks of the wave, and not all the minor peaks (if any):
	wM <- N/P/2
	# Set the seed if missing:
	if(length(seed)==0){
		seed <- runif(1, 0, 1e9)
	}
	# Check for consistancy of 'olpn' and the dimension specifyers 'J' and 'I':
	if(length(dim(olpn))>2 && any(J!=dim(olpn)[2], I!=dim(olpn)[3])){
		stop(paste("The dimension specifyers 'J' (",J,") and 'I' (",I,") do not match the last two dimensions of 'olpn' (",paste(dim(olpn)[-1],collapse=", "),")",sep=""))
	}	
	
	##### Execution #####
	if(length(dim(olpn))==3){
		U <- .C("rexp_MultSines_VaribleCorrelation", as.integer(J), as.integer(I), as.integer(L), as.integer(N), as.double(P), as.double(w), as.integer(ci), as.double(olpn), as.double(out), as.integer(wM), as.integer(seed), PACKAGE="echoIBM")
	}
	else{
		U <- .C("rexp_MultSines", as.integer(J), as.integer(I), as.integer(L), as.integer(N), as.double(P), as.double(w), as.integer(ci), as.double(olpn), as.double(out), as.integer(wM), as.integer(seed), PACKAGE="echoIBM")
	}
	# Insert the simulated noise data:
	out <- U[[9]]
	dim(out) <- c(J,I)
	
	# Get the scale used to assure expected value equal to 1:
	scale_exp1 <- echoIBM.w2scale(w)
		
	# Constant overlap for all beams and all voxels. Dimension of 'olpn': [number of adjacet beams correlated]:
	if(length(dim_all(olpn))<=1){
		scale_exp1 <- scale_exp1 * sum(olpn^2)
	}
	# Constant overlap for all beams and all voxels within fans, but between fans. Dimension of 'olpn': [number of adjacet beams correlated, number of beams]:
	else if(length(dim_all(olpn))==2){
		scale_exp1 <- scale_exp1 * colSums(olpn^2)
	}
	# Variable overlap for all beams and constant overlap for all voxels within beams. Dimension of 'olpn': [number of adjacet beams correlated, number of voxels along the beams, number of beams]:
	else if(length(dim_all(olpn))==3){
		scale_exp1 <- scale_exp1 * colSums(olpn^2)
	}
	# Variable overlap for all voxels, only used with the periodic noise. Dimension of 'olpn': [number of adjacet beams correlated, number of voxels along the beams, number of beams in each fan of constant frequency, number of fans of constant frequency]:
	else if(length(dim_all(olpn))==4){
		scale_exp1 <- scale_exp1 * apply(olpn[,seq_len(J),]^2, 2:3, sum)
	}
	else{
		stop("Only overlap matrices up to 4 dimensions are supported")
	}
	
	
	##### Output #####
	# Divide by the expected value of the exponential variables to obtain mean=1, and return:
	out <- out/scale_exp1
	
	# Include the shape when transforming to Weibull with scale parameter = 1:
	if(!identical(shape,1)){
		out <- out^(1/shape) / gamma(1 + 1/shape)
	}
	else{
		out <- out
	}
	# Apply the mean by multiplying by 'mean' both for exponential and Weibull:
	if(!identical(mean, 1)){
		if(length(mean)<2){
			mean <- matrix(mean, byrow=TRUE, ncol=ncol(out), nrow=nrow(out))
		}
		out <- out * mean
	}
	return(out)
	##################################################
	##################################################
	}