
"""
    gcve_pdf(epochlength::Real,fs::Int;filter=nothing,itr=1000000,pruned=0.5, σ = 1)

Calculate the distribution of the coefficient of variation of the envelope (CVE) of discrete-time zero-mean white Gaussian noise using random samples from N(0,1).

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `filter`: A filter function
* `itr`: The number of epochs used to compute the distribution
* `pruned`: The number of seconds cut from both tails of the envelope before calculating the CVE
* `σ`: Standard deviation of the Normal distribution

# Notes
* Notice that the pdf for the CVE of white Gaussian noise can also be bootstrapped using Fourier transform phase randomization.

# Examples
```julia
julia> gcve_pdf(10,1000,itr=10)
10-element Array{Float64,1}:
 0.5246024892582911
 0.52295322056966
 0.5256115270351591
 0.5263961505330913
 ⋮
 0.5237027693791765
 0.5181673343579238
 0.5260079659619668
```
"""
function gcve_pdf(epochlength::Real,fs::Int; filter = nothing,itr = 1000000,pruned = 0.5,σ = 1)
    epochlength>2*pruned ? nothing : error("Epoch length too short")
	isinteger(epochlength*fs) ? nothing : error("Number of points must be integer")
	isinteger(pruned*fs) ? nothing : error("Prunning points must be integer")
    println("The effective epoch length after prunning is $(epochlength-2*pruned) seconds")

	points = Int(epochlength*fs)
	d = Normal(0,σ)
    gcve_expdistro = zeros(Float64,itr)

    if filter != nothing
    	println("\n","Calculating CVE pdf for filtered white Gaussian noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            gaussian_sample = rand(d, points)
            gaussian_sample = filter(gaussian_sample,fs)
            gaussian_sample = abs.(DSP.hilbert(gaussian_sample))
            gaussian_sample = gaussian_sample[Int(fs*pruned)+1:end-Int(fs*pruned)]
            gcve_expdistro[i] = std(gaussian_sample)/mean(gaussian_sample)
        end
    else
	println("\n","Calculating CVE pdf for non-filtered white Gaussian noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            gaussian_sample = rand(d, epochlength*fs)
            gaussian_sample = abs.(DSP.hilbert(gaussian_sample))
            gaussian_sample = gaussian_sample[Int(fs*pruned)+1:end-Int(fs*pruned)]
            gcve_expdistro[i] = std(gaussian_sample)/mean(gaussian_sample)
        end
    end
    return gcve_expdistro
end


"""
    rcv_pdf(points::Int; itr=1000000, σ = 1)

Calculate the distribution of the coefficient of variation of Rayleigh noise using random samples from Rayleigh(σ).

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `itr`: The number of epochs used to compute the distribution
* `σ`: Scale parameter of the Rayleigh distribution

# Examples
```julia
julia> rcv_pdf(10000,itr=10,σ = 2)
10-element Array{Float64,1}:
 0.5208939217095829
 0.523914916275298
 0.5242368643462019
 0.519126955367746
 ⋮
 0.5229001421212949
 0.5239581722102657
 0.5275819378236367
```
"""
function rcv_pdf(points::Int;itr=1000000, σ = 1)

	d = Rayleigh(σ)
    rcv_expdistro = zeros(Float64,itr)

    println("\n","Calculating CV pdf for Rayleigh noise...")
    @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            rayleigh_sample = rand(d, points)
            rcv_expdistro[i] = std(rayleigh_sample)/mean(rayleigh_sample)
    end
    return rcv_expdistro
end


"""
    ucve_pdf(epochlength::Int,fs::Int;filter=nothing,itr=1000000,prunning=1)

Calculate the distribution of the coefficient of variation of the envelope (CVE) of discrete-time zero-mean uniform noise.

# Arguments
* `epochlength`: The length of the epoch in seconds
* `fs`: The sampling frequency
* `filter`: A filter function
* `itr`: The number of epochs used to compute the distribution
* `prunning`: The number of seconds cut from both tails of the envelope before calculating the CVE

# Notes
* Notice that the pdf for the CVE of white Gaussian noise can also be bootstrapped using Fourier transform phase randomization.

# Examples
```julia
julia> ucve_pdf(10,1000,itr=10)
10-element Array{Float64,1}:
 0.42808995848724163
 0.4353783310060235
 0.4341374928673676
 0.43346989277220405
 ⋮
 0.4317732242362664
 0.42726595981167065
 0.4327565016764864
```
"""
function ucve_pdf(epochlength::Int,fs::Int;filter=nothing,itr=1000000,prunning=1)
    epochlength>2*prunning ? nothing : error(" Epoch length too short")
    println("The effective epoch length after prunning is $(epochlength-prunning) seconds")
    d = Uniform(-1,1)
    ucve_expdistro = zeros(Float64,itr)


    if filter != nothing
    	println("\n","Calculating CVE pdf for filtered uniform noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
            uniform_sample = rand(d, epochlength*fs)
            uniform_sample = filter(uniform_sample,fs)
            uniform_sample = abs.(DSP.hilbert(uniform_sample))
            uniform_sample = uniform_sample[(fs*prunning)+1:end-(fs*prunning)]
            ucve_expdistro[i] = std(uniform_sample)/mean(uniform_sample)
        end
    else
	println("\n","Calculating CVE pdf for non-filtered uniform noise...")
        @inbounds for i in 1:itr
            println("Iteration $(i)/$(itr)")
			uniform_sample = rand(d, epochlength*fs)
            uniform_sample = abs.(DSP.hilbert(uniform_sample))
            uniform_sample = uniform_sample[(fs*prunning)+1:end-(fs*prunning)]
            ucve_expdistro[i] = std(uniform_sample)/mean(uniform_sample)
        end
    end
    return ucve_expdistro
end
