Vector
obj = ysrcdest.interpolate(xdest, xsrc)obj = ydest.interpolate(xdest, xsrc, ysrc)//... objref xs, ys, xd, yd xs = new Vector(10) xs.indgen() ys = xs.c.mul(xs) ys.line(g, xs, 1, 0) // black reference line xd = new Vector() xd.indgen(-.5, 10.5, .1) yd = ys.c.interpolate(xd, xs) yd.line(g, xd, 3, 0) // blue more points than reference xd.indgen(-.5, 13, 3) yd = ys.c.interpolate(xd, xs) yd.line(g, xd, 2, 0) // red fewer points than reference
Vector
obj = vdest.deriv(vsrc)obj = vdest.deriv(vsrc, dx)obj = vdest.deriv(vsrc, dx, method)obj = vsrcdest.deriv()obj = vsrcdest.deriv(dx)obj = vsrcdest.deriv(dx, method)vec
is placed in vdest.
The variable dx gives the increment of the independent variable
between successive elements of vec.
method = 1 = Euler derivative:
vec1[i] = (vec[i+1] - vec[i])/dx vec is lost since i cannot equal -1. Therefore, since the
integral function performs an Euler
integration, the integral of vec1 will reproduce vec minus the first
element.
method = 2 = Central difference derivative:
vec1[i] = ((vec[i+1]-vec[i-1])/2)/dx vec1. The central difference method maintains the
same number of elements in vec1
as were in vec and is a more accurate method than the Euler method.
A vector differentiated by this method cannot, however, be integrated
to reproduce the original vec.
objref vec, vec1
vec = new Vector()
vec1 = new Vector()
vec.indgen(0, 5, 1)
func sq(){
return $1*$1
}
vec.apply("sq")
vec1.deriv(vec, 0.1)
vec1 with elements:
10 20 40 60 80 90Since
dx=0.1, and there are eleven elements including 0,
the entire function exists between the values of 0 and 1, and the derivative
values are large compared to the function values. With dx=1,the vector
vec1 would consist of the following elements:
1 2 4 6 8 9
The Euler method vs. the Central difference method:
Beginning with the vector vec:
0 1 4 9 16 25
vec1.deriv(vec, 1, 1) (Euler) would go about
producing vec1 by the following method:
1-0 = 1 4-1 = 3 9-4 = 5 16-9 = 7 25-16 = 9whereas
vec1.deriv(vec, 1, 2) (Central difference) would go about
producing vec1 as such:
1-0 = 1 (4-0)/2 = 2 (9-1)/2 = 4 (16-4)/2 = 6 (25-9)/2 = 8 25-16 = 9
Vector
obj = vdest.integral(vsrc)obj = vdest.integral(vsrc, dx)obj = vsrcdest.integral()obj = vsrcdest.integral(dx)dx sets the size of the discretization.
vdest[i+1] = vdest[i] + vsrc[i+1] and the first element of vdest is always
equal to the first element of vsrc.
objref vec, vec1
vec = new Vector()
vec1 = new Vector()
vec.indgen(0, 5, 1) //vec will have 6 values from 0 to 5, with increment=1
vec.apply("sq") //sq() squares an element
//and is defined in the example for .deriv
vec1.integral(vec, 1) //Euler integral of vec elements approximating
//an x-squared function, dx = 0.1
vec1.printf()
vec1 to the screen:
0 1 5 14 30 55In order to make the integral values more accurate, it is necessary to increase the size of the vector and to decrease the size of dx.
objref vec2
vec2 = new Vector(6)
vec.indgen(0, 5.1, 0.1) //vec will have 51 values from 0 to 5, with increment=0.1
vec.apply("sq") //sq() squares an element
//and is defined in the example for .deriv
vec1.integral(vec, 0.1) //Euler integral of vec elements approximating
//an x-squared function, dx = 0.1
for i=0,5{vec2.x[i] = vec1.x[i*10]} //put the value of every 10th index in vec2
vec2.printf()
vec2 (which are the elements of
vec1 corresponding to the integers 0-5) to the screen:
0 0.385 2.87 9.455 22.14 42.925The integration naturally becomes more accurate as dx is reduced and the size of the vector is increased. If the vector is taken to 501 elements from 0-5 and dx is made to equal 0.01, the integrals of the integers 0-5 yield the following (compared to their continuous values on their right).
0.00000 -- 0.00000 0.33835 -- 0.33333 2.6867 -- 2.6666 9.04505 -- 9.00000 21.4134 -- 21.3333 41.7917 -- 41.6666
Vector
median = vsrc.median()vec.
Vector
obj = vdest.medfltr(vsrc)obj = vdest.medfltr(vsrc, points)obj = vsrcdest.medfltr()obj = vsrcdest.medfltr( points)
Vector
obj = vsrcdest.sort()vec1 in place, putting them in numerical order.
Vector
vdest = vsrc.sortindex()vdest = vsrc.sortindex(vdest)objref a, r, si r = new Random() r.uniform(0,100) a = new Vector(10) a.setrand(r) a.printf si = a.sortindex si.printf a.index(si).printf
Vector
obj = vsrcdest.reverse()vec in place.
Vector
obj = vsrcdest.rotate(value)obj = vsrcdest.rotate(value, 0)orders the elements ofvec.indgen(1, 10, 1) vec.rotate(3)
vec as follows:
8 9 10 1 2 3 4 5 6 7whereas,
orders the elements ofvec.indgen(1, 10, 1) vec.rotate(-3)
vec as follows:
4 5 6 7 8 9 10 1 2 3execute following example
objref vec vec = new Vector() vec.indgen(1,5,1) vec.printf vec.c.rotate(2).printf vec.c.rotate(2, 0).printf vec.c.rotate(-2).printf vec.c.rotate(-2, 0).printf
Vector
obj = vdest.rebin(vsrc,factor)obj = vsrcdest.rebin(factor)vsrc by an integer factor. The sum of
elements is conserved, unless the factor produces a remainder,
in which case the remainder values are truncated from vdest.
producesvec.indgen(1, 10, 1) vec1.rebin(vec, 2)
vec1:
3 7 11 15 19where each pair of
vec elements is added together into one element.
But,
adds triosvec.indgen(1, 10, 1) vec1.rebin(vec, 3)
vec elements and gets rid of the value 10, producing
vec1:
6 15 24
Vector
obj = vdest.pow(vsrc, power)obj = vsrcdest.pow(power)
Vector
obj = vdest.sqrt(vsrc)obj = vsrcdest.sqrt()
Vector
obj = vdest.log(vsrc)obj = vsrcdest.log()
Vector
obj = vdest.log10(vsrc)obj = vsrcdest.log10()
Vector
obj = vdest.tanh(vsrc)obj = vsrcdest.tanh()
Vector
obj = vdest.abs(vsrc)obj = vsrcdest.abs()objref v1 v1 = new Vector() v1.indgen(-.5, .5, .1) v1.printf() v1.abs.printf()
Vector
obj = vdest.index(vsrc, indices)vsrc indexed by the vector indices are collected
into vdest.
makesobjref vec, vec1, vec2, vec3 vec = new Vector() vec1 = new Vector() vec2 = new Vector() vec3 = new Vector(6) vec.indgen(0, 5.1, 0.1) //vec will have 51 values from 0 to 5, with increment=0.1 vec1.integral(vec, 0.1) //Euler integral of vec elements approximating //an x-squared function, dx = 0.1 vec2.indgen(0, 50,10) vec3.index(vec1, vec2) //put the value of every 10th index in vec2
vec3 with six elements corresponding to the integrated integers from
vec.
Vector
x = vec.min()x = vec.min(start, end)
Vector
i = vec.min_ind()i = vec.min_ind(start, end)
Vector
x = vec.max()x = vec.max(start, end)
Vector
i = vec.max_ind()i = vec.max_ind(start, end)
Vector
x = vec.sum()x = vec.sum(start, end)
Vector
x = vec.sumsq()x = vec.sumsq(start, end)
Vector
x = vec.mean()x = vec.mean(start, end)
Vector
x = vec.var()x = vec.var(start, end)
Vector
vec.stdev()vec.stdev(start,end)
Vector
x = vec.stderr()x = vec.stderr(start, end)
Vector
x = vec.dot(vec1)vec and vec1.
Vector
x = vec.mag()
Vector
obj = vsrcdest.add(scalar)obj = vsrcdest.add(vec1)vsrcdest.
vsrcdest and vec1 must have the same size.
Vector
obj = vsrcdest.sub(scalar)obj = vsrcdest.sub(vec1)vsrcdest.
vsrcdest and vec1 must have the same size.
Vector
obj = vsrcdest.mul(scalar)obj = vsrcdest.mul(vec1)vsrcdest either by either a scalar or the
corresponding elements of vec1. vsrcdest
and vec1 must have the same size.
Vector
obj = vsrcdest.div(scalar)obj = vsrcdest.div(vec1)vsrcdest either by a scalar or by the
corresponding elements of vec1. vsrcdest
and vec1 must have the same size.
Vector
scale = vsrcdest.scale(low, high)
Vector
boolean = vec.eq(vec1)
Vector
x = vec.meansqerr(vec1)x = vec.meansqerr(vec1, weight_vec)vec and
the corresponding elements of vec1. vec and vec1 must have the
same size.
If the second vector arg is present, it also must have the same size and the return value is sum of w[i]*(v1[i] - v2[i])^2 / size
Vector convlv fft spctrm correl filterThe following routines are based on the fast fourier transform (FFT) and are implemented using code from Numerical Recipes in C (2nd ed.) Refer to this source for further information.
Fourier
obj = vdest.correl(src)obj = vdest.correl(src, vec2)
Fourier
obj = vdest.convlv(src,filter)obj = vdest.convlv(src,filter, sign)t=0..t=n/2 followed by countdown t=n..t=n/2. The size of filter
should be an odd <= the size of v1>.
objref v1, v2, v3 v1 = new Vector(16) v2 = new Vector(16) v3 = new Vector() v1.x[5] = v1.x[6] = 1 v2.x[3] = v2.x[4] = 3 v3.convlv(v1, v2) v1.printf() v2.printf() v3.printf()
Fourier
obj = vdest.spctrm(vsrc)
Fourier
obj = vdest.filter(src,filter)obj = vsrcdest.filter(filter)
Fourier
obj = vdest.fft(vsrc, sign)obj = vsrcdest.fft(sign)If vsrc.size() is not an integral power of 2, it is padded with 0's to the next power of 2 size.
The complex frequency domain is represented in the vector as pairs of numbers --- except for the first two numbers. vec.x[0] is the amplitude of the 0 frequency cosine (constant) and vec.x[1] is the amplitude of the highest (N/2) frequency cosine (ie. alternating 1,-1's in the time domain) vec.x[2, 3] is the amplitude of the cos(2*PI*i/n), sin(2*PI*i/n) components (ie. one whole wave in the time domain) vec.x[n-2, n-1] is the amplitude of the cos(PI*(n-1)*i/n), sin(PI*(n-1)*i/n) components. The following example of a pure time domain sine wave sampled at 16 points should be played with to see where the specified frequency appears in the frequency domain vector (note that if the frequency is greater than 8, aliasing will occur, ie sampling makes it appear as a lower frequency) Also note that the forward transform does not produce the amplitudes of the frequency components that goes up to make the time domain function but instead each element is the integral of the product of the time domain function and a specific pure frequency. Thus the 0 and highest frequency cosine are N times the amplitudes and all others are N/2 times the amplitudes.
... //define a gui for this example
N=16 // should be power of 2
c=1 // 0 -> sin 1 -> cos
f=1 // waves per domain, max is N/2
setup_gui() // construct the gui for this example
proc p() {
v1 = new Vector(N)
v1.sin(f, c*PI/2, 1000/N)
v1.plot(g1)
v2 = new Vector()
v2.fft(v1, 1) // forward
v2.plot(g2)
v3 = new Vector()
v3.fft(v2, -1) // inverse
v3.plot(g3) // amplitude N/2 times the original
}
p()
The inverse fft is mathematically almost identical to the forward transform but often has a different operational interpretation. In this case the result is a time domain function which is merely the sum of all the pure sinusoids weighted by the (complex) frequency function (although, remember, points 0 and 1 in the frequency domain are special, being the constant and the highest alternating cosine, respectively). The example below shows the index of a particular frequency and phase as well as the time domain pattern. Note that index 1 is for the higest frequency cosine instead of the 0 frequency sin.
Because the frequency domain representation is something only a programmer could love, and because one might wish to plot the real and imaginary frequency spectra, one might wish to encapsulate the fft in a function which uses a more convenient representation.
Below is an alternative FFT function where the frequency values are spectrum amplitudes (no need to divide anything by N) and the real and complex frequency components are stored in separate vectors (of length N/2 + 1).
Consider the functions
FFT(1, vt_src, vfr_dest, vfi_dest)FFT(-1, vt_dest, vfr_src, vfi_src)The forward transform (first arg = 1) requires a time domain source vector with a length of N = 2^n where n is some positive integer. The resultant real (cosine amplitudes) and imaginary (sine amplitudes) frequency components are stored in the N/2 + 1 locations of the vfr_dest and vfi_dest vectors respectively (Note: vfi_dest.x[0] and vfi_dest.x[N/2] are always set to 0. The index i in the frequency domain is the number of full pure sinusoid waves in the time domain. ie. if the time domain has length T then the frequency of the i'th component is i/T.
The inverse transform (first arg = -1) requires two freqency domain source vectors for the cosine and sine amplitudes. The size of these vectors must be N/2+1 where N is a power of 2. The resultant time domain vector will have a size of N.
If the source vectors are not a power of 2, then the vectors are padded with 0's til vtsrc is 2^n or vfr_src is 2^n + 1. The destination vectors are resized if necessary.
This function has the property that the sequence
leaves vt unchanged. Reversal of the order would leave vfr and vfi unchanged.FFT(1, vt, vfr, vfi) FFT(-1, vt, vfr, vfi)
The implementation is:
execute following example
proc FFT() {local n, x
if ($1 == 1) { // forward
$o3.fft($o2, 1)
n = $o3.size()
$o3.div(n/2)
$o3.x[0] /= 2 // makes the spectrum appear discontinuous
$o3.x[1] /= 2 // but the amplitudes are intuitive
$o4.copy($o3, 0, 1, -1, 1, 2) // odd elements
$o3.copy($o3, 0, 0, -1, 1, 2) // even elements
$o3.resize(n/2+1)
$o4.resize(n/2+1)
$o3.x[n/2] = $o4.x[0] //highest cos started in o3.x[1
$o4.x[0] = $o4.x[n/2] = 0 // weights for sin(0*i)and sin(PI*i)
}else{ // inverse
// shuffle o3 and o4 into o2
n = $o3.size()
$o2.copy($o3, 0, 0, n-2, 2, 1)
$o2.x[1] = $o3.x[n-1]
$o2.copy($o4, 3, 1, n-2, 2, 1)
$o2.x[0] *= 2
$o2.x[1] *= 2
$o2.fft($o2, -1)
}
}
...
N=128
delay = 0
duration = N/2
setup_gui()
proc p() {
v1 = new Vector(N)
v1.fill(1, delay, delay+duration-1)
v1.plot(g1)
v2 = new Vector()
v3 = new Vector()
FFT(1, v1, v2, v3)
v2.plot(g2)
v3.plot(g3)
v4 = new Vector()
FFT(-1, v4, v2, v3)
v4.plot(g4)
}
p()