class GSL::SpectralAnalysis::Lomb

A Lomb periodogram is a method of spectrally analysing a set of data which are not evenly spaced, and thus cannot be Fourier transformed. The Lomb periodogram is something akin to a probability distribution function for a given set of frequencies.

Attributes

frequencies[RW]
periodogram[RW]

Public Class Methods

new(times, data) click to toggle source

Create a new Lomb object. times and data should be GSL::Vectors.

# File lib/gsl_extras.rb, line 986
def initialize(times, data)
        @times = times; @data = data
        raise "Times #{times.size} and data #{data.size} do not have the same size" unless @times.size == @data.size
        @n = data.size
        @dmean = data.mean
        @dvar = data.variance
end

Public Instance Methods

calculate_periodogram(frequency_factor=2.0, oversampling=4.0, frequency_indexes=nil) click to toggle source

Calculate the Lomb periodogram. Without studying the Lomb analysis, it's best to leave the defaults alone. frequency_factor is how far above the estimated Nyquist frequency (calculated by dividing the net time interval by the number of data points) the spectrum should be calculated.

# File lib/gsl_extras.rb, line 998
        def calculate_periodogram(frequency_factor=2.0, oversampling=4.0, frequency_indexes=nil)
                @frequency_factor = frequency_factor
                @oversampling = oversampling
                @nout = (@n * 0.5 * frequency_factor * oversampling).to_i # (nout or @n * 0.5 * frequency_factor * 4.0).to_i
                t_window = @times.max - @times.min
                delta_f = 1.0 / t_window / oversampling # / 2.0 / Math::PI #(@nout / @n / 0.5 / frequency_factor)
#                       data_min, data_max = @data.minmax
                
                @frequencies = GSL::Vector.linspace(delta_f, delta_f*@nout, @nout)
#               p @nout, delta_f, @frequencies
                if frequency_indexes 
                        @frequencies = GSL::Vector.alloc(frequency_indexes.map{|i| @frequencies[i]})
                end
                @periodogram = @frequencies.collect do |freq|
                        p_n(freq)
                end
#               @frequencies = @frequencies / Math::PI / 2.0
                [@frequencies, @periodogram]
        end
confidence(pn, frequency_factor = @frequency_factor) click to toggle source

Equal to 1.0 - the probability that the value of pn could have been generated by gaussian noise

# File lib/gsl_extras.rb, line 1039
def confidence(pn, frequency_factor = @frequency_factor)
        (1.0 - Math.exp(-pn)) ** (@n  * frequency_factor)
end
graphkit() click to toggle source
# File lib/gsl_extras.rb, line 1056
def graphkit
        CodeRunner::GraphKit.autocreate(x: {title: "Frequency", data: @frequencies}, y: {title: "P_N", data: @periodogram})
end
p_n(freq) click to toggle source

Proportional to the probability that the given frequency was present in the data. Roughly akin to p(k) for a Fourier transform.

# File lib/gsl_extras.rb, line 1020
def p_n(freq)
        omega = freq * Math::PI * 2.0
        twoomt = @times * 2.0 * omega
        tau = Math.atan(
                twoomt.sin.sum / twoomt.cos.sum
                )/ 2.0 / omega
        omttau = ((@times - tau) * omega)
        c = omttau.cos
        s = omttau.sin
        ddmean = @data - @dmean
        pn = 1 / 2.0 / @dvar * (
                (ddmean * c).sum ** 2.0 / c.square.sum +
                (ddmean * s).sum ** 2.0 / s.square.sum
        )
        pn
end
p_n_from_confidence(confidence, frequency_factor = @frequency_factor) click to toggle source

Find a

# File lib/gsl_extras.rb, line 1051
def p_n_from_confidence(confidence, frequency_factor = @frequency_factor)
        - Math.log(1.0 - confidence ** (1.0 / @n / frequency_factor))
end
pnull(pn, frequency_factor = @frequency_factor) click to toggle source

The probability that the value of pn could have been generated by gaussian noise.

# File lib/gsl_extras.rb, line 1045
def pnull(pn, frequency_factor = @frequency_factor)
        1.0 - confidence(pn, frequency_factor)
end