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
Public Class Methods
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 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
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
# File lib/gsl_extras.rb, line 1056 def graphkit CodeRunner::GraphKit.autocreate(x: {title: "Frequency", data: @frequencies}, y: {title: "P_N", data: @periodogram}) end
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
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
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