Manual Matlab scripts Laboratory of Fluid Mechanics

15 juni 2005

This manual is combination of all manuals made for different Matlab scripts. The setup may differ per script. An overview of the usage has been given in Usage.doc.

Table of Contents

Manual Matlab scripts Laboratory of Fluid Mechanics. 1

Asc2wav, Converting Ascii files and using the DELFT-AUKE program Waves. 2

Decomp, Separate spectrum of irregular wave in incident and reflected part. 3

Calculate distance of wave gauges needed. 3

Program Decomp. 4

Usage of the program: 6

How to handle if a file cannot be read. 8

Lvmview, view measurements. 9

Refreg, calculation of the reflection of a regular wave. 11

Short description of the program Refreg: 13

Program used to test Refreg. 15

Things to be investigated. 15

Additional program points in the future. 16

Refreg and Spect, remarks about spectra. 17

References. 19

Asc2wav, Converting Ascii files and using the DELFT-AUKE program Waves.

If you want to run Waves after decomposition, run Decomp to calculate a time signal containing the incident and reflected wave. Otherwise use the original measurements as input in the conversion program.

Run the Matlab program Asc2wav; this program converts the Ascii-file to a .dat-file and a .seq-file, which can be used as input in the DELFT AUKE program Waves, and a .pcf-file containing commands to instruct Waves.

De header in the ASCII-file produced by Decomp is 11 lines long.

==============================================================

REMARK: the scale factor from measurement to meters is:

1 if the output of Decomp is used as input in waves,

depending on the settings of the wave gauges if the measurement file is input to waves.

==============================================================

On the Windows98-pc with DELFT-AUKE/Process:

Copy *.dat, *.seq, and the .pcf-file to a folder where you want to store the results; make a file Waves.bat with line

waves.exe in=name1.pcf ou=name2.out er=waves.err

name1 is the name of the .pcf-file: waves.pcf if you confirmed to overwrite the file, otherwise: the name of the .dat-file,

name2 is a name chosen by yourself, e.g. the name of the .dat-file.

If wanted, change the percentages of exceedance, specified in the lines

exceedance,perc=...

27 April 2005

Decomp, Separate spectrum of irregular wave in incident and reflected part.

18 April 2005

Calculate distance of wave gauges needed

To be able to separate the waves into incident and reflected wave, 2 or 3 wave gauges are needed (3 is preferable). The signals from the wave gauges lead to a system of equations, which will enable the calculation of incident and reflected wave. If the distance is a multiple of the wavelength, the system of equations will be singular.

If three wave gauges are used, distances of 0.3 m between gauge 1 and gauge 2, and 0.4 m between gauge 2 and gauge 3 will do in most cases, where gauges 1 is supposed to be nearest to the wave board. This will work if the peek period is between 1 to 3 seconds.

If two wave gauges are used, the program Distance of WL | Delft Hydraulics can be used:

DelftAUKE, program Distance (DOS), 2 wave gauges.

Input to Distance:

T_p peak period

h water depth

Output:

f_min- f_max frequency range of valid solution

L_min – L_max wave lengths according to frequencies at given water depth

dx_low - dx_highrange of distances between which a valid solution is possible.

The minimum distance of the wave gauges specified in the manual is 0.2 m. If this distance is not in the range found by the program Distance, specify 0.2 m as input to the program Distance. In that case Distance calculates the frequncy range of valid solution if the distace is 0.2 m.

Measurements -> Dasylab file

Input parameters to Decomp (input from .txt-file):

About measurements:

parameter	description
dt (s)	time step used in data file
T_p (s)	approximated peak period (only used to calculate the frequency resolution in the variance spectra to be calculated)
f_res	desired frequency resolution of spectrum/peak frequency
f_min- f_max(Hz)	frequency range where the spectrum of incident and reflected wave must be calculated (f_max=-1 to calculate this interval automatically from the distances in case of two gauges, and takes the default 0 – 2.5 Hz in case of 3 gauges). Some times, f_max=-1 will select a part that doesn’t cover the whole Jonswap spectrum. Decomp will show a message box on the screen in that case.
thresh	threshold to determine the part of the variance spectrum where the decomposition is valid, E>=thresh • E_max. Most of the times, the frequency range of an ideal Jonswap spectrum with the same threshold will extend this part, caused by noise in the calculated spectrum. The program will print both ranges in the output file.
col_t	column number in the data file where the time has been registrated(0 if not registrated); col_t is used by the input function only; it will not be used by this program, but it is needed to locate the colum numbers of the wave gauges correctly.

