Class: GSL::SpectralAnalysis::Lomb

Inherits:
Object
  • Object
show all
Defined in:
lib/gsl_extras.rb

Overview

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.

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize(times, data) ⇒ Lomb

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



986
987
988
989
990
991
992
 
# 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

Instance Attribute Details

#frequencies ⇒ Object

Returns the value of attribute frequencies.



994
995
996
 
# File 'lib/gsl_extras.rb', line 994

def frequencies
  @frequencies
end

#periodogram ⇒ Object

Returns the value of attribute periodogram.



994
995
996
 
# File 'lib/gsl_extras.rb', line 994

def periodogram
  @periodogram
end

Instance Method Details

#calculate_periodogram(frequency_factor = 2.0, oversampling = 4.0, frequency_indexes = nil) ⇒ Object

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.



998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
 
# 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) ⇒ Object

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



1039
1040
1041
 
# File 'lib/gsl_extras.rb', line 1039

def confidence(pn, frequency_factor = @frequency_factor)
	(1.0 - Math.exp(-pn)) ** (@n  * frequency_factor)
end

#graphkit ⇒ Object 



1056
1057
1058
 
# 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) ⇒ Object

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



1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
 
# 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) ⇒ Object

Find a



1051
1052
1053
 
# 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) ⇒ Object

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



1045
1046
1047
 
# File 'lib/gsl_extras.rb', line 1045

def pnull(pn, frequency_factor = @frequency_factor)
	1.0 - confidence(pn, frequency_factor)
end