home *** CD-ROM | disk | FTP | other *** search
-
- 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.
-
- -------------------
-
-
-