about channels to be analysed:

column number

position of gauge (m)

scale from registered data value to meters

about the water depth as a function of the position:

x (m), h (m) position and water depth in at least P+1 points, where P is the number of wave gauges

Program Decomp

Method: based on method of Zelt en Skjelbreia (1992)

Matlab version: Henk Jan Bakkenes (2002)

Suppose N points, measured at time t=n∆t, n=1,2,…N. In this case a wave can be described by N spectral terms.

The data registered by wave gauge p can be represented by a sum of incident and reflected wave:

(1)

where a_i,j and a_r,j are complex (i.e. containing the phase shift of the incident and reflected wave respectively). The mean value of η is zero, so j=0 disappears.

Furthermore, ω_j=2π/T, whereT=N∆t.

Standard Fourier analysis of η(x_p,t ) gives:

(2)

Combination of (1) and (2) -> system of P linear equations and two unknown complex values (a_i,j and a_r,j) for each spectral component j (P : number of gauges).

A least-squares solution can be calculated using a standard Matlab operation.

Equations must be well conditioned: it is not allowed to put wave gauges in points of equal phase.

To calculate the distances needed, a DOS program Distance is available (from the Delft Hydraulics package DELFT-AUKE).

Steps of Decomp:

(1) Ask which output the user desires (see list on p.5).

(2) Read data file, parameter file, and file containing profile of water depth.

(3) If asked by user: calculate frequency range of valid decomposition according to the distances of the wave gauges. Otherwise, use the values supplied by the user.

(4) Ask part of input data to be analysed (begin and end point).

(5) Calculate phase shift at positions of the wave gauges, and if needed, calculate the change of the amplitude caused by changes in water depth.

(6) For each wave gauge, calculate FFT A , variance spectrum E, and averaged spectrum Es, averaged over frequency bands of length specified by input.

Points (7) through (9) describe parameter calculations per wave gauge. The parameters are chosen according to a program Reflecmf from the DELFT-AUKE package (manual July 2000).

(7) From averaged variance spectrum: peak period T_p, period T_pD in middle of connected part of the spectrum where Es >= 0.8 Es_max , and the frequency range of a connected part of the spectrum where Es >= thresh • E_max. This frequency range will be used to calculate the parameters in points (8) and (9).

(8) From "raw" variance spectrum in frequency range Es >= thresh • E_max: spectral moments (m₀, m_2,and m₄).

(9) From spectral moments: mean period between successive zero crossings and between successive maxima, and correlation between successive wave heights.

Continuation of program Decomp:

As a valid frequency range, the frequencies must be valid for the used distances (point (3)), and the spectrum must be above the threshold (point (7)), otherwise the equations will be of the type

0.x = 0.

Most of the times, an ideal Jonswap spectrum will have a wider frequency range above the threshold. This range will be given in the output as well. Compare the range used by Decomp with the frequecy range of an ideal Jonswap spectrum to evaluate the results.

In the figures, more frequencies than the above range are shown. This will help the user to inspect if the valid frequency range contains a good representation of the signal. Vertical lines with arrows are used to show which part of the spectrum is valid.

If the coefficient matrix of the equations is nearly singular, the program will stop, and show an error message. Sometimes the conditions in the equations are reasonable, but the spectrum will show an unexpected peak. If this peak is in the valid part of the spectrum, the distances may be chosen incorrectly. Check if the values are in accordance with this program.

More information about distances: DELFT-AUKE Process Document (2000).

(10) Calculate amplitude spectra of incident and reflected wave from FFT-values A for each spectral component.

(11) Variance spectrum and spectrum averaged over frequency bands of incident and reflected wave in specified frequency range.

