The Fred Fish Collection 1.5/ffcollection-1-5-1992-11.iso/ff_disks/300-399/ff364.lzh/DPFFT/readme

Text File | 1990-08-12 | 16KB | 374 lines

This code is an extension of the program on FISH-disk 322(?) and is
considered as freeware in the usual sense.
The following has been added:

    -A serial correlation
    -An EWK (=Einstein-Wiener-Khinchine) transform
    -Normalised-linear and log-log spectra
    -A 3-D display of various FFT-spectra

The 3-D display is on the basis of FFT. The EWK-transform is too slow.
The main reason for the use of the latter is the window-closing principle
described in Ch.7 of G.M.Jenkins and D.G.Watts, Spectral Analysis,
Holden-Day, 1968. For more information, see below.

Any suggestions or information about scientific software on the Amiga
is appreciated. I am particularly interested in a comfortable and
interactive display of theoretical concepts and of experimental data.
This program is used by me as a tool for a characterisation of animated
sequences by means of topological invariants.

            A.A.Walma
            Ziegelmattenstrasse 5
            7800 FREIBURG
            W-Germany

    --------------------------------------------------------------------

                               DPFFT 2.2

                          (c) 1990 Alle Anne Walma


1.DATA-DISPLAY
  -----------

    The program is started either from CLI by invoking "DPFFT <cr>"
or from the Workbench, where the icons are selfexplanatory. Double
klicking in the main window allows you to change the screen to PAL height.

    DPFFT works within an integer-range of 0-10000. This
suffices for the usual 8-Bit and 12-Bit AD converters.
Negative binary values are accepted (sound-file for instance) and shown
in a different color in the plotwindow (vertical scale).
Two examples of experimental data are included:

    - eeg    = A sample of a 12-Bit brainwave digitalisation (1500
                   data at a sampling rate of 250 Hz)
    - sinus  = A sinusoid with 10 periods and 100 samples per period

    The filerequester asks you for the total number of data to be read in.
With "starting index" is meant the first datum in the display (default
is one). If you have lots of data and not so much RAM, you may want to start
at 78000 for instance and put in 90000 for the "total number". RAM is only
needed for 12000 then.
With "redundancy" is meant the number of points you wish to ignore in case
of a very high sampling rate. Two points means that every third point is
taken into account.

    NOK (="Not OK") brings you back to the main menu and with OK the first
datasamples are being displayed. The display is "pixeloriented". In
other words, each horizontal pixel represents a sampled value and each
vertical pixel an amplitude-unit. Reading in "eeg" shows what this
means since this signal varies in real amplitude from ca 2000-7000, whereas
the displayheight amounts to 200. You can measure this by rightklicking
the mouse if the cursor is in the displayfield ("klick and move").


The displayed number, in the title bar left, indicates the degree
of horizontal expansion. One stands for one horizontal pixel per datapoint.
The triangular gadgets at the right of it allow a modification of the
number of pixels per point.

With the horizontal arrows you can "page" through the data. A leftklick
on the next gadget changes the window allowing a faster "paging". More
to that, see below. RST means a reset of the data. If they disappear,
the horizontal and even more so the vertical potentiometer may bring
them back. Keeping the left-mouse button pressed and moving the
slider or just  single klicks in the proportional gadgets will make this
clear to you.

