EXAMPLES

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()

spctrm

Fourier

SYNTAX

DESCRIPTION

filter

Fourier

SYNTAX

DESCRIPTION

fft

Fourier

SYNTAX

DESCRIPTION

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.

execute following example

...	//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

SYNTAX

DESCRIPTION

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

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()

SEE ALSO

neuron/general/classes/vector/vect2.hel : 20988 Oct 18