(12) Spectral parameters from "raw" variance spectra of incident and reflected wave: moment m₀, significant wave height, and mean factor of reflection.

(13) Reconstruction of the incident and reflected wave as a time series by inverse FFT. Only the "valid" part of the spectrum will be used for this computation.

Parameters calculated in steps (8) and (9):

parameter	description
	spectral moments
T_0,1= m₀/m₁	average period
	average period between two zero crossings
	average period between two successive local maximum values
	narrowness parameter
	broadness parameter
κ	correlation value, defined in DELFT-AUKE manual (2000), p.4-6
γ	coefficient of linear correlation, defined in DELFT-AUKE manual (2000), p.4-7 as ρ_HH
H_m₀ = 4√m₀	significant wave height

Usage of the program:

N.B.: This program uses dialog boxes. In Matlab 7, dialog boxes will not be deleted immediately after clicking OK, but the program will continue; please don’t click OK again.

Input files:

· data file (Dasylab file),

· parameter file: .txt file

Contents of parameter file:

time step used in data file

0.02

Tp (s), f-resol/Tp fmin (Hz), fmax (Hz), thresh

1.27, 0.025, 0.25, 1.65, 0.005

column with time (0 if none)

columno gauge, pos (m), scale fact to m

2, 0, 0.025

3, 0.25, 0.025

Line number 1, 3, 5, and 7 are text lines; the other lines are data lines; if more values are stored in one line, the values must be separated by commas. The values are examples, and have to be replaced by the user. If fmax = -1, fmin and fmax will be calculated by the program (2 gauges), or replaced by 0.02 and 2.5 respectively (3 or more gauges).

· profile of depth

Contents of file containing profile of depth:

pos (m), depth (m)

0.0 0.6

0.125 0.6

0.25 0.6

0.5 0.6

Line 1 contains text; the other lines contain data, separated by blanks or tabs. The values are examples.

If P gauges are used, at least P+1 data lines are needed (this restriction has not yet been checked).

N.B.: In the file containing the depth profile, spaces or tabs ar allowed as a separation. In the parameter file, commas are allowed only, since the data are read in a different way. The program will send a message if a comma is missing.

Usage:

· Open Matlab and go to the folder where Decomp and subfolder .\Functions are stored.

Type Decomp.

Remark: the user must make the parameter file and the file containing the profile of depth in advance. The contents are described above.

· Desired output:

The program shows a list of possible outputs:

plot of raw variance spectrum of original time signal

plot of averaged variance spectrum of original time signal

plot of averaged variance spectrum of incident and reflected wave

plot of reflection factor

file of statistical results (spectral moments, refl. coeff)

file of variance spectra of inc. and refl. waves

file of time signal of inc. and refl. wave

For each item of the list, fill in 1 for yes or 0 for no.

· The input files are asked from dialog boxes.

· The time interval to be used in the calculations is asked by the program. A plot of the measurements will be put on the screen, which will help the user to select suitable data.

The begin time and end time in seconds must be given in a dialog box.

Close the figure to start the calculation.

· If the frequency range selected in the input file covers less than 90% of the frequency range of an ideal Jonswap spectrum with peak period TpD (mean value of all wave gauges), a message box will be shown on the screen.

The frequency range of an ideal Jonswap spectrum will be printed to the output file of statistical results. The algorithme to calculated this ideal Jonswap spectrum has been copied from Ap van Dongeren, WL | Delft Hydraulics (see function FN_jonswap).

· The desired figures will be produced, and dialog boxes will ask the names of the desired output files. If the spectra of incident and reflected waves show peeks, the equations probably have a singularity at those frequencies.

The incident and reflected wave as time series are calculated using the inverse FFT of the spectra of incident and reflected wave. Only frequencies for which the spectrum has been detected as "valid" will be used for the reconstruction. If the output file of time series of incident and reflected wave is asked by the user (see "Desired Output"), these series can be used as input for the DELFT-AUKE program Waves, which determines special wave parameters. A special Matlab script Asc2wav has to be run to make the files readable by Waves (and other DELFT-AUKE programs). Information from the programmer of Stevin III, room 0.07, tel. 85974.