Other gadgets appear by clicking twice the right mousebutton.
The first two (triangular) gadgets can be used to compress the amplitude
of the data. This  maybe useful for very large amplitudes. More useful
is the next gadget. The data are now centered according to the largest
and smallest available amplitude. All these actions are taken into
account by the value in the titlebar ("klick left on OK and klick
right in the window") indicating the vertical amplitude.

With PRT a rastportdump can be made if a printer is connected.
A small requester is fired up. For printers with 8 needles
RDC reduces the picture in size by about a half. NOK allows
a retreat without action.

There are a number of ways to reduce the size of the printout.
In the first place there is RDC of course. Going back to the main
window, however, the size of the window can be custimized by
doubleklicking the right mousebutton and changing the windowsize
in the appearing requester. Finally, you may wish to use
preferences (WorkBench 1.3).

With NOK you leave the plotwindow and return to the main window.
OK lets the requester disappear.

The inactivated gadgets are used in the expanded window where a fast
"paging" becomes possible. To that purpose the two horizontal
bars in the titlebar are klicked with the left mouse-button. De
abovementioned gadgets are now activated and using them, more data
will be  displayed. The more data,the larger the "page" and the faster
you will get an overall picture by means of the horizontal arrows in the
titlebar. If you have 4 fields of data in the window for instance, one klick
shows the next 4*600=2400 data. This can be checked by removing the
gadget (klick on OK) and klicking once with the right mouse-button somewhere
in the display field ("klick and move"). The display in the titlebar
shows you where you are.

You can freeze the vertical line somewhere by klicking the
right button one more time. Going back to the plotwindow (activate
horizontal bars in titlebar) shows that the plotted data start at this
chosen position. For an accurate datacut the freezing should be carried out
carefully (no fast moving around and keep the button pressed for a second).


                            ------------

2.FOURIER (FFT)
  -------------
The general idea of the program is to read in arbitrary experimental
data, visualize them and carry out a Fouriertransform of an interesting
detail. The start of the transform is determined by the first datum in
the plotwindow.

The defaultvalue of 6 in the FFT requester means that 2^6 points will be
used for the transform. The reason behind the power of two is simply
velocity. Clicking on OK ,using default, yields 2^(6-1)=32 harmonics.
With redundancy you can skip a variable number of points for the
transform in case of a long record. Typing in 1 means that every third
point is taken into account for the transformation.

The appearing window shows in the upper half the amplitudes and
below you see the appropriate phases.
A double-klick with the right mousebutton brings up the requester.
The arrows manipulate the Amplitudes of the Harmonics. OK brings you
back into the main window and with NOK into the plotwindow. With PRT
a little requester is fired up where RDC means a reduction of the
printed picture (can be useful for 9 needles). The crosshair for the
amplitudes (leftklick) can only be activated without the requester.
IDCMP input of Intuition is blocked in case of a request. To that
purpose the CRH gadget is activated which lets the requester
dissappear.


EXAMPLE
-------

To see the full power of a fouriertransform as far as periodic signals
is concerned (the interpretation is more involved for random signals),
the following step by step procedure may be useful (for the usage of the
program as well):

1.

If you are on the workbench, double-klick with left on the dpfft icon.
(If you are in CLI, type <dpfft> and return). It is possibly simplest to
boot with the FISH-disk itself so you can use the default values of dpfft.

2.

Activate with the right mousebutton the ASCII-file option in the menu

3.

Search for and Klick on signals(dir) and after that on "eeg"

4.

Read in, say, 1024 data (leftklick on the 'Number of data:' stringadget
and type in 1024)

6.

Left-klick with the mousebutton on OK and choose the 'plot' option in
the menu

7.

Double-klick the right mousebutton to fire up the requester for the options
Center the data with the second gadget from the left. Observing the signal
(expand it horizontally, for instance, by means of a leftklick on OK and
using the small triangles in the titlebar) shows that it is not possible
to detect hidden periodicities in the signal

8.

Let the requester reappear by means of a double-klick with the right
mousebutton and activate  FFT

9.

Since we read in 1024 data, the maximum N that can be used is 10.
(The program works with max 2^10=1024 points since the windowresolution
amounts to 640 sothat 512 harmonics can be displayed maximally, if we want
to retain a pixel to pixel resolution). Klick on the stringadget and
type in 10.

10.

With no window or prewhitening (=first difference filter), which are
the WND and FLT gadgets, a direct transform is carried out by activating
OK. This takes a few seconds. A transform with N<10 will be faster of course.

11.

Double-klick with the right mousebutton for the requester. Use the
arrow which is pointing downwards for reducing the amplitudes of the
harmonics. Do this until all the harmonics are displayed. A pronounced
peak should appear at right.

12.

Let the requester dissappear by means of CRH and klick once with the
right mousebutton in the harmonics area. A crosshair appears. Go to
this peak and note the number of the harmonic (which should be 206).



DISCUSSION
----------

The pronounced peaks are the hidden periodicities in the signal.
Remember that the recordlength of 1024 points with a 4 ms sampling period
totals to ca 4 seconds. The question is as to whether the peak at right
belongs to the signalsource itself or not. Placing the crosshair on that
particular harmonic shows that it carries the number 206. It is now clear
where it comes from, since 204/(4 seconds) represents the frequency of
the linevoltage (50 Hz in Europe). The other peaks can be identified
with the socalled alpha and beta waves of an electroencephalogram.


                    ---------------------


3. THE EWK-TRANSFORM
   -----------------

This one is based on a Fourier-transform of the autocorrelation function
of the data. I prefer this one for smoothing reasons. The window-closing
lets you choose yourself the degree of smoothing of the spectra. The
socalled "window-carpentry" (Parzen, Welch, Bartlett, Hanning etc) is
not so flexible.


A SAMPLE SESSION

We know now how to FFT the eeg signal, so the final results can be
compared.

1.  Read in the eeg again (you can change the default path
    by double klicking the right mousebutton in the main window)

2.  Choose "plot" in the menu and double klick for the control gadgets

3.  Use the FFT gadget and choose the correlation image

4.  Fill in the number of points which should be at least 200 pts
    if you want a Fourier transform afterwards, since it is worked
    with 100 frequencies (why this is so, see further)

5.  A serial correlation takes time (though it has been coded in assembler)
    Since the calculation time goes up quadratically, a large number
    of points should be avoided (not of course, if a long transform is needed).
    The appearing window has been scaled to +/- 1000, sothat a value of
    340, for instance, actually represents a correlation of 0.34. You can
    check this by klicking once with the right mousebutton in order
    to obtain the crosshair. The arrows of the controlgadgets (double klick)
    allow a horizontal expansion.

6.  Choose the FFT gadget. You are asked now for the number of lags. This
    means that you force the autocorrelation to zero over that number
    in a linear way. This is the Windowclosing idea (See Jenkins and Watts
    Ch.7). The more points you choose ,the more detailed the spectrum
    will be. Type in 50 for instance. Choose the LIN gadget.

7.  Again you have to wait. A real cosinetransform is carried out (no FFT).
    The spectrum is amplified by the proportional gadget at right with
    automatic scaling. Compare the obtained form with a FFT spectrum
    for N=8 (128 freq. vs. 100 freq. here) for instance.

8.  Get the requester back by a doubleklick and choose LOG. If you want
    to pull the spectrum down, go back to the linear window (double klick
    again) and reduce the amplification. Note that the harmonic component
    has now the dimension of a real frequency, in contrast to the FFT.



DISCUSSION
----------

If you are not only interested in a certain harmonic but also in the spectral
power as such, the FFT should be used. The cosine transform is normalised.
I did this because I wanted to compare different signals with different
variances. The relative power (power/variance) is then more conclusive.
For all sorts of signals, the area below the powerspectrum will then be
the same (=0.5, if correct).

Since it takes a lot of time for a long correlation, the final result is put
in a file on the root directory (double klick in main window for the
path) named "corr.dat". If you want to have a look at this file
(with ffp data), choose in the main menu the float option. Getting
back to the main window is always possible by klicking on NOK-NOK-....

The choice of 100 freq for the display has a special reason. I am interested
in relatively low frequencies ( < 100 Hz). For a 250 data record with a
sampling rate of 4 ms, the recordlength is one second. The max. freq. of
the transform is then the 125. harmonic, which means in this case 125 Hz.
The 100. harmonic is then 100 Hz and the first harmonic 1 Hz.
For a data record with 2500 points, this will be 10 Hz and 0.1 Hz respectively.
The trick is therefore to modify the frequency-range with the length of
the data record.


                        ------------

4. 3-D SPECTRA
   -----------

This is a timeconsuming process. For a first look take the defaultvalue
N=6. This means that 2^(N-1)=32 spectra wil be displayed. To achieve this,
you have to klick on the "3-D" image of the FFT requester. In detail:
   -Read in 1000 data of the EEG signal for instance
   -Choose "plot" in the main menue
   -Fire up the control gadgets (double klick with right)
   -Choose FFT
   -Now you have the FFT requester
After choosing the "3-D" Image a small requester appears. "Step" represents
the number of points between the subsequent FFT spectra. All spectra are
smoothed with the Bartlett window. Alpha and beta are numbers between 0-360
with which you are able to rotate the picture. Just use the defaultvalues for
the moment and activate OK.
    With a double-klick another requester comes up. Change first the
value of alfa, for instance. This represents a rotation about the z-axis.
With the  rdce gadget the rotation about the x-axis (which is beta) can
be reduced. There are cases where a value of 1 for beta is already too large.
The image is scaled in such a way, that the largest amplitude
automatically covers all vertical available pixels. It therefore depends
on the spectrum whether a full rotation about the x-axis is possible.
Particularly with a large low power (such as is the case with the eeg),
the rotation with beta is restricted.

EXAMPLE
-------

1. Read in the sinus signal

2. Choose "plot" in the main menu

3. Fire up the controlgadgets

4. Choose "FFT"

5. Activate the "3-D" image

6. Use a "step" of ten and klick OK

Since the sinusoid has 100 samples per period, and the default value for
the FFT is N=6, only part of the sinus is transformed, which explains
the sidebands. The sinusoidal appearance of the 3-D picture is due to the
fact that a running transform depends on the first and last data-point
of the record and these evolve in a sinusoidal way.

                        -------------------