Investment Studio > Expressions > Functions > DSP > MEM

float array[1][poles + 1] mem(float array data, integer poles)

Returns a row vector {a0, a1, ..., apoles} containing the 1 + poles coefficients of the order poles MEM (Maximum Entropy Method) estimate of the data array's spectrum (see formula below). The result can be used with memdom, memf and memp. It can also be used to obtain linear predictions of future data elements.

The number of poles (the order of the MEM estimate) should be >= 1 and < the number of elements in data (but will be automatically capped if larger). In practice, it should be limited to a few times the number of expected sharp spectral features (more poles will always yield more peaks, ultimately by splitting actual peaks).

All elements in data are converted to floats (with exclusion if conversion fails) and interpreted as a single, one-dimensional sequence. If the data array is two-dimensional, it's read in the usual order (i.e. row by row: left to right, top to bottom).

Properties of the MEM spectral estimate

In principle, the amplitude spectrum of a function is the amplitude of the function's Fourier transform. When dealing with a sequence of sampled values, the amplitude spectrum of the underlying function can be estimated by computing the amplitude of the sequence's Discrete Fourier Transform (DFT), as done by fftp. But this is just one out of an infinite number of possible estimates, and not necessarily the best.

The Maximum Entropy Method (MEM) is an alternate procedure for obtaining a spectral estimate from a sequence of sampled function values. In mathematical terms, the power (= amplitude squared) spectrum P(f) is decomposed as

          a0    
P(f) =  ¾¾ ¾¾ ¾¾¾ ¾¾ ¾¾¾¾¾¾ ¾¾ ¾
  ½   poles   i (2 p / N) k f ½ 2
  ½ 1 + å ak e   ½  
  ½   k = 1     ½  

where N is the number of samples in the input sequence, e is Euler's number and i is the imaginary unit.

The main advantage of this decomposition over the DFT is that it doesn't depend on a fixed grid of equally spaced frequency values. The DFT samples the spectrum at integer multiples of the input sequence's fundamental frequency; the MEM fits the dominating peaks at their actual locations, while smoothing the rest of the spectrum. This makes it especially suited to finding and fitting sharp spectral features, corresponding to well-defined periodic patterns in the data.

There are several disadvantages, too. MEM power spectrum estimation is slower than DFT estimation, since it can't take advantage of the FFT (Fast Fourier Transform) algorithm. Large numbers of poles and/or data points can lead to significant roundoff errors. Data with very sharp spectral features can lead to artificially split peaks even at low orders, and noisy data can induce spurious peaks, especially for large orders. When possible, peaks found with MEM should therefore be verified with more conservative methods, like an ordinary DFT.

MEM coefficients can also be used for so-called linear prediction filtering. The scalar product of the last poles elements in data and the vector {apoles, ..., a1} is the linear prediction of the next element in the data sequence. In this context, a0 is the mean square error of the fit; its square root can therefore be used as an estimate of the prediction error. A predicted sequence can be obtained by iteration: pop data (move back all elements one position, discarding the oldest one) and append the predicted element at the other end, then take the scalar product to get the next predicted element.

Example

The 1024-point sine wave with normalized frequency 6 / 100 = 0.06 and amplitude 1,

_data = mop("sin()", makevector(1024, 0, 12 * PI / 100))

has the 10-pole MEM coefficients

_mem = mem(array(_data), 10)

To obtain the actual 10-pole spectrum, we feed the MEM coefficients to memf:

_memf = memf(array(_mem), FALSE, 0, 1/1024, 513)

where we have chosen amplitude (rather than power) output in order to facilitate comparison with the amplitude spectrum returned by fftp:

_fftp = fftp(array(_data))

Putting it all together, the graph sources

=index(array(_data), x, 1)

=index(array(_memf), x + 1, 1)

=index(array(_fftp), x + 1, 1)

look like this (black lines for input sequence, blue dots for MEM spectrum estimate, gray crosses for DFT spectrum estimate):

DFT and MEM are in excellent agreement about the position of the amplitude spectrum's peak, but note how the the DFT spectrum is smeared out over many bins, while the MEM peak is effectively all in a single bin.

We can pinpoint the peak using memdom:

=memdom(array(_data), FALSE, 0, 0.001, 101)

returns {0.06, 1.70741234308845, 99.5048103133754}. This tells us that the spectral maximum in the normalized frequency range [0, 0.1] occurs at normalized frequency 0.06 (± < 0.001, with the chosen step length), with an amplitude of » 1.71, and a signal/noise ratio (max amplitude over average amplitude in the scanned frequency range) of » 99.5.

See also fftc, fftp, memdom, memf, memp.