How to handle if a file cannot be read

The next errors may occur:

If the number of header lines specified in case of a non-Dasylab file is less than the actual number, the read function used in the scripts gives a message. However, in some cases Matlab gives the next error:

Number of colums in line -- of ASCII file ... must be the same as the previous line.

In such a case, remove all spaces, signs, or numbers from the first position of each header line.

If the number of header lines specified is more than the actual number, no message will occur, but less data lines will be read.

If a Dasylab file gives an error like above, check if the right file has been chosen. If yes, copy this file, delete the first line, and read the copy specifying the remaining number of header lines.

The maximum number of samples to be read is exceeded. This number is set to 500000/nchan,

where nchan is the number of channels used. If more data are needed, change the value of maxnum in \Functions\FN_indata7, line 86, or read the file in parts, using different sample numbers. Sample numbers to be read can be specified via the dialog box:

read all samples?

Specify No in this case.

In a Dasylab file, it is preferable to store the time in decimal format. The time format with a ':' in Dasylab is not read as a time in Matlab.

Lvmview, view measurements

March 16, 2005

The Matlab function Lvmview has been developed to view ASCII files containing measurement data. If times are recorded, the recommended format is decimal. The ":"-format may differ from the Matlab time format (e.g. Dasylab files have a different time format).

Specifications to be fulfilled by the ASCII files are given in Help, General, or in the document Usage.doc. Some incidental occuring error message are discussed at the end of this manual.

The data can be viewed directly in the program Lvmview. After closing Lvmview, the data are kept in the workspace, and commands as mean(A) and plot(A(:,2)) can be given.

The use of Lvmview will be explained here.

1. Start Matlab.

2. Go to folder D:\Matlab_progs\Labprogs.

3. Type: Lvmview.

4. In the menu bar, click on File, followed by Open.

5. The box Input File is shown; find the file to be viewed, and select it.

6. The box Read all data? will be shown. In most cases, yes is the best choice. Choose no in case of a file of 10 MB or more only.

7. The number of header lines will be asked. Specify 0 if the input file has been created by Dasylab, because the program detects the end of the header of Dasylab files.

If no is specified in pint 6, begin sample, end sample, and step size are asked as well. Lvmview will send a warning if there are too many data to be read in once.

8. The parameter screen appears. Fill in the values asked by the screen:

· time step in the data file (this is not read from the header, because all types of ASCII files can be read),

· number of the column where the time has been stored (0 if no times are written to the data file, or if ":" is present in the time format).

9. The window channels appears. Select the channels to be plotted by clicking the mouse and holding the <control> key.

10. A plot of all samples of the first and the last channel will be shown, combined with a dialog box asking the time interval to be plotted on the screen. The times can be read from the horizontal axes. Fill in the desired interval, and click on OK. The samples chosen are displayed on the next screen. Close the figure to continue.

11. Click on Plot on the menu bar. Choose one of the next possibilities:

· plot channels separated in subplots, 3 channels per figure,

· plot all channels in one figure,

· choose other samples,

· choose other channels.

12. If one of the plot options has been chosen, one or more figures appear. The figures can be edited (see the manual "Using Matlab Graphics"), and/or printed. If wanted, the figure can be saved by clicking File, save as on the menu bar of the figure, and afterwards the figure can be opened by clicking File, Open on the menu bar of Lvmview.

If another time interval must be viewed, close the figures and press Plot, option samples on the menu bar.

13. After closing the figures, File, Open can be used to open a new file. This can be a figure (.fig), or an ASCII file. A .fig file always will be opened as a picture, and all other files will be opened as a data file. The old file will be kept in the workspace as long as no new file has been read, even after closing Lvmview. At the next restart of Lvmview, the file last read can be plotted again without reopening it.

To close Lvmview, always click the upper right cross, or click on File, Exit. Before closing Lvmview, check if all figures are closed (the figures are not closed automatically).

The next errors may occur:

If the number of header lines specified in case of a non-Dasylab file is less than the actual number, the read function used in the scripts gives a message. However, in some cases Matlab gives the next error:

Number of colums in line -- of ASCII file ... must be the same as the previous line.

In such a case, remove all spaces, signs, or numbers from the first position of each header line.

If the number of header lines specified is more than the actual number, no message will occur, but less data lines will be read.

If a Dasylab file gives an error like above, check if the right file has been chosen. If yes, copy this file, delete the first line, and read the copy specifying the remaining number of header lines.

The maximum number of samples to be read is exceeded. This number is set to 500000/nchan,

read all samples?

Specify No in this case.

Refreg, calculation of the reflection of a regular wave

April 02

Feb 03

June 05

To calculate the reflection of a regular wave, a Matlab program Refreg has been written in the Laboratory of Fluid Mechanics. The method used has been described by Goda and Suzuki (1976), see Goda (1985). In this method two wave gauges are used at a distance of about one fourth of the wave length.

Basic equations in the case of a regular wave with wave gauges at positions x=x₁ en x=x₂ are:

(1a)

(1b)

where

η the water-surface elevation relative to the mean water level,

t the time,

a_i,n , a_r,n the amplitude of the n-th harmonic of the incoming and the reflected wave,

k_n the wave number of the n-th harmonic,

ω_n the angular wave frequency of the n-th harmonic,

φ_i,n,φ_r,n the phase of the n-th harmonic of the incoming and the reflected wave.

In the Refreg program, the first harmonic is used only. Higher harmonic components and free or bound harmonics are not taken into account.

Equations (1a) and (1b) for the first harmonic are

(2a)

(2b)

(2a) can be written as

(3a)

In the same way, (2b) can be written as

(3b)

where

(4a)

(4b)

(4c)

(4d)

(4a) through 4(d) lead to the complex equations:

(5a)

, (5b)

where i=√-1.

(5a) and (5b) can be written as matrices:

(6)

The A and B in the right hand side of (6) can be found from a harmonic analysis of η(x₁,t) and η(x₂,t) in (3a) and (3b), e.g. by using a Fast Fourier Transform (FFT). In the program Refreg two zero crossings with the same sign, one at the begin and one at the end of the first data series, are used to determine the length of the series to be analysed. In that case the data series can be regarded as cyclic. The only error is a cut off error if the wave period does not fit on the time step. The FFT of Matlab is used on the two data series from the wave gauges under consideration, where the number of points used fits to the time between the zero crossings as meant above. The period with the maximum modulus of the FFT-coefficients is used as the base period.

Short description of the program Refreg:

N.B.: This program uses dialog boxes. In Matlab 7, dialog boxes will not be deleted immediately after clicking OK, but the program will continue; please don’t click OK again.

N.B. After closing a (sub)menu, allways check the number of Matlab Windows (see Matlab icon in lower bar of monitor), and close all windows left open (except Labmenu, and the main program being in use).

Don’t forget to jump back to a dialog box if help is consulted, and the box is left open.

Remarks:

To get a reliable value of the wave period, the time step in the measurements must be small enough (at maximum 0.05 times the wave period).
In a later version of this program, the moving average over the length of a period, using the time step of the data file as step, is substracted before calculation of the reflection. This version is called Refreg19, and is stored on the cd in subfolder Lab_process\Extra_s.

· 1. Input parameters

Click on File, Open, Data file,
Search the measurements file to be processed,
Click on Open.

Parameters to be typed via dialog boxes :

read all samples? choose yes (no only in case of files of 10 Mb or more),

number of header lines 0 in case of a Dasylab file (prog. determines length of header); otherwise: the number of header lines,

At this point an error message may occur, see point 7.

time step in original data file,

number of harmonic components to be printed,

average water depth in m.

For each of the two wave gauges:

column number of the wave gauge series in the file,

position of the wave gauge in m,

scale factor of measurements from volts to meters.

Next, a graph of the read samples will appear, together with a dialog box to choose a suitable time interval. The graph can be helpful to choose this interval.

After making your choise, click on OK; a new graph will be displayed, showing the time interval chosen; close this graph to start the calculation.

Next an output file must be chosen. Choose a file, and click on OK, or click on Cancel to stop.

The output file will be opened after the calculation, and the user can view the results.

· 2. Subtraction of mean value

In Refreg, the mean value of the time signal is substracted. In Refreg19 (folder Lab_process\Extra_s), a moving average over the length of the basic period is substracted. The step used in the moving average is the time step in the data file.

· 3. Determination of the length of the signal

In the input part of the program, the user specified the part of the data that can be used in the calculation. From the begin and endpoint of this part, two zero crossings with equal sign are searched. This is done for the first wave gauge. The length of the series between the two zero crossings enables the program to find the correct base period from the FFT.

· 4. Checking the time step

Generally, a period does not contain an integer number of samples. To reduce the error caused by this fact, the Matlab function Resample is used if there less than 20 points per period are present. Note that the phase may be inaccurate in that case.

· 5. Calculating the first harmonic using the FFT of Matlab

The frequency having the maximum C-coefficient will be taken, where C = √(A² + B²).

To coeffficients are determined by the Matlab function fft, and multiplied by 2/N as a scaling factor, where N is the number of points used in the FFT (see "using Matlab, v. 6, p. 6-32).

· 6. Calculating the amplitude and the phase of the incoming and the reflected wave

Equation (6) has the form La = b, where a and b are vectors and L is a matrix. This equation system can be solved in Matlab directly.

The wave number k is determined by the dispersion relation in case of free gravitation surface waves is used:

where ω is the angular frequency and g the gravitational acceleration constant. The calculation is performed by the Matlab-function Disper, written by Gert Klopman.

· 7. How to handle if a file cannot be read

The next errors may occur:

If the number of header lines specified in case of a non-Dasylab file is less than the actual number, the read function used in the scripts gives a message. However, in some cases Matlab gives the next error:

Number of colums in line -- of ASCII file ... must be the same as the previous line.

In such a case, remove all spaces, signs, or numbers from the first position of each header line.

If the number of header lines specified is more than the actual number, no message will occur, but less data lines will be read.

If a Dasylab file gives an error like above, check if the right file has been chosen. If yes, copy this file, delete the first line, and read the copy specifying the remaining number of header lines.

The maximum number of samples to be read is exceeded. This number is set to 500000/nchan,

read all samples?

Specify No in this case.

In a Dasylab file, it is preferable to store the time in decimal format. The time format with a ':' in Dasylab is not read as a time in Matlab.

Program used to test Refreg

To test Refreg, a program ‘testsig’ has been written. In this program, an incoming and a reflected wave with one frequency are generated. The user must specify the following parameters:

name of the output file,

number of wave periods to be generated,

time step,

wave period,

water depth,

distance between water gauges,

reflection coefficient,

phase of the incoming wave,

phase of the reflected wave,

amplitude of normally distributed random noise.

The program has been tested with an amplitude of 1 for the incoming wave (To get an answer in meters, in the program Refreg a scaling factor may be used).

The random noise has been used only to investigate the influence of noise on the frequency of the first harmonic found. If the time step is small, this kind of noise can influence the frequency found, since the zero crossings with equal sign can be shifted by half a period if the frequency of the noise has the order of half the sampling frequency. This kind of noise will normally not be found in waves. It is better to use another kind of disturbance if testing the program. It may be useful to add higher harmonic components to the test signals in the future.

Things to be investigated

1. The influence of the time step, the duration of the measurements, and the wave period on the results.

2. A method to detect the wave period from the measurements, e.g. a regression method.

3. The influence of perturbations, e.g. due to the equipment or caused by undesired objects in the measurement system.

4. A sensitivity analysis of the results. If some noise of a period of about two times the time step is present, an error may be found in the first harmonic component found from the zero crossings. This kind of noise is not expected in surface waves. On the other hand a small shift of the wave gauges can change the results.

A sensitivity analysis can be done by cross correlation methods.

Additional program points in the future

1. Te water velocity as a parameter (now: velocity 0).

2. Adding more (higher) harmonics to the calculation, and splitting into bounded and free components. To be able to find the reflection of more harmonics, more than two wave gauges are needed.

Refreg and Spect, remarks about spectra

June, 2002

Determination of the amplitude spectrum of a regular wave from measurements

1. Sample time: at least 20 points per period are needed.

2. Duration: at least 20 periods are needed, starting from a point where the signal is stationary.

3. After the data has been stored, start the computation. First of all, subtract the mean value from the measurements to avoid a high value in the spectrum at frequency zero.

4. The harmonic components of the signal can be found from a Fast Fourier Transform (FFT)-spectrum. The dominating component in the spectrum is the basic wave period. The frequencies in the spectrum are at discrete values k * 1/T for k = 0,1,…,N/2, where T is the record time and N the total number of samples. As a concequence, having the basic wave period as a spectral component, the part of the signal used in the calculation of the FFT must be an integer number of periods.

This can be obtained as follows:

Method 1:

Determine a zero-crossing at the begin and at the end of the stationary part of the data. The two zero-crossings must have the same sign:

Use the points nearest to the begin and the end of the stationary part of the measured data as begin and end point of the analyses.

Method 2:

If the basic wave period is known beforehand, the multiple of the number of samples per period nearest to the total number of samples can be calculated:

The number of samples per period is

N_p = floor (T_w/ ∆t),

where

floor( ): a function returning the nearest integer <= the argument,

T_w: the wave period,

∆t: the time step in the measurements,

N_p: the rounded-off number of periods.

The number of periods in the used part of the signal is:

N_c = floor (N_s / N_p),

where

N_s: the total number of samples in the analysed part of the signal (= stationary part),

N_c: the rounded-off number of periods in the analysed part of the signal.

The number of points to be used in the FFT is

N_fft = N_c * N_p (1)

If N_fft> N_s, use (N_c-1) * N_p points.

The methods 1 and 2 described above give an integer number of periods. Generally, method 1 gives a better result than method 2.

In some applications, the number of points in the FFT must be a power of two. In that case other methods are needed to find the basic wave period.

5. Calculate the FFT, using the number of points calculated as described above.

6. Determine the maximum absolute value in the FFT-spectrum.

The FFT of a time signal x(n∆t), 0 ≤ n ≤ N-1

can be written as

where i = √-1 and N must be replaced by N_fft as found from eqation (1).

X(k) can be written as

X(k) = A(k) + i B(k),

and the absolute value of X(k) is

Find the maximum of |X(k)|.

If k₁ is the index of the first maximum, the basic harmonic is at frequency k₁ / N_fft∆t.

The second harmonic component can be found at k₂ ≈ 2k₁. Since the round-off error in the basic frequency found is about 0.5 times the frequency step, search between the indices 2k₁-1 and 2k₁+1 for a lokal maximum of |X(k)|.

The same method can be used to find more harmonic components (at 3k₁, 4k₁, …). Because of the round-off error mentioned above, search in the interval [round(i*(k₁-1)), round(i*(k₂-1))].

contact person:

Hetty Klaasman, k. 0.07, tel. (27) 85974.

References

Bakkenes, Henk Jan (2002), Observation and separation of bound and free low-frequency waves in the nearshore zone, MSc-Thesis, Delft University of Technology, Faculty of Civil Engineering and Geosciences, Fluidmechanics section, June 2002.

Battjes,J.A., Windgolven. Collegehandleiding b78, uitgave 1984, 5e herdruk 1992 (Dutch version).

Goda, Y. en Y. Suzuki (1976). "Estimation of Incident and Reflected Waves in Random Wave

Experiments". Proc. 15^th Int Conf. Coastal Engrg., Honolulu, ASCE, New York, N.Y., 1,

pp. 828-845.

Goda, Y. (1985). Random seas and design of maritime structures. Univ. Tokyo Press, Tokyo, Japan.

WL Delft Hydraulics (2000), Delft-Auke Process Description, July 2000.

Zelt, J.A. en J.E. Skjelbreia (1992), Estimating Incident and Reflected Wave Fields Using an Arbitrary Number of Wave Gauges, Proc. 23^rd Int. Conf. Coastal Eng., ASCE, pp. 777-789.