Aalborg Universitet

Analysis of Waves
Technical documentation for WaveLab 3
Frigaard, Peter; Andersen, Thomas Lykke

Publication date:
2014

Document Version
Publisher's PDF, also known as Version of record

Link to publication from Aalborg University

Citation for published version (APA):

Frigaard, P., & Andersen, T. L. (2014). Analysis of Waves: Technical documentation for WaveLab 3. Department
of Civil Engineering, Aalborg University. DCE Lecture notes No. 33

General rights
Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners
and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights.

- Users may download and print one copy of any publication from the public portal for the purpose of private study or research.
- You may not further distribute the material or use it for any profit-making activity or commercial gain
- You may freely distribute the URL identifying the publication in the public portal -
Take down policy
If you believe that this document breaches copyright please contact us at [email protected] providing details, and we will remove access to
the work immediately and investigate your claim.

Downloaded from vbn.aau.dk on: septembre 02, 2025

 Analysis of Waves
Technical documentation for WaveLab 3

Peter Frigaard
Thomas Lykke Andersen

ISSN 1901-7286
DCE Lecture Notes No. 33 Department of Civil Engineering
 Aalborg University
Department of Civil Engineering
Water and Soil

DCE Lecture Notes No. 33

Analysis of Waves
Technical documentation for WaveLab 3

by

Peter Frigaard
Thomas Lykke Andersen

October 2014


c Aalborg University
Scientific Publications at the Department of Civil Engineering
Technical Reports are published for timely dissemination of research re-
sults and scientific work carried out at the Department of Civil Engineering
(DCE) at Aalborg University. This medium allows publication of more de-
tailed explanations and results than typically allowed in scientific journals.
Technical Memoranda are produced to enable the preliminary dissemi-
nation of scientific work by the personnel of the DCE where such release is
deemed to be appropriate. Documents of this kind may be incomplete or
temporary versions of papers—or part of continuing work. This should be
kept in mind when references are given to publications of this kind.
Contract Reports are produced to report scientific work carried out un-
der contract. Publications of this kind contain confidential matter and are
reserved for the sponsors and the DCE. Therefore, Contract Reports are
generally not available for public circulation.
Lecture Notes contain material produced by the lecturers at the DCE for
educational purposes. This may be scientific notes, lecture books, example
problems or manuals for laboratory work, or computer programs developed
at the DCE.
Theses are monograms or collections of papers published to report the sci-
entific work carried out at the DCE to obtain a degree as either PhD or
Doctor of Technology. The thesis is publicly available after the defence of
the degree.
Latest News is published to enable rapid communication of information
about scientific work carried out at the DCE. This includes the status of
research projects, developments in the laboratories, information about col-
laborative work and recent research results.

Published 2014 by
Aalborg University
Department of Civil Engineering
Sohngaardsholmsvej 57,
DK-9000 Aalborg, Denmark

Printed in Denmark at Aalborg University

ISSN 1901-7286 DCE Lecture Notes No. 33

Contents

1 Introduction 5

2 Time-Domain Analysis of Waves 7

2.1 Definition of the individual wave : Zero-down crossing . . . . . 7

2.2 Characteristic wave heights and periods . . . . . . . . . . . . . 8

2.3 Distribution of individual wave heights . . . . . . . . . . . . . 10

2.4 Maximum wave height Hmax . . . . . . . . . . . . . . . . . . . 14

2.5 Distribution of wave periods . . . . . . . . . . . . . . . . . . . 18

2.6 Skewness, kurtosis and atiltness . . . . . . . . . . . . . . . . . 18

3 Frequency-Domain Analysis 21

3.1 Basic concepts of linear waves . . . . . . . . . . . . . . . . . . 21

3.2 Example of variance spectrum . . . . . . . . . . . . . . . . . . 23

3.3 Fourier series . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.4 Discrete signal analysis . . . . . . . . . . . . . . . . . . . . . . 29

3.5 Characteristic wave height and period . . . . . . . . . . . . . . 34

4 Reflection Analysis of Long-Crested Waves 37

4.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

4.2 Wave Reflection . . . . . . . . . . . . . . . . . . . . . . . . . . 38

4.3 Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

1
 4.4 Goda & Suzuki’s method . . . . . . . . . . . . . . . . . . . . . 40

4.5 Mansard & Funke’s method . . . . . . . . . . . . . . . . . . . 44

4.6 Zelt & Skjelbreia’s method . . . . . . . . . . . . . . . . . . . . 49

4.7 Robustness to possible errors . . . . . . . . . . . . . . . . . . . 52

4.7.1 Random noise on signals . . . . . . . . . . . . . . . . . 52

4.7.2 Non-linear effects in the waves . . . . . . . . . . . . . . 54

4.7.3 Numerical problems . . . . . . . . . . . . . . . . . . . . 55

5 Seperation of Incident and Reflected Long-Crested Waves Using Digital Filters

5.1 Principle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.2 Design of Filters . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5.2.1 Numerical Test Examples . . . . . . . . . . . . . . . . 64

5.2.2 Physical Model Tests . . . . . . . . . . . . . . . . . . . 67

5.3 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6 Reflection Analysis of Oblique Long-Crested Waves 71

7 Methods for estimation of directional spectra 75

7.1 The Maximum Likelihood Method . . . . . . . . . . . . . . . . 79

7.2 MLM utilising standard spectra . . . . . . . . . . . . . . . . . 83

7.3 The Bayesian Directional spectrum estimation Method . . . . 89

7.4 Tests of the presented methods . . . . . . . . . . . . . . . . . 93

7.4.1 MLM utilising standard spectra on numerical data . . 93

7.4.2 BDM on numerical data . . . . . . . . . . . . . . . . . 97

7.4.3 BDM on Model Test Data . . . . . . . . . . . . . . . . 98

7.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

8 Wave Groups 101

 8.1 Description of Wave Groups . . . . . . . . . . . . . . . . . . . 102

8.2 Hilbert Transform Technique . . . . . . . . . . . . . . . . . . . 105

8.2.1 Computation of half the square envelope . . . . . . . . 107

8.3 Groupiness Factor . . . . . . . . . . . . . . . . . . . . . . . . . 110

8.4 Conclusions and Further Use . . . . . . . . . . . . . . . . . . . 112

9 References 115
Chapter 1

Introduction

The present book describes the most important aspects of wave analysis tech-
niques applied to physical model tests. Moreover, the book serves as technical
documentation for the wave analysis software WaveLab 3, cf. Aalborg Uni-
versity (2012). In that respect it should be mentioned that supplementary to
the present technical documentation exists also the online help document de-
scribing the WaveLab software in detail including all the inputs and output
fields. In addition to the two main authors also Tue Hald, Jacob Helm-
Ptersen and Morten Møller Jakobsen have contributed to the note. Their
input is highly acknowledged.

The outline of the book is as follows:

• Chapter 2 and 3 describes analysis of waves in time and frequency

domain.

• Chapter 4 and 5 describes the separation of incident and reflected waves

for the two-dimensional case.

• Chapter 6 describes the estimation of the directional spectra which

also includes the separation of incident and reflected waves for multi-
directional waves.

• Chapter 7 describes the analysis of wave groups and calculation of the

groupiness factor.

5
Chapter 2

Time-Domain Analysis of
Waves

2.1 Definition of the individual wave : Zero-

down crossing

The individual wave is defined by two successive zero-down crossings, cf.

2.1. For many years it was common to use zero-up crossings to define a
wave, but due to the asymmetry of natural waves, the greatest wave forces
are often experienced when the wave front hits a structure. That’s one of
the main reasons why IAHR (1986) recommended that the height of a wave
is defined as the height from a trough to the following crest in a time series.
Fig. 2.2 is an example of surface elevation recordings. The application of

Figure 2.1: Individual waves defined by zero-down crossing.

zero-downcrossing gives 15 individual waves (N=15). In Table 4.1 the data

are arranged according to the descending order of wave height.

7
 Figure 2.2: Application of zero-down crossing.

Table 4.1. Ranked individual wave heights and corresponding periods in Fig. 2.2.

rank i 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

H (m) 5.5 4.8 4.2 3.9 3.8 3.4 2.9 2.8 2.7 2.3 2.2 1.9 1.8 1.1 0.23

T (s) 12.5 13.0 12.0 11.2 15.2 8.5 11.9 11.0 9.3 10.1 7.2 5.6 6.3 4.0 0.9

wave no.
in 2.2 7 12 15 3 5 4 2 11 6 1 10 8 13 14 9

2.2 Characteristic wave heights and periods

Usually a surface elevation recording, exemplified in Fig. 2.2, contains more

than several hundred individual waves.

Both wave height and wave period can be considered as random variables,
which have certain probability distributions.

Before these distributions are discussed, some definitions of characteristic

waves will be given.

Mean wave: H, T

H and T are the mean values of the heights and periods , respectively, of all

8
individual waves. Table 4.1 yields

1  1 
15 15
H = Hi = 2.9 m T = Ti = 9.25 s
15 i=1 15 i=1

Root-mean-square wave: Hrms

This wave has a height defined as



 1 N
Hrms =  H2
N i=1 i

From Table 4.1 is found



 1  15
Hrms =  H 2 = 3.20 m
15 i=1 i

Significant wave: Hs , Ts or H1/3 , TH1/3

The significant wave height is the average of the wave heights of the one-third
highest waves. The significant wave period is the average of the wave periods
of the one-third highest waves. From Table 4.1 one finds

1  1 
5 5
Hs = Hi = 4.44 m Ts = Tj = 12.8 s i, j are the rank no.
5 i=1 5 j=1

The significant wave is very often used as the design wave. The reason might
be that in old days structures were designed an a basis of visually observed
waves. Experiences show that often the wave height and period reported
by visual observation correspond approximately to the measured significant
wave. Therefore the choice of significant wave as design wave can make use
of the existing engineering experience.

9
Maximum wave: Hmax , THmax

This is the wave, which has the maximum wave height. In Table 4.1,

Hmax = 5.5 m THmax = 12.5 s

Note, however, that Hmax is a random variable which depends on the number
of individual waves in the time series.

The maximum wave from a long time series corresponding to a storm with a
return period of e.g. 100 years is often chosen as the design wave for struc-
tures which are very important and very sensitive to wave loads.

Highest one-tenth wave: H1/10 , TH1/10

H1/10 is the average of the wave heights of the one-tenth highest waves. TH1/10
is the average of the wave periods of the one-tenth highest wave.

Wave height with exceedence probability of α%: Hα%

It is often practical to denote a wave height according to the probability of

exceedence. Examples are H0.1% , H1% , H2% etc. In many situations H0.1% in
the 100 year storm is used as the design wave.

2.3 Distribution of individual wave heights

Histogram of wave heights

Instead of showing all individual wave heights, it is easier to use the wave
height histogram which gives information about the number of waves in var-
ious wave height intervals. Fig. 2.3 is the histogram of wave heights corre-
sponding to Table 4.1.

Non-dimensionalized histogram

In order to compare the distribution of wave heights at different locations,

the histogram of wave heights is non-dimensionalized, cf. Fig. 2.4.

When Δ(H/H) approaches zero, the probability density becomes a smooth

curve. Experience and theory have shown that this curve is very close to

10
 Figure 2.3: Wave height histogram.

Figure 2.4: Non-dimensionalized wave height histogram.

the Rayleigh distribution in case of deep water waves. In other words, the
individual wave heights follow the Rayleigh distribution.

Rayleigh distribution

The Rayleigh probability density function is defined as

π  π  H
f (x) = x exp − x2 where x= (2.1)
2 4 H
11
The Rayleigh distribution function is
 π 
F (x) = Prob{X < x} = 1 − exp − x2 (2.2)
4

Relation between characteristic wave heights

Figure 2.5: Relation between Hs and H.

If we adopt the Rayleigh distribution as an approximation to the distribution

of individual wave heights, then the characteristic wave heights H1/10 , H1/3 ,
Hrms and Hα% can be expressed by H through the manipulation of the
Rayleigh probability density function.

H1/10 = 2.03 H
H1/3 = 1.60 H

Hrms = 1.13 H (2.3)

H2% = 2.23 H
H0.1% = 2.97 H

Fig. 2.5 illustrates how to obtain the relation between Hs and H.

The Rayleigh distribution function given by Hs instead of H reads


2
H
F (H) = 1 − exp −2 (2.4)
Hs

12
Individual wave height distribution in shallow water

Only in relatively deep water, the Rayleigh distribution is a good approxima-

tion to the distribution of individual wave heights. When wave breaking takes
place due to limited water depth, the individual wave height distribution will
differ from the Rayleigh distribution.

Stive (1986), proposed the following empirical correction to the Rayleigh

distribution based on model tests and also roughly checked against some
prototype data.

1
2
− 13
ln100 Hm0
H1% = Hm0 1+ (2.5)
2 h
1
2
− 12
ln1000 Hm0
H0.1% = Hm0 1+ (2.6)
2 h

where h is the water depth, H1% means the 1% exceedence value of the wave
height determined by zero-down-crossing analysis, whereas the significant
wave height Hm0 is determined from the surface elevataion spectrum, cf.
chapter 3. The correction formulae are very useful for checking the wave
height distribution in small scale physical model tests, cf. Fig. 2.6.

Figure 2.6: Comparison of the expression by Stive, 1986, for shallow water
wave height distribution with model test results. Aalborg University Hy-
draulics Laboratory 1990.

13
Klopmann et al. (1989) proposed also a semi-empirical expression for the
individual wave height distribution. Worth mentioning is also the very much
used method of Battjes and Groenendijk (2000) which was calibrated against
a lot of physical model test data, but the numerical procedure is somewhat
more complicated.

Chapter 3 gives a more detailed discussion on the validity of the Rayleigh

distribution, based on energy spectrum width parameter.

2.4 Maximum wave height Hmax

As mentioned above Hmax is a random variable that depends on the number

of waves in the timeseries. Below some facts about the distribution of Hmax
are given.

Distribution of Hmax

The distribution function of X = H/H is the Rayleigh distribution

 π 
FX (x) = Prob{X < x} = 1 − exp − x2 (2.7)
4

If there are N individual waves in a storm1 , the distribution function of

Xmax = Hmax /H is

FXmax (x) = Prob{Xmax < x} = ( FX (x) )N

  π  N
= 1 − exp − x2 (2.8)
4

Note that FXmax (x) can be interpreted as the probability of the non-occurrence
of the event ( X > x ) in any of N independent trials. The probability density
function of Xmax is

dFXmax
fXmax (x) =
dx
π  π    π  N −1
= N x exp − x2 1 − exp − x2 (2.9)
2 4 4
1
A storm usually lasts some days. The significant wave height is varying during a
storm. However we are more interested in the maximum significant wave height in a short
period of time. In practice, N is often assumed to be 1000.

14
The density function of X and the density function of Xmax are sketched in
Fig. 2.7.

Figure 2.7: Probability density function of X and Xmax .

Mean, median and mode of Hmax

Mean, median and mode are often used as the characteristic values of a
random variable. Their definitions are given in Fig. 2.8.

 +∞
Mean xmean = E[X] = xfX (x)dx
−∞


Median xmedian = x
FX (x)=0.5


Mode xmode = x
fX (x)=max

Figure 2.8: Probability density function of X and Xmax .

By putting eqs (2.8) and (2.9) into the definitions, we obtain


ln N 0.577
(Hmax )mean ≈ + √ Hs (2.10)
2 8 ln N

ln N
(Hmax )mode ≈ Hs (2.11)
2

15
Furthermore, (Hmax )μ , defined as the maximum wave height with exceedence
probability of μ (cf. Fig. 2.9), is

 ⎛ ⎞

1
 N
(Hmax )μ ≈  ln ⎝   ⎠ Hs (2.12)
2 ln 1
1−μ

Obviously (Hmax )median = (Hmax )0.5 .

We get for the commonly used N = 1000 the following maximum wave
heights:

(Hmax )mode ≈ 1.86Hs = 2.97H (2.13)

(Hmax )mean ≈ 1.94Hs = 3.10H (2.14)
(Hmax )10% ≈ 2.14Hs = 3.42H (2.15)

Figure 2.9: Definition of (Hmax )μ .

16
Monte-Carlo simulation of Hmax distribution

The distribution of Hmax can also be studied by the Monte-Carlo simulation.

Individual wave heights follow the Rayleigh distribution

2
H
F (H) = 1 − exp −2 (2.16)
Hs

The storm duration corresponds to N individual waves.

1) Generate randomly a data between 0 and 1. Let the non-exceedence

probability F (H) equal to that data. One individual wave height
H is obtained by (cf. Fig. 2.10)

− ln(1 − F (H))
H = F −1 ( F (H) ) = Hs (2.17)
2

2) Repeat step 1) N times. Thus we obtain a sample belonging to the

distribution of eq (2.16) and the sample size is N .

3) Pick up Hmax from the sample.

4) Repeat steps 2) and 3), say, 10,000 times. Thus we get 10,000
values of Hmax .

5) Draw the probability density of Hmax .

Figure 2.10: Simulated wave height from the Rayleigh distribution.

17
2.5 Distribution of wave periods

It is summarized as

• When we talk about the distribution of wave periods, we often mean

the joint distribution of significant wave height and significant wave
period. Until now there is no general theoretical expression for the
joint distribution, even though there are some so-called scatter dia-
grams based on wave recording. Such a diagram is valid only for the
measurement location. The relation between Hs and Ts is often sim-
plified as Ts = αHsβ , e.g. in Canadian Atlantic waters α = 4.43 and
β = 0.5 (Neu 1982).

• The distribution of wave periods is narrower than that of wave height.

• The empirical relation Tmax ≈ T1/10 ≈ T1/3 ≈ 1.2 T (Goda 2010).

2.6 Skewness, kurtosis and atiltness

Non-linear waves do not have a gaussian distributed surface as the waves are
assymetric with higher and narrower crests and wider and less deep troughs.
This distortion of the distribution of the surface can be described by the
skewness (β1 ) and kurtosis (β2 ) defined as:

 1 
N
1
β1 = 3
· (ηi − η)3 (2.18)
ηrms N i=1

1 
N
1
β2 = 4
· (ηi − η)4 (2.19)
ηrms N i=1

For a Gaussian surface the skewness is zero. Positive skewness indicate that
the median value is below the mean value and a heavy positive tail which
will be the case for non-linear waves. Huang and Long (1980) found that in
deep water the skewness is proportional to the wave skewness. The kurtosis
describes the peakedness of the statistical distribution. A gaussian process
has a kurtosis of 3 and higher values indicate high narrow peak and heavier
tails than given by the gaussian distribution.

Goda (1986) introduced the overall atiltness paramter (β3 ) which describes

18
the assymetry in the direction of propagation.

N
−1
1
N −1
(η̇i − η)3
i=1
β3 = 3/2
(2.20)
N
−1
1
N −1
(η̇i − η)2
i=1

Waves in the near shore zone as the waves tend to have a steep front slope
and a gentle rear slope. Such waves would have a positive atiltness and the
atiltness paramter can exceed unity within the surf zone.

19
Chapter 3

Frequency-Domain Analysis

The concept of a spectrum can be attributed to Newton, who discovered

that sunlight can be decomposed into a spectrum of colors from red to vio-
let, based on the principle that white light consists of numerous components
of light of various colors (wave length or wave frequency).

Energy spectrum means the energy distribution over frequency. Spectral

analysis is a technique of decomposing a complex physical phenomenon into
individual components with respect to frequency.

Spectral analysis of irregular waves is very important for the design of struc-
tures. For example, in the oil-drilling platform design where wave forces
plays an important role, it is of importance to design the structure in such a
way that the natural frequency of the structure is rather far away from the
frequency band where the main part of wave energy concentrates. In this
way resonance phenomenon and the corresponding dynamic amplification of
force and deformation can be avoided.

3.1 Basic concepts of linear waves

Surface elevation
The surface elevation of a linear wave is:

H
η(x, t) = cos(ωt − kx + δ) = a cos(ωt − kx + δ) (3.1)
2
21
 where H wave height
a amplitude, a = H/2
ω angular frequency, ω = 2π/T
T wave period.
k wave number, k = 2π/L
L wave length
δ initial phase

We can also define the observation location to x = 0 and obtain

η(t) = a cos(ωt + δ) (3.2)

The relation between wave period and wave length (dispersion relationship)
is

g T2 2πh
L = tanh (3.3)
2π L

where h is water depth.

Wave energy

The average wave energy per unit area is:

1 1
E = ρ g H 2 = ρ g a2 (Joule/m2 in SI unit) (3.4)
8 2

Variance of surface elevation of a linear wave

The variance of the surface elevation of a sinus wave is:

 2 
Var[η(t)] = E η(t) − η(t) (E: Expectation)
 
=E ( η(t) )2

1
T
= T 0
η 2 (t) dt (T: wave period)

1
= 2
a2

22
Superposition of linear waves

Since the governing equation (Laplace equation) and boundary conditions are
linear in small amplitude wave theory, it is known from mathematics that
small amplitude waves are superposable. This means that the superposition
of a number of linear waves with different wave height and wave period will
be:

superposition wave 1 wave 2 ··· wave N

velocity potential ϕ = ϕ1 + ϕ2 + ··· + ϕN

surface elevation η = η1 + η2 + ··· + ηN
particle velocity u = u1 + u2 + ··· + uN
dynamic pressure p = p1 + p2 + ··· + pN

3.2 Example of variance spectrum

First we will make use of an example to demonstrate what a variance spec-

trum is.

Surface elevation of irregular wave

Fig. 3.1 gives an example of an irregular wave surface elevation which is
constructed by adding 4 linear waves (component waves) of different wave
height and wave period. The superposed wave surface elevation is


4 
4
η(t) = ηi (t) = ai cos(ωi t + δi ) (3.5)
i=1 i=1

23
Figure 3.1: Simulation of irregular waves by superposition of linear waves.

Variance diagram
Instead of Fig. 3.1, we can use a variance diagram, shown in Fig. 3.2, to
describe the irregular wave.

Figure 3.2: Variance diagram.

In comparison with Fig. 3.1, the variance diagram keeps the information on
amplitude (ai ) and frequency (fi , hence Ti and Li ) of each component, while
the information on initial phase (δi ) is lost. This information loss does not

24
matter because the surface elevation of irregular wave is a random process.
We can simply assign a random initial phase to each component.

Variance spectral density Sη (f )

The variance diagram can be converted to variance spectrum, The spectral
density is defined as

1 2
2
a
Sη (f ) = (m2 s) (3.6)
Δf

where Δf is the frequency band width1 , cf. Fig. 3.3.

Figure 3.3: Stepped Variance spectrum.

In reality an irregular wave is composed of infinite number of linear waves

with different frequency. Fig. 3.4 gives an example of stepped variance
spectrum.

Figure 3.4: Continuous variance spectrum (wave energy spectrum).

1
we will see later that Δf depends on signal recording duration. In the figure it is
assumed that Δf = 0.01Hz

25
When Δf approaches zero, the variance spectrum becomes a continuous
curve. A variance spectrum is also called energy spectrum. But strictly
speaking, the energy spectral density should be defined as

1
ρ g a2
S(f ) = 2
( Jm2s ) (3.7)
Δf

Construction of time series from variance spectrum

We can also construct time series of surface elevation from variance spectrum.
In Fig. 3.4 the known variance spectral density Sη (f ) is divided into N parts
by the frequency band width Δf . This means that the irregular wave is
composed of N linear waves


N 
N
η(t) = ηi (t) = ai cos(ωi t + δi ) (3.8)
i=1 i=1

The variance of each linear wave is

1 2
Sη (fi ) Δf = a i = 1, 2, · · · , N (3.9)
2 i

Therefore the amplitude is


ai = 2 Sη (fi ) Δf i = 1, 2, · · · , N (3.10)

The angular frequency is

2π
ωi = = 2πfi i = 1, 2, · · · , N (3.11)
Ti

The initial phase δi is assigned a random number between 0 and 2π. Hence
by use of Eq. (3.8) we can draw the time-series of the surface elevation of the
irregular wave which has the variance spectrum as shown in Fig. 3.4.

26
3.3 Fourier series

Conversion of irregular surface elevation into variance spectrum is not as

simple as the above example, where the linear components of the irregular
wave are pre-defined (cf. Fig. 3.1). We need to decompose the irregular wave
into its linear components. First let’s see how it can be done with a known
continuous function x(t).

Fourier series is used to represent any arbitrary function2 .

Figure 3.5: Arbitrary periodic function of time.


∞
2πi 2πi
x(t) = a0 + ai cos t + bi sin t
i=1
T0 T0


∞
= (ai cos ωi t + bi sin ωi t) (3.12)
i=0

2
Not all mathematicians agree that an arbitrary function can be represented by a
Fourier series. However, all agree that if x(t) is a periodic function of time t, with period
T0 then x(t) can be expressed as a Fourier series. In our case x(t) is the surface elevation
of irregular wave, which is a random process. if T0 is large enough, we can assume that
x(t) is a periodic function with period T0 .

27
where ai and bi are Fourier coefficients given by

1
 T0
a0 = T0 0
x(t) dt and b0 = 0
2
 T0 ⎫
ai = T0 0
x(t) cos ωi t dt ⎬
 T0 i = 1, 2, 3, · · · , ∞ (3.13)
2 ⎭
bi = T0 0
x(t) sin ωi t dt

Physical interpretation

Now we say that the continuous function x(t) is the surface elevation of
irregular wave. η(t) can be expanded as a Fourier series.


∞
η(t) = (ai cos ωi t + bi sin ωi t)
i=0

∞ 
 
 
= ci sin δi cos ωi t + ci cos δi sin ωi t
i=0


∞
 
= ci (sin δi cos ωi t + cos δi sin ωi t)
i=0


∞

= ci sin(ωi t + δi )
i=0


∞
= ci cos(ωi t + δi ) (3.14)
i=0


where δi = δi − π
2
and sin(x + y) = sin(x) cos(y) + cos(x) sin(y) have been
used.

That is to say, any irregular wave surface elevation, expressed as a continues

28
function, is composed of infinite number of linear waves with

 
amplitude ci = a2i + b2i
i = 0, 1, · · · , ∞ (3.15)
2π T0
period Ti = ωi
= i

{ai , bi }, i = 0, 1, 2, · · · , ∞, are given in Eq. (3.13).

3.4 Discrete signal analysis

The measurement of surface elevation is carried out digitally. We do not have,

neither necessary, a continuous function of the surface elevation. Instead we
have a series of surface elevation measurement equally spaced in time, cf.
Fig. 3.6.

Figure 3.6: Sampling of surface elevation at regular intervals.

If the sampling frequency is fs , then the time interval between two succeeding
points is Δt = 1/fs . Corresponding to the sample duration T0 the total
number of sample is N = T0 /Δt. Thus we obtain a discrete time series of
surface elevation

η0 , η1 , ··· , ηN −1

The Fourier coefficients

(a0 , b0 ), (a1 , b1 ), ··· , (aN −1 , bN −1 )

29
can be obtained by Fast Fourier Transforms (FFT)3 . That is to say, the
irregular wave surface elevation, expressed by digital time series, is composed
of N linear waves


N −1 
N −1 
η(t) = ηi (t) = a2i + b2i cos(ωi t + δi ) (3.16)
i=0 i=0

 ⎫
amplitude a2i + b2i ⎪
⎪
⎪
⎪
angular frequency ωi = 2πi ⎬
T0
i = 0, 1, · · · , N − 1 (3.17)
period Ti = 2π
= Ti0 ⎪
⎪
ωi ⎪
⎪
1 ⎭
frequency fi = Ti
= Ti0

Therefore we obtain the variance spectrum

1
frequency width Δf = fi+1 − fi =
T0
1 1 2
2
(amplitude)2 (a
2 i
+ b2i )
spectral density Sη (fi ) = = (3.18)
Δf Δf

An example of variance spectrum is shown in Fig. 3.7.

Figure 3.7: Variance spectrum.

3
FFT is a computer algorithm for calculating DFT. It offers an enormous reduction in
computer processing time. For details of DFT and FFT, refer to Newland (1975)

30
Nyquist frequency fnyquist

Nyquist frequency fnyquist is the maximum frequency which can be detected

by the Fourier analysis.

Fourier analysis decomposes N digital data into N linear components. The

frequency of each component is

i
fi = i = 0, 1, · · · , N − 1 (3.19)
T0

The Nyquist frequency is

N 1 T0
2 2 Δ fs
fnyquist = f N = = = (3.20)
2 T0 T0 2

where fs sample frequency

Δ time interval between two succeeding sample points, Δ = 1/fs
T0 sample duration
N total number of sample, N = T0 /Δ

The concept of Nyquist frequency means that the Fourier coefficients { ai , bi

}, i = 0, 1, · · · , N − 1, contains two parts, the first half part below the
Nyquist frequency (i = 0, 1, · · · , N/2 − 1) represents true components while
the second half part (i = N/2, N/2 + 1, · · · , N − 1) is the folding components
(aliasing).

Fig. 3.8 gives an example on aliasing after the Fourier analysis of discrete
time series of a linear wave.

The solution to aliasing is simple: let { ai , bi }, i = N/2, N/2 + 1, · · · , N − 1,

equal to zero, cf. Fig. 3.9. That is the reason why fnyquist is also called
cut-off frequency. In doing so we are actually assuming that irregular wave
contains no linear components whose frequency is higher than fnyquist . This
assumption can be assured by choosing sufficiently high sample frequency fs
in combination with an analog low-pass filter, cf. Eq. 3.20.

31
 Figure 3.8: Aliasing after Fourier analysis.

Figure 3.9: Variance spectrum after cut-off (refer to Fig. 3.7).

32
Taper data window

Fourier analysis requires that η(t) is a periodic function with period T0 , it

may be desirable to modify the recorded time series before Fourier analysis,
so that the signal looks like a periodic function. The modification is carried
out with the help of taper data window.

The widely-used cosine taper data window reads

⎧  
⎪
⎪ 1
1 − cos 10πt 0≤t≤ T0
⎪
⎪ 2 T0 10
⎪
⎪
⎨
d(t) = 1 T0
10
≤t≤ 9T0
10 (3.21)
⎪
⎪
⎪
⎪ 9T0
⎪
⎪ 1 10π (t− 10 ) 9T0
≤ t ≤ T0
⎩ 2
1 + cos T0 10

Figure 3.10: Taper data window.

33
3.5 Characteristic wave height and period

The variance spectrum, illustrated in Fig. 3.11, says nothing about how
high the individual waves will be. Now we will see how to estimate the
characteristic wave height and period based on the variance spectrum.

Figure 3.11: Variance spectrum.

n order moment mn

mn is defined as
! ∞
mn = f n Sη (f ) df (3.22)
0

The zero moment is

! ∞
m0 = Sη (f ) df (3.23)
0

which is actually the area under the curve, cf. Fig. 3.11.

Spectrum width parameter and validity of the Rayleigh distribution

From the definition of mn , it can be seen that the higher the order of moment,
the more weight is put on the higher frequency portion of the spectrum. With
the same m0 , a wider spectrum gives larger values of the higher order moment
(n ≥ 2). Longuet-Higgins has defined a broadness parameter:

"
m22
ε4 = 1 − (3.24)
m0 m4

34
An alternative parameter is the narrowness paramter:

m0 m2
ε2 = −1 (3.25)
m21

spectrum width parameter wave height distribution

It has been proven theoretically that: ε=0 narrow spectrum Rayleigh distribution
ε=1 wide spectrum Normal distribution

In reality ε lies in the range of 0.4-0.5. It has been found that Rayleigh dis-
tribution is a very good approximation and furthermore conservative, as the
Rayleigh distribution gives slightly larger wave height for any given proba-
bility level.

Significant wave height Hm0 and peak wave period Tp

When wave height follows the Rayleigh distribution, i.e. ε = 0 , the signifi-
cant wave height Hs 4 can theoretically be expressed as

√
H s ≈ H m0 ≡ 4 m0 (3.26)

In reality where ε = 0.4 − 0.5, a good estimate of significant wave height

from energy spectrum is

√
Hs = 3.7 m0 (3.27)

Peak frequency is defined as (cf. Fig. 3.11)

#
#
fp = f# (3.28)
Sη (f )=max

Wave peak period (Tp = 1/fp ) is approximately equal to significant wave

period defined in time-domain analysis.

4
Hm0 denotes a wave height determined from spectrum while Hs or H1/3 is significant
wave height determined from time-domain analysis. They are equal to each other when
wave height follows the Rayleigh distribution.

35
36
Chapter 4

Reflection Analysis of
Long-Crested Waves

4.1 Introduction

Coastal and offshore engineering cover a wide spectrum in the field of civil
engineering. Some of the major subjects involved are the construction of
harbours and breakwaters.

Design of breakwaters and breakwater lay-outs are based on the need for an
acceptable wave climate in the harbour and in the harbour entrance. How-
ever, the complex behaviour of waves propagating into a more or less closed
basin makes the design of the quaywalls and breakwaters important.

When a physical or numerical model is used for wave climate investigations

it is important that the model will provide reliable results. Naturally this
depends upon whether the boundary conditions in the model are correct.

Establishing the correct boundary conditions is a problem due to scale ef-

fects. However the important properties of a boundary are the absorption,
transmission and reflection of waves.

In order to provide correct reflection characteristics in the models one must

know the true reflection characteristics for the structure under study.

The wave field in front of a structure consist of both incident waves and re-
flected waves. As one cannot first measure the incident waves and then the
reflected waves, methods are required to separate or distinguish between the
incident and reflected waves.

37
4.2 Wave Reflection

Wave reflection occurs when waves are propagating onto breakwaters and
quaywalls.

Separating the incident waves and the reflected waves from a recorded wave
elevation time series has several purposes. Either:

• The reflection characteristic for a structure is wanted, or

• The incident wave characteristic for a given site is wanted

The reflection characteristics for a structure are mainly wanted in sheltered

areas, e.g. inside harbours, where the level of wave disturbance plays an
important role in the design situation. Though, also the reflection character-
istics of breakwaters are wanted, because it is important for the wave climate
in the harbour entrance.

Generally a low reflection is wanted.

In order to apply physical models and numerical models the reflection char-
acteristics must be modeled correct.

The incident waves are wanted where response of a structure (breakwater) is

wanted as function of the incomming waves in front (seaward) of the reflect-
ing structure.

In the following several different methods for separating incident waves and
reflected waves will be presented. After separation of the incident waves and
the reflected waves it is easy to calculate the reflection characteristics for the
structure.

The following methods for separating the waves will be divided into two main
groups:

• 2D methods (treated in this lection)

The waves are assumed to be 2-dimensional (long crested)

• 3D methods (treated in lection 14)

The waves are assumed to be 3-dimensional (short crested)

There is a wide need for 2D methods due to the amount of experiments car-
ried out in wave flumes. In wave flumes reflection and re-reflection influence
upon the validity of measurements if it is not possible to extract the incident

38
waves. Further when operating in the near shore environment the wave-field
will often be approximately two-dimensional as the waves refract, i.e. they
bend towards orthogonals to the shore.

The 2D methods presented herein are all based on surface elevation measure-
ments at a number of positions.

Three of the 2D methods are derived from the same principles but for vari-
ous number of probes. The simplest method requires two probes. The other
methods need three or more probes. The three first methods work in fre-
quency domain and give the incident wave spectrum as well as the reflected
wave spectrum. The last 2D method works in time domain and gives the
incident wave timeserie.

There is a considerable step from 2D methods to 3D methods. This is due to

the circumstance that 3D methods for estimation of random sea states even
without reflection are tedious and not as reliable as 2D methods. In princi-
ple, all methods for estimating directional wavespectra may be considered as
a tool to determine reflection coefficients. However, only a few methods are
capable of handling reflected seas as the reflected waves propagate with the
same frequency as the incident waves, and thereby makes it impossible or at
least difficult to distinguish. One should also recall that the direction of the
waves is unknown in contrast to a 2D case.

As an intermediate case there is oblique long crested waves. This may be

considered as a 3D case having a very narrow peak in the directional spread-
ing function. This justifies the assumption that only one direction is present.
Furthermore in laboratories real oblique waves can be generated. This may
be used to determine the 3D reflection characteristics for e.g. a breakwater
by use of physical models.

4.3 Methods

In the following four methods for separation of incident and reflected waves
in a two-dimensional wave field will be presented.

The first three methods are all working in the frequency domain and have
the same basic principle. They assume the wave elevation to be a sum of
regular waves travelling with different frequency and phase. The first method
by Goda & Suzuki needs measurements of the wave elevation in two distinct
points. Hence by use of Fourier analysis the amplitude of the incident and
reflected waves for a given frequency can be estimated. This procedure does
not account for the noise which probably is contaminated in the measured

39
wave signals, i.e. the measured wave elevations.

The method presented by Mansard & Funke takes this into account, but
also requires the wave elevation to be measured in three distinct points. By
applying Fourier analysis this should result in the same waves except that the
measured noise varies from wave gauge to wave gauge. That is, 10 Fourier
coefficients are needed but only 6 are available. Instead the noise is expressed
in terms of the Fourier coefficients for the incident and reflected waves and
the measured coefficients. The best estimate of the Fourier coefficients is
then found by minimising the noise.

The method by Zelt & Skjelbreia extents the method to apply to an arbitrary
number (though larger than one) of wave gauges. Further a weighting of
each frequency component is introduced. This is to control the possible
influence of singularities which occur for some geometric relations between
the distances between the probes. Grønbech et. al (1996) presented a revised
version of the Zelt & Skjelbreia method to separate cross-modes from the
incident and reflected waves. The method is good to quantify the amount of
cross-mode activity and if neglectable the standard method can be applied.
However, due to the more degrees of freedom more gauges are needed to give
reliable results.

Finally a method by Frigaard & Brorsen is presented in the following chapter.

This method is based on another approach than the previously mentioned
methods. The method has the advantage that it works in the time domain
and can be applied in real time as it is based on digital filters. The other
methods can in post analyses also be used to get the time domain results
by the use of InvFFT on the Fourrier coefficients. However, they cannot be
applied in real time.

Other methods have been published to account for sloping bottoms and non-
linear waves, but they will not be presented herein. The methods by Goda
& Suzuki and Mansard & Funke are the methods being used most often
as they are relatively simple to apply and yield reliable results within most
applications.

4.4 Goda & Suzuki’s method

The method presented by Goda & Suzuki (1976) is the most simple method.
It is based on the assumption that the wave elevation can be considered as a
sum of waves travelling with different frequency, amplitude and phase. Fur-
ther for each wave a reflected wave will travel in the opposite direction. The
method makes use of Fourier analysis and will due to singularities put con-

40
straint to the distance between the waveprobes. The method is very easy to
implement but has a lack of accuracy as no measuring errors are accounted
for.

The surface elevation in a two-dimensional wave field is assumed to be a

summation of a number of waves, say N waves, i.e.


N
η(x, t) = an cos(kn x − ωn t + Φn ) (4.1)
n=1

where η is the surface elevation relative to MWL

x is the position of the wave gauge in a predefined coordinatesystem
t is time
an is the amplitude
kn is the wavenumber
ωn is the angular frequency of the waves
Φn is the phase

If reflection happens 4.1 may be expanded to


N 
N
η(x, t) = aI,n cos(kn x − ωn t + ΦI,n ) + aR,n cos(kn x + ωn t + ΦR,n )(4.2)
n=1 n=1

where indices I and R denotes incident and reflected respectively.

In the following the index n will be omitted for simplicity, that is, only one
frequency is considered. Hence eq. (4.2) will consist of only two terms:

η(x, t) = aI cos(kx − ωt + ΦI ) + aR cos(kx + ωt + ΦR ) (4.3)

Sampling the surface elevation at two positions within a distance of x1,2

measured in the direction of propagation of the waves, as illustrated in fig.

41
 Figure 4.1: Definition sketch showing x1,2 .

(4.1), yields:

η1 = η(x1 , t) = aI cos(kx1 − ωt + ΦI ) + aR cos(kx1 + ωt + ΦR )

(4.4)
η2 = η(x2 , t) = aI cos(kx2 − ωt + ΦI ) + aR cos(kx2 + ωt + ΦR )

Decomposing the trigonometric terms in eq. (4.4) by use of

cos A cos B ± sin A sin B = cos(A ∓ B) (4.5)

leads to

η1 = aI (sin(kx1 + ΦI ) sin(ωt) + cos(kx1 + ΦI ) cos(ωt))

+aR (cos(kx1 + ΦR ) cos(ωt) − sin(kx1 + ΦR ) sin(ωt)) (4.6)
η2 = aI (sin(kx2 + ΦI ) sin(ωt) + cos(kx2 + ΦI ) cos(ωt))
+aR (cos(kx2 + ΦR ) cos(ωt) − sin(kx2 + ΦR ) sin(ωt)) (4.7)

which is rearranged to

η1 = A1 cos(ωt) + B1 sin(ωt)
(4.8)
η2 = A2 cos(ωt) + B2 sin(ωt)

42
where

A1 = aI cos(kx1 + ΦI ) + aR cos(kx1 + ΦR )
B1 = aI sin(kx1 + ΦI ) − aR sin(kx1 + ΦR )
(4.9)
A2 = aI cos(kx2 + ΦI ) + aR cos(kx2 + ΦR )
B2 = aI sin(kx2 + ΦI ) − aR sin(kx2 + ΦR )

In eq. (4.8) the elevation is seen to be a composite signal of a sine and cosine
signal having different time-constant amplitude, i.e. the LHS of eq. (4.9) must
correspond to the Fourier coefficients which can be obtained from Fourier
analysis of the recorded time series.

Thus eq. (4.9) contains four equations with four unknowns, i.e. aI , aR , ΦI
and ΦR . The solution giving the amplitudes is given by Goda & Suzuki
(1976) as:


aI = 1
2| sin(kx1,2 )|
(A2 − A1 cos(kx1,2 ) − B1 sin(kx1,2 ))2

+(B2 + A1 sin(kx1,2 ) − B1 cos(kx1,2 ))2

 (4.10)
aR = 1
2| sin(kx1,2 )|
(A2 − A1 cos(kx1,2 ) + B1 sin(kx1,2 ))2

+(B2 − A1 sin(kx1,2 ) − B1 cos(kx1,2 ))2

where x1,2 = x2 − x1 .

It is seen from eq. (4.10) that singularities will occur for sin(kx1,2 ) = 0.
Hence it should be avoided that
x1,2 n
= where n = 0, 1, 2, ...
L 2
Further Goda & Suzuki (1976) suggest to avoid values in the range ±0.05 xL1,2
at the singularity points. For a wide wave spectrum this will be impossible
for all frequencies, but of course values applying for the peak frequencies
should be weighted highest.

For regular waves the method is quite accurate but for irregular waves the
confidence of the FFT analysis plays a significant role. However in both cases
noise may be the dominant error as it cannot be accounted for.

43
4.5 Mansard & Funke’s method

As a natural extension to the method by Goda & Suzuki, Mansard & Funke
(1980) presented a three-points method. Here an additional probe is taken
into use which makes it possible to add an error to the measurements and
hence minimise it in a least squares sense.

The general equation for a progressive wave field is obtained as in the previous
method, i.e. eq. (4.1):


N
η(x, t) = an cos(kn x − ωn t + Φn ) (4.11)
n=1

The wave elevation given by eq. (4.1) and eq. (4.11) is separated into incident
waves and reflected waves, also as previously, but now a noise function is
added. This leads to:

N 
N
η(x, t) = aI,n cos(kn x − ωn t + ΦI,n ) + aR,n cos(kn x + ωn t + ΦR,n )
n=1 n=1

or


N
η(x, t) = aI,n cos(kn x − ωn t + Φn )
n=1
N
+ aR,n cos(kn (x + 2xR ) + ωn t + Φn + θs ) + Ω(t) (4.12)
n=1

Here in contrary to Goda & Suzuki (1976) the distance from the point of
observation to the reflecting structure is assumed known. However as the
distance may be impossible to determine (e.g. for a slope) a phase θs is
introduced to compensate herefore. This has the consequence that the phase
Φn remains the same for the incident and reflected wave. Ω(t) is the noise
function and expresses all kinds of errors. xR is the distance from the wave
probe to the reflecting structure.

Mansard & Funke allow a phaseshift θs to occur at the reflecting structure.

For three probes placed as indicated on fig 4.2 eq. (4.12) can be applied for

44
 Figure 4.2: Definition sketch.

each probe.

Hence

ηp = η(xp , t)
N
= aI,n cos(kn xp − ωn t + Φn )
n=1
N
+ aR,n cos(kn (xp + 2xR,p ) + ωn t + Φn + θs ) + Ωp (t) (4.13)
n=1

where index p refer to probe number.

Inserting x1,p , which is the distance from the 1st probe to the p’th probe,
eq. (4.13) can be modified to


N
ηp = aI,n cos(kn (x1 + x1,p ) − ωn t + Φn )
n=1
N
+ aR,n cos(kn (x1 + 2xR,1 − x1,p ) + ωn t + Φn + θs ) + Ω(4.14)
p (t)
n=1

In the following until eq. (4.25) only one frequency will be considered, i.e. in-

45
dex n will be omitted. Thus

ηp = aI cos(k(x1 + x1,p ) − ωt + Φ)
+aR cos(k(x1 + 2xR,1 − x1,p ) + ωt + Φ + θs ) + Ωp (t) (4.15)

The following manipulations have the purpose to derive an algebraic solution,

which is based on a minimisation of the noise function.

Fourier transformation of eq. (4.15) yields

Ap + iBp = aI exp (ik(x1 + x1,p ) + iΦ)

+aR exp (ik(x1 + 2xR,1 − x1,p ) + i(Φ + θs ))
+Yp exp(iρp ) (4.16)

where Ap and Bp are the Fourier coefficients.

Let

ZI = aI exp(ikx1 + iΦ) (4.17)

ZR = aR exp(ik(x1 + 2xR,1 ) + i(Φ + θs )) (4.18)
ZN,p = Yp exp(iρp ) (4.19)

where index N refer to noise.

Now, eq. (4.16) can be applied to each probe yielding

A1 + iB1 = ZI + ZR + ZN,1
A2 + iB2 = ZI exp(ikx1,2 ) + ZR exp(−ikx1,2 ) + ZN,2
A3 + iB3 = ZI exp(ikx1,3 ) + ZR exp(−ikx1,3 ) + ZN,3

or in general form

Ap + iBp = ZI exp (ikx1,p ) + ZR exp (−ikx1,p ) + ZN,p (4.20)

One is interested in solving 4.20 with regard to ZI and ZR as they contain

information of the reflection coefficients.

46
Now eq. (4.20) can be rearranged to

ε1 = ZI + ZR − (A1 + iB1 )
ε2 = ZI exp(ikx1,2 ) + ZR exp(−ikx1,2 ) − (A2 + iB2 ) (4.21)
ε3 = ZI exp(ikx1,3 ) + ZR exp(−ikx1,3 ) − (A3 + iB3 )

where

εp = −ZN,p + fe (ZI , ZR )

Minimising the noise function Ωp (t) introduced in eq. (4.12) correspond to

minimising the sum of squares of εp for all p, i.e.


3 
3
2
(εp ) = (ZI exp(ikx1,p ) + ZR exp(−ikx1,p ) − (Ap + iBp ))2 (4.22)
p=1 p=1

should be minimised.

Assuming that the minimum of eq. (4.22) is achieved when both partial
derivatives are zero, i.e.

3 3
∂ p=1 ε2p ∂ p=1 ε2p
= =0 (4.23)
∂ZI ∂ZR

Hence one obtains:

(4.24)

3
0 = 2 p=1 (ZI exp(ikx1,p ) + ZR exp(−ikx1,p ) − (Ap + iBp )) exp(ikx1,p )
3
0 = 2 p=1 (ZI exp(ikx1,p ) + ZR exp(−ikx1,p ) − (Ap + iBp )) exp(−ikx1,p )

47
When written out and again including index n eq. (4.24) leads to

ZI,n (1 + exp(i2kn x1,2 ) + exp(i2kn x1,3 )) + 3ZR,n =

(A1,n + iB1,n ) + (A2,n + iB2,n ) exp(ikn x1,2 ) + (A3,n + iB3,n ) exp(ikn x1,3 )

ZR,n (1 + exp(−i2kn x1,2 ) + exp(−i2kn x1,3 )) + 3ZI,n =

(A1,n + iB1,n ) + (A2,n + iB2,n ) exp(−ikn x1,2 ) + (A3,n + iB3,n ) exp(−ikn x1,3 )
(4.25)

The solution is given by Mansard & Funke (1980) as:

1
ZI,n = Dn ((A1,n + iB1,n )(R1,n + iQ1,n )

+(A2,n + B2,n )(R2,n + iQ2,n ) + (A3,n + iB3,n )(R3,n + iQ3,n ))

1
ZR,n = Dn ((A1,n + iB1,n )(R1,n − iQ1,n )

+(A2,n + iB2,n )(R2,n − iQ2,n ) + (A3,n + iB3,n )(R3,n − iQ3,n ))

(4.26)

where

Dn = 2(sin2 (kn x1,2 ) + sin2 (kn x1,3 ) + sin2 ((kn x1,3 ) − kn x1,2 ))
R1,n = sin2 (kn x1,2 ) + sin2 (kn x1,3 )
Q1,n = sin(kn x1,2 ) cos(kn x1,2 ) + sin(kn x1,3 ) cos(kn x1,3 )
R2,n = sin(kn x1,3 ) sin(kn x1,3 − kn x1,2 )
Q2,n = sin(kn x1,3 ) cos(kn x1,3 − kn x1,2 ) − 2 sin(kn x1,2 )
R3,n = − sin(kn x1,2 ) sin(kn x1,3 − kn x1,2 )
Q3,n = sin(kn x1,2 ) cos(kn x1,3 − kn x1,2 ) − 2 sin(kn x1,3 )

The only unknowns in eq. (4.26) are the Fourier coefficients of the measured
wave elevations. These are obtained by use of FFT analysis of the measure-
ments.

Compared to the previously discussed method this method has the advan-
tage that it minimises the noise contaminated to the elevation measurements.

48
Thus instead of obtaining an exact solution a fitted solution is obtained.

Problems due to singularities, however, occur as well. Singularity will occur

when
Dn = 0
i.e.
sin2 (kn x1,2 ) + sin2 (kn x1,3 ) + sin2 (kn x1,3 − kn x1,2 ) = 0
As all the terms are positive the solution is solution to

sin(kn x1,2 ) = sin(kn x1,3 ) = sin(kn x1,3 − kn x1,2 ) = 0

but if sin(kn x1,2 ) = sin(kn x1,3 ) = 0 then sin(kn x1,3 − kn x1,2 ) = 0 whereas the
solution is reduced to

sin(kn x1,2 ) = sin(kn x1,3 ) = 0

which is obtained for

kn x1,2 = mπ ∧ kn x1,3 = lπ m = 1, 2, 3, ...
l = m + 1, m + 2, ...
kn x1,2 = mπ ∧ kn x1,3 = lπ
x1,2 = m2 Ln ∧ x1,3 = 2l Ln = l
x
m 1,2

as by definition 0 < x1,2 < x1,3 . Mansard & Funke (1980) suggest that

Ln
x1,2 =
10
Ln Ln
< x1,3 <
6 3
Ln
x1,3 =
5
3Ln
x1,3 =
10

4.6 Zelt & Skjelbreia’s method

Zelt and Skjelbreia (1992) introduced a method based on the same principles
as the previously described methods.

This method applies for an arbitrary number p of probes and further intro-
duces a weighting of the measurements from the gauges. The latter facility
takes into account the spacing between each pair of wave gauges. However

49
determining the weighting coefficients is a rather subjective process and re-
quires some additional work.

For p = 2 the method corresponds to the method by Goda & Suzuki. For
p = 3 and a uniform weighting, the method corresponds to the method by
Mansard & Funke.

Once again the wave elevation is considered as a sum of incident and reflected
waves, i.e. for the pth probe at position xp utilising eq. (4.2)

N
ηp = η(xp , t) = n=1 aI,n cos(kn xp − ωn t + ΦI,n ) +
N
n=1 aR,n cos(kn xp + ωn t + ΦR,n )

or by use of eq. (4.14)


N
ηp = aI,n cos(kn (x1 + x1,p ) − ωn t + Φn )
n=1
N
+ aR,n cos(kn (x1 + 2xR,1 − x1,p ) + ωn t + Φn + θs ) + Ωp (t)
n=1

As in the method of Mansard & Funke (1980) this is by use of the Fourier
coefficients rearranged to eq. (4.21) which has the general form

εp,n = ZI,n exp(ikn x1,p ) + ZR,n exp(−ikn x1,p ) − (Ap,n + iBp,n ) (4.27)

That is, if the estimated coefficients are correct there will be no error,
i.e. εp,n = 0. A function is chosen to be a weighted sum of the squares
of εp,n , i.e.


P
En = Wp,n εp,n ε∗p,n (4.28)
p=1

where Wp,n > 0 are weighting coefficients.

It is assumed that the minimum of eq. (4.28) occurs where the partial deriva-

50
tives with regard to ZI,n and ZR,n are zero. Thus


P
Wp,n εp,n exp(±ikn x1,p ) = 0 (4.29)
p=1

Inserting eq. (4.27) into eq. (4.29) gives:

P P
p=1 Wp,n ZI,n exp(i2kn x1,p ) + p=1 Wp,n ZR,n =
P
p=1 Wp,n (Ap,n + iBp,n ) exp(ikn x1,p )
(4.30)
P P
p=1 Wp,n ZI,n + p=1 Wp,n ZR,n exp(−i2kn x1,p ) =
P
p=1 Wp,n (Ap,n + iBp,n ) exp(−ikn x1,p )

This, eq. (4.30) can be solved with regard to ZI,n and ZR,n as a linear system
of equations. Simply isolate aI,n and aR,n in both equations in eq. (4.30) and
substitute into the other equation respectively. Hence the following solution
is obtained

ZI,n = 1
Sn Pp=1 Wp,n (Ap,n + iBp,n ) exp(−ikn x1,p )
Dn
  
− Pp=1 Wp,n (Ap,n + iBp,n ) exp(ikn x1,p ) Pq=1 Wq,n exp(−i2kn x1,q )
 
1 P
ZR,n = Dn Sn p=1 Wp,n (Ap,n + iBp,n ) exp(ikn x1,p )
P  
− p=1 Wp,n (Ap,n + iBp,n ) exp(−ikn x1,p ) Pq=1 Wq,n exp(i2kn x1,q )

(4.31)

where


P 
P
Dn = Sn2 − Wp,n exp(i2kn x1,p ) Wq,n exp(−i2kn x1,q )
p=1 q=1


P
Sn = Wp,n
p=1

It is seen from eq. (4.31) that Dn = 0 will lead to singularity whereas it

should be avoided that xq,p = 0 for all p and q where p > q.

51
4.7 Robustness to possible errors

In order to evaluate the performance of the proposed methods for estimating

reflection, the methods are tested upon numerical data, where possible errors
are introduced.

The possible errors are:

• Random noise on signals

• Non-linear effects in the waves

• Phase locked waves

• Three-dimensional waves in cases where a two-dimensional method is

used

• Divergens in the estimation due to numerical problems in solving the

system of equations

4.7.1 Random noise on signals

The following section describes a numerical test, where the elevations are
calculated for monochromatic waves according to the shown equation:

η(x, t) = aI cos(kx − ωt) + αaI cos(kx + ωt) + random · β · aI

where α is reflection coefficient

β is noise coefficient
random is a random number in [−1; 1]

52
 Goda & Suzuki
β ω Δt aI α α̃12 α̃13 α̃23 α̃average
% rad/sec sec metre % % % % %
10 2π 0.0625 0.10 50.00 50.88 50.92 50.24 50.68
10 2π 0.0625 0.10 10.00 11.03 10.98 10.14 10.72

Goda & Suzuki

β ω Δt aI α α̃12 α̃13 α̃23 α̃average
% rad/sec sec metre % % % % %
30 2π 0.0625 0.10 50.00 52.26 52.34 51.13 51.91
30 2π 0.0625 0.10 10.00 12.78 13.65 11.57 12.67

Mansard & Funke

β ω Δt aI α α̃123 α̃average
% rad/sec sec metre % % %
10 2π 0.0625 0.10 50.00 50.50 50.50
10 2π 0.0625 0.10 10.00 10.42 10.42

Mansard & Funke

β ω Δt aI α α̃123 α̃average
% rad/sec sec metre % % %
30 2π 0.0625 0.10 50.00 51.49 51.49
30 2π 0.0625 0.10 10.00 11.32 11.32

where Δt is sampling interval

α is reflection coefficient
˜ indicates estimated

All tests were performed with water depth d = 0.50 metre with a 160 sec.
long time series. Distances between wave gauges were x1,2 = 0.25 metre and
x1,3 = 0.60 metre .

The tests show good restistance to random noise. Though, the method pro-
posed by Mansard & Funke gives sligthly better results than the method
proposed by Goda & Suzuki.

53
The tests were repeated for 15 different wave periods and all tests showed
the same tendency as described in the examples above.

4.7.2 Non-linear effects in the waves

All the described methods assume linear waves. In order to examine the
effect from non-linearites a 2-order wave were generated. First with a free
2-harmonic eq. (4.32), and then with a bounded 2-harmonics eq. (4.33):

η(x, t) = aI cos(kx − ωt) + 0.2 · aI cos(2k ∗ x − 2ωt) +

αaI cos(kx + ωt) + α · 0.2 · aI cos(2k ∗ x − 2ωt) (4.32)
η(x, t) = aI cos(kx − ωt) + 0.2 · aI cos(2kx − 2ωt) +
αaI cos(kx + ωt) + α · 0.2 · aI cos(2kx − 2ωt) (4.33)

where k is wave number corresponding to ω

k∗ is wave number corresponding to 2 · ω

Goda & Suzuki

Free 2-harmonics Bounded 2-harm
2−harm 1−harm 2−harm 1−harm 2−harm
ω Δt aI aI α α̃average α̃average α̃average α̃average
rad/sec sec metre metre % % % % %
π 0.0625 0.10 0.02 50.00 50.04 49.78 50.12 61.10
π 0.0625 0.10 0.02 10.00 10.22 10.18 10.18 26.33

Mansard & Funke

Free 2-harmonics Bounded 2-harm
2−harm 1−harm 2−harm 1−harm 2−harm
ω Δt aI aI α α̃average α̃average α̃average α̃average
rad/sec sec metre metre % % % % %
π 0.0625 0.10 0.02 50.00 50.00 53.14
π 0.0625 0.10 0.02 10.00 10.00 19.42

54
4.7.3 Numerical problems

Numerical problems due to the discretisation of the signal and the succeeding
FFT-analysis.

Timeseries were calculated from:

η(x, t) = aI cos(kx − ωt) + αaI cos(kx + ωt) (4.34)

where α is reflection coefficient

Goda & Suzuki

ω Δt aI α α̃12 α̃13 α̃23 α̃average
rad/sec sec metre % % % % %
2π 0.0625 0.10 50.00 50.21 50.32 49.84 50.12
2π 0.0625 0.10 10.00 10.28 10.45 9.80 10.18
0.98 · 2π 0.0625 0.10 50.00 48.92 49.78 53.00 49.20
0.95 · 2π 0.0625 0.10 50.00 53.88 34.64 53.32 42.98

Mansard & Funke

ω Δt aI α α̃123 α̃average
rad/sec sec metre % % %
2π 0.0625 0.10 50.00 50.00 50.00
2π 0.0625 0.10 10.00 10.00 10.00
0.98 · 2π 0.0625 0.10 50.00 49.66 49.66
0.95 · 2π 0.0625 0.10 50.00 48.70 48.70

All tests were performed with water depth d = 0.50 metre with a 80 sec.
long time serie. The time serie were divided into 10 sub-series, which were
cosine tapered. Distances between wave gauges were x12 = 0.25 metre and
x13 = 0.60 metre .

In order to examine errors from comming from a wrong calibration constant

on one or two of the wave gauges.

55
Timeseries were calculated from:

η(x, t) = aI cos(kx − ωt) + αaI cos(kx + ωt)

η(x1 , t) = 1.00 · aI cos(kx1 − ωt) + 1.00 · αaI cos(kx1 + ωt) (4.35)
η(x2 , t) = 1.05 · aI cos(kx2 − ωt) + 1.05 · αaI cos(kx2 + ωt)
η(x3 , t) = 1.10 · aI cos(kx3 − ωt) + 1.10 · αaI cos(kx3 + ωt)

Goda & Suzuki

ω Δt aI α α̃12 α̃13 α̃23 α̃average
rad/sec sec metre % % % % %
2π 0.0625 0.10 50.00 49.61 55.53 51.30 52.15
2π 0.0625 0.10 10.00 9.80 17.57 11.80 13.06

Mansard & Funke

ω Δt aI α α̃123 α̃average
rad/sec sec metre % % %
2π 0.0625 0.10 50.00 51.40 51.40
2π 0.0625 0.10 10.00 18.90 18.90

56
Chapter 5

Seperation of Incident and

Reflected Long-Crested Waves
Using Digital Filters

In the hydraulic laboratory environment a separation of an irregular wave

field into incident waves propagating towards a structure, and reflected waves
propagating away from the structure is often wanted. This is due to the fact
that the response of the structure to the incident waves is the target of the
model test.

Goda and Suzuki (1976) presented a frequency domain method for estimation
of irregular incident and reflected waves in random waves. Mansard and
Funke (1980) improved this method using a least squares technique.

In the following a time-domain method for Separating the Incident waves

and the Reflected Waves (SIRW-method) is presented. The method is based
on the use of digital filters and can separate the wave fields in real time.

5.1 Principle

To illustrate the principle of the SIRW-method the set-up shown in Fig. 5.1
will be considered. The surface elevation η(x, t) at a distance x from the wave
generator may be written as the sum of the incident and reflected waves: the
incident wave propagating away from the wave generator, and the reflected
wave propagating towards the wave generator. Even though the method
works for irregular waves it will be demonstrated in the following pages for

57
 Figure 5.1: Wave channel with piston-type wave generator.

the case of monochromatic waves.

η(x, t) = ηI (x, t) + ηR (x, t)

= aI cos(2πf t − kx + φI ) + aR cos(2πf t + kx + φR ) (5.1)

where
f : frequency
a = a(f ) : wave amplitude
k = k(f ) : wave number
φ = φ(f ) : phase

and indices I and R denote incident and reflected, respectively.

At the two wave gauges we have:

η(x1 , t) = aI cos(2πf t − kx1 + φI ) + aR cos(2πf t + kx1 + φR ) (5.2)

η(x2 , t) = aI cos(2πf t − kx2 + φI ) + aR cos(2πf t + kx2 + φR )
= aI cos(2πf t − kx1 − kΔx + φI ) +
aR cos(2πf t + kx1 + kΔx + φR ) (5.3)

where x2 = x1 + Δx has been substituted into Eq. 5.3.

It is seen that the incident wave is phase shifted Δφ = kΔx from signal
η(x1 , t) to signal η(x2 , t), and the reflected wave is phase shifted Δφ = −kΔx
due to opposite travel directions. These phase shifts are called the physical
phase shifts and are denoted φphys I and φphys
R , respectively.

58
The idea in the following manipulations of the elevation signals is to phase-
shift the signals from the two wave gauges in such ways that the incident
parts of the wave signals are in phase while the reflected parts of the signals
are in mutual opposite phase. In this case the sum of the two manipulated
signals is proportional to and in phase with the incident wave signal.

An amplification C and a theoretical phase shift φtheo are introduced into

the expressions for η(x, t). The modified signal is denoted η ∗ . For the i’th
wave gauge signal the modified signal is defined as:

η ∗ (xi , t) = CaI cos(2πf t − kxi + φI + φtheo

i )+
theo
CaR cos(2πf t + kxi + φR + φi ) (5.4)

This gives at wave gauges 1 and 2:

η ∗ (x1 , t) = CaI cos(2πf t − kx1 + φI + φtheo

1 )+
theo
CaR cos(2πf t + kx1 + φR + φ1 ) (5.5)

η ∗ (x2 , t) = CaI cos(2πf t − kx2 + φI + φtheo

2 )+
theo
CaR cos(2πf t + kx2 + φR + φ2 )

= CaI cos(2πf t − kx1 − kΔx + φI + φtheo

2 )+
theo
CaR cos(2πf t + kx1 + kΔx + φR + φ2 ) (5.6)

The sum of η ∗ (x1 , t) and η ∗ (x2 , t), which is denoted η calc (t), gives:

η calc (t) = η ∗ (x1 , t) + η ∗ (x2 , t)

= CaI cos(2πf t − kx1 + φI + φtheo

1 )+
theo
CaR cos(2πf t + kx1 + φR + φ1 ) +
CaI cos(2πf t − kx1 − kΔx + φI + φtheo
2 )+
theo
CaR cos(2πf t + kx1 + kΔx + φR + φ2 )

= 2CaI cos(0.5(−kΔx − φtheo

1 + φtheo
2 ))
cos(2πf t − kx1 + φI + 0.5(−kΔx + φtheo1 + φtheo
2 )) +
2CaR cos(0.5(−kΔx + φtheo
1 − φtheo
2 ))
cos(2πf t + kx1 + φR + 0.5(kΔx + φtheo
1 + φtheo
2 )) (5.7)

59
It is seen that η calc (t) and ηI (x1 , t) = aI cos(2πf t − kx1 + φI ) are identical
signals when the following three conditions are met:

2C cos(0.5(−kΔx − φtheo
1 + φtheo
2 )) = 1 (5.8)
0.5(−kΔx + φ1 + φ2 ) = n · 2π
theo theo
n ∈ (0, ±1, ±2, ..) (5.9)
π
0.5(−kΔx + φtheo
1 − φtheo
2 ) = +m·π m ∈ (0, ±1, ±2, ..) (5.10)
2

Solving Eqs. 5.8 - 5.10 with respect to φtheo

1 , φtheo
2 and C gives Eqs. 5.11 -
5.13. n and m can still be chosen arbitrarily.

φtheo
1 = kΔx + π/2 + mπ + n2π (5.11)
theo
φ2 = −π/2 − mπ + n2π (5.12)
1
C = (5.13)
2 cos(−kΔx − π/2 − mπ)

All the previous considerations and calculations were done in order to find
an amplification and a phase shift for each of the two elevation signals η1
and η2 . Eqs. 5.11 - 5.13 give the result of our efforts, i.e. ηI (x1 , t) = η calc (t).
Remembering that φtheo 1 = φtheo
1 (f ), φtheo
2 = φtheo
2 (f ) and C = C(f ), it is
seen that the goal is already reached in the frequency domain. However, the
implementation of the principle will be done in the time domain using digital
filters. It is seen that singularities may occur. The consequences and the
handling of the singularities will be treated later on in the chapter. Here
it should just be mentioned that one way to bypass the singularities is to
use a velocity meter instead of one of the two wave gauges. Such system is
dealt with in the following chapter while the present chapter deals with two
surface elevation signals only.

η(x1, t) −→ FILTER 1


+ −→ ηI (x1, t)

η(x2, t) −→ FILTER 2


Figure 5.2: Flow diagram for signals in the SIRW-method.

60
The purposes of the filters shown in Fig. 5.2 are exactly a frequency depen-
dent amplification and a frequency dependent phase shift on each of the two
elevation signals.

Taking n = 0 and m = 0 the frequency response functions H1 (f ) for filter 1

and H2 (f ) for filter 2 calculated due to Eqs. 5.11 - 5.13 are given below in
complex notation:

1
Re{H1 (f )} = · cos(kΔx + π/2)
2 cos(−kΔx − π/2)
1
Im{H1 (f )} = · sin(kΔx + π/2) (5.14)
2 cos(−kΔx − π/2)
1
Re{H2 (f )} = · cos(−π/2)
2 cos(−kΔx − π/2)
1
Im{H2 (f )} = · sin(−π/2) (5.15)
2 cos(−kΔx − π/2)

Based on Eqs. 5.14 and 5.15 it is straightforward to design the time domain
filters. The design of the filters will be given on the next pages.

5.2 Design of Filters

The impulse response of the filters is found by an inverse discrete Fourier

transformation, which means that N discrete values of the complex frequency
response are used in the transformation, see Fig. 5.3. This gives an impulse
response of finite duration, i.e. the impulse response hj or the filter coeffi-
cients are found by:

1  r
N −1
2πr j
h = h(j · Δtf ilter ) =
j
H · exp i (5.16)
N r=0 N

where
r = 0,. . . , N -1
j = 0,. . . , N -1

and H r is the complex frequency response given by Eqs. 5.14 and 5.15 at the
frequency f = r · Δf . Note that H r should be hermitian, i.e. H(fn + fi ) =
H ∗ (fn − fi ).

61
 |H(f)|
10.0

8.0

6.0

4.0

2.0

0.0
0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0
frequency (Hz)

Figure 5.3:
Magnitude (gain) of the frequency responses of a discrete filter.
N = 64 , h = 0.5 m, Δf =0.10 Hz, Δtf ilter = 0.16 sec.,Δx =
0.2m.

The frequency increment, Δf , in the frequency response is found by

1
Δf = (5.17)
N · Δtf ilter

where Δtf ilter is the time increment of the filter.

Fig. 5.4 gives an example on the filters. The price paid for handling only N
frequencies in this transformation, is a minor inaccuracy in the performance
of the filter at input frequencies, which do not coincide with one of the cal-
culated frequencies in the discrete filter.

If the length of the filter (N ) is increased, more frequencies are included,

and in principle the overall accuracy of the filter is improved. In practice,
however, there is a limit beyond which the accuracy of the filter starts to
decrease due to other effects in the model.

The convolution integral (summation), Eq. 5.18, describes the input-output

relationship for the filters. Notice that the output η ∗ (x, t) is delayed (N -
1)/2 time steps relative to the input η(x, t).

62
 h(t), Filter 1
2.0

1.0

0.0

-1.0

-2.0
h(t), Filter 2
2.0

1.0

0.0

-1.0

-2.0
0.0 2.0 4.0 6.0 8.0 10.0
time (sec.)

Figure 5.4:
Filter coefficients corresponding to Filter 1 and Filter 2. N = 64,
water depth h = 0.5 m, Δf =0.10 Hz, Δtf ilter = 0.16 sec.,Δx =
0.2m.


N −1
∗p
η = hj · η p−j (5.18)
j=0

where
j,p = 0,. . . , N -1
η p−j : elevation at time t = (p − j) · Δtf ilter
η ∗p : output from filter at time t = p · Δtf ilter
hj : the filter coefficient corresponding to time t = j · Δtf ilter

Fig. 5.3 indicates that in the present example, singularities are present at
frequencies of about 2.0 Hz and 2.8 Hz. The figure also shows that due to

63
the fact that the frequency response is calculated only at discrete frequencies
in the filters, the singularities will not destroy the calculations. However,
it is recommended to cut off the frequency responses whenever the value is
larger than around 5. For practical use this means that, if |H(f )| ≥ 5 when
calculated, then |H(f )| should be valued 5. Furthermore, it is recommended
to place the singularities in a frequency range where the wave spectrum is
without significant energy, for example 3 times the peak frequency of the
spectrum. This can always be done by choosing appropriate values of Δx
and Δtf ilter , i.e. Δx smaller than a quarter of the shortest wave lengths.

5.2.1 Numerical Test Examples

In order to evaluate the SIRW-method we will look at two numerical examples

with known incident and reflected waves. The error is described by the
difference between the calculated incident wave signal η calc and the actual
incident wave signal ηI .

In the examples the total elevation due to two superimposed sine waves is
described by Eq. 5.19, corresponding to 50 % reflection of the incident waves.

η(x, t) = 0.01 · cos(2πf1 t − k1 x) + 0.01 · cos(2πf2 t − k2 x) +

0.01 · 0.5 · cos(2πf1 t + k1 x) +
0.01 · 0.5 · cos(2πf2 t + k2 x) (5.19)

The signals are sampled with a frequency of 6.4 Hz. Fig. 5.5 illustrates the
functionality of the method, when f1 and f2 are both coinciding with some
frequencies of the discrete filter, i.e. n · Δf .

64
 η (cm) ηI
η calc
3.0 ηx1

1.5

0.0

-1.5

-3.0

Error = η calc − ηI (cm)

3.0
Error
1.5

0.0

-1.5

-3.0
0.0 5.0 10.0 15.0 20.0
time (sec.)

Figure 5.5: A comparison between ηI , η calc and ηx1 . f1 = 4Δf , f2 = 7Δf .

As expected the method is exact for signals consisting only of energy placed
at the discrete frequencies (Fig. 5.5), though it is seen that errors are present
during warm up of the filters.

The second example (Fig. 5.6) is identical to the first example except that f1
and f2 are not coinciding with frequencies in the digital filter i.e. f1 = 4.2Δf ,
f2 = 7.5Δf . It must be stressed that the output signal shown in Fig. 5.6
corresponds to the worst case situation, where the wave frequencies are placed
midway between filter frequencies.

65
 η (cm) ηI
η calc
3.0 ηx1

1.5

0.0

-1.5

-3.0

Error = η calc − ηI (cm)

3.0
Error
1.5

0.0

-1.5

-3.0
0.0 5.0 10.0 15.0 20.0
time (sec.)

Figure 5.6: A comparison between ηI , η calc and ηx1 . f1 = 4.2Δf , f2 =

7.5Δf .

One way to improve the results is to apply a tapering of the filter coefficients,
because the output from a digital filter is more stable in case the absolute
values of the filter coefficients are almost zero in both ends of the filter, Karl
(1989). Cosine tapering of the filter coefficients improves the accuracy of the
SIRW method as demonstrated in Fig. 5.7.

66
 η (cm) ηI
η calc
3.0 ηx1

1.5

0.0

-1.5

-3.0

Error = η calc − ηI (cm)

3.0
Error
1.5

0.0

-1.5

-3.0
0.0 5.0 10.0 15.0 20.0
time (sec.)

Figure 5.7: A comparison between ηI , η calc and ηx1 . The filters have been
cosine tapered. f1 = 4.2Δf , f2 = 7.5Δf .

5.2.2 Physical Model Tests

The SIRW-method previously described was also tested in a laboratory flume

at the Hydraulics and Coastal Engineering Laboratory, Aalborg University,
cf. Fig 5.8.

First, the waves (incident part of the timeseries) were generated and sent
towards a spending permeable beach (slope 1:8) with low reflection (app. 5
%) in order to obtain a good estimate of the incident waves. Next, a reflecting
wall was mounted in the flume giving a fairly high reflection (app. 50 %) and
the same incident waves were reproduced by play back of the same digital
steering signal to the wave maker. Notice, that the incident wave fields are
identical only until re-reflection occurs. In Fig. 5.9 the output from the
SIRW-filters is compared with the incident waves measured in the case of
very low reflection.

67
 Figure 5.8: Set-up for physical model tests.
η (cm) ηImeasured
η calc
3.0 ηxmeasured
1

1.5

0.0

-1.5

-3.0

Error (cm)

3.0
η calc − ηImeasured
1.5 η calc − ηxmeasured
1

0.0

-1.5

-3.0
10.0 15.0 20.0 25.0 30.0
time (sec.)

Figure 5.9: A comparison between ηImeasured , η calc and ηxmeasured

1
. f1 = 4.2Δf ,
f2 = 7.5Δf .

The specific part of the signals, where reflection is present but re-reflection

68
from the wave paddle is still not present, is shown. Two different estimates
of the incident waves are used, namely the measured elevation at gauge no
1 (ηxmeasured
1
) and the calculated elevation at gauge no 1 (η calc ). In the the
specific example the SIRW-method reduces the error (variance of the differ-
ence η calc −ηx1 ) from 30 % of the incident energy to 3 % of the incident energy.

5.3 Conclusions

A time-domain method for Separating Incident and Reflected Irregular Waves

(The SIRW-method) has been presented.

By numerical and physical model tests it is demonstrated that the method

is quite efficient in separating the total wave field into incident and reflected
waves. Please note, that all the tests shown were done with fairly small filters
(few filter components), and that longer filters will improve the efficiency of
the method. Taking the example shown in Fig. 5.6 and doubling the number
of filter coefficients the error (variance) will decrease to 2/3 of the shown
example.

The accuracy of the SIRW-method is comparable with the accuracy of the

method proposed by Goda and Suzuki (1976), but the SIRW-method has the
advantage that where the incident wave signal is wanted in time domain (i.e.
for zero-crossing analysis) the singularity points are treated more properly
than in the Goda-method. The SIRW-method can easily be extended to give
the same accuracy as the method proposed by Mansard and Funke (1980).

The greatest advantage of the SIRW-method is that it works in real time.

Brorsen and Frigaard (1992) previously used digital filters to make an open
boundary condition in a Boundary Element Model, based on a filtering of
the surface elevation. That boundary condition accumulated errors, because
separation of the surface elevation into incident and reflected waves was not
possible in real time at that moment and, consequently, the Boundary El-
ement Model became unstable and could only run for a limited time. The
SIRW-method will make it possible to use digital filters as boundary condi-
tion in these models.

At the moment the SIRW-method is implemented at Aalborg Hydraulics

Laboratory, Aalborg University and the method is used in active absorption.

69
Chapter 6

Reflection Analysis of Oblique

Long-Crested Waves

In the case of oblique waves, i.e. 2D-waves travelling along a line not perpen-
dicular to the reflecting structure, the previous derived methods can all can
be applied by assuming that the waves can be decomposed into two vectorial
components, being respectively ⊥ and || to the structure.

N
η(x, t) = n=1 aI,n cos(kn x − ωn t + ΦI,n ) +
N
n=1 aR,n cos(kn x + ωn t + ΦR,n ) (6.1)

where η is the surface elevation relative to MWL

x is the vectorial position of the wave gauge
t is time
an is the amplitude
kn is the vectorial wavenumber
ωn is the angular frequency of the waves
Φn is the phase

71
Omitting index n and rewriting eq. (6.1) in cartesian coordinates:

η(x, y, t) = aI cos( kx cos(θI ) + ky sin(θI ) − ωt + ΦI ) +

aR cos( kx cos(θR ) + ky sin(θR ) + ωt + ΦR ) (6.2)

where x, y is the vectorial position of the wave gauge

k is the wavenumber
θ is direction of wave

Figure 6.1: Placement of wave gauges in front of structure.

It is seen that if the angle of the incident wave is know, and the angle of the
reflected wave is calculated using Snells law (incident angle = reflected an-
gle) eq. (6.3) is very similar to the original expression for the wave elevation,
eq(4.2).

The easiest way to solve the problem is to place the wave gauges on a line
perpendicular to the reflecting structure (y-coordinate ≡ 0 for all gauges) as
shown in fig 6. Eq. (6.2) will then simplify to:

η(x, 0, t) = aI cos(kx cos(θI ) − ωt + ΦI ) +

aR cos(kx cos(θR ) + ωt + ΦR ) (6.3)

72
Which shows the with this positioning of the gauges the normal methods
presented in the previous two chapters can be used when the geometrical
gauge distances are modified by the factor cos(θ).

Please notice that near the structure edge-waves will exist. These waves
travel along the structure and can destroy the calculations because they are
not included in the modelling of the waves.

73
Chapter 7

Methods for estimation of

directional spectra

A short crested wave field is normally referred to as a three-dimensional wave

field. This implies introduction of an additional parameter namely the di-
rection of travel. A two dimensional wave field is commonly described in the
frequency domain by use of the wave spectrum, i.e. the auto spectrum of
the wave elevation process. Now in the three dimensional case a directional
spreading function depending on frequency and direction of travel is intro-
duced. The combination of the wave spectrum and the spreading function
is the directional wave spectrum. Having determined the directional wave
spectrum the wave field is fully described in the frequency domain.

For a two dimensional sea state the elevation was assumed to be a summation
of a number of wavelets as stated in e.g. eq. (4.1). In a three dimensional
case the corresponding expression is


N 
M
η(x, t) = amn cos(kmn x − ωn t + Φmn ) (7.1)
n=1 m=1

where k is the wavenumber vector. Eq. (7.1) implies that for each pair of
discrete values of frequency and direction of travel there exist a long crested
wave propagating with these properties.

The amplitude amn can alternatively be described by use of the wavespectrum

75
as the variance of a sinusoidal wave is half the amplitude squared. Hence

1 2
a = S(ωm , θn )ΔθΔω
2 mn 
amn = 2S(ωm , θn )ΔθΔω (7.2)

Inserting eq. (7.2) into eq. (7.1) yields


N 
M

η(x, t) = 2S(ωm , θn )ΔθΔω cos(kmn x − ωn t + Φmn ) (7.3)
n=1 m=1

Letting Δθ → dθ and Δω → dω eq. (7.3) turns into a double integral

! ∞ ! π 
η(x, t) = cos(kx − ωt + Φ(ω, θ)) 2S(ω, θ)dθdω (7.4)
0 −π

The aim of the following is to establish an expression which relates known

and measured numbers to the directional wave spectrum or the directional
spreading function. These latter functions are related through

S(ω, θ) = H(ω, θ)S(ω) (7.5)

Two wave gauges are considered. These are positioned as shown in fig. 7.

6 Gauge 2
e

r12
Direction of travel
 
β12 
*

 
 θ
 

`
e - x
Gauge 1

Figure 7.1: Definition of geometric parameters.

76
The elevations at gauge 1 and gauge 2, respectively, are
! ∞ ! π 
η1 = η(x1 , t) = cos(kx1 − ωt + Φ(ω, θ)) 2S(ω, θ)dθdω
−π
! 0
∞! π 
η2 = η(x2 , t) = cos(kx2 − ωt + Φ(ω, θ)) 2S(ω, θ)dθdω
−π
! 0
∞! π
= cos(kx1 − ωt + kr12 cos(θ − β12 ) + Φ(ω, θ)) ·
0 −π

2S(ω, θ)dθdω

The cross-correlation function is obtained as

! T
1
Rη1 η2 (τ ) = η1 (t)η2 (t + τ )dt (7.6)
T 0

Here it is necessary to specify the time argument, thus η1 (t) = η1 and η2 (t) =
η2 . Inserting into eq. (7.6) yields

! ! !
1 T ∞ π
Rη1 η2 (τ ) = cos(kx1 − ωt + Φ(ω, θ)) ·
T 0 0 −π
cos(kx1 − ω(t + τ ) + kr12 cos(θ − β12 ) + Φ(ω, θ)) ·
2S(ω, θ)dθdωdt (7.7)

Applying the trigonometric relations

2 cos α cos β = cos(α − β) + cos(α + β)

cos(α ± β) = cos α cos β ∓ sin α sin β

leads to
! ! !
1 T ∞ π
Rη1 η2 (τ ) = [cos(kr12 cos(θ − β12 ) − ωτ )
T 0 0 −π
+ cos(2kx1 + kr12 cos(θ − β12 ) − 2ωt − ωτ + 2Φ(ω, θ))] ·
S(ω, θ)dθdωdt
! ∞! π
= cos(kr12 cos(θ − β12 ) − ωτ )S(ω, θ)dθdω
0 −π

77
which can be written as
! ∞ ! π
Rη1 η2 (τ ) = [cos(ωτ ) cos(kr12 cos(θ − β12 ))
0 −π
+ sin(ωτ ) sin(kr12 cos(θ − β12 ))] S(ω, θ)dθdω (7.8)

Because the cross-correlation is the inverse Fourier transform of the cross-

spectrum, the cross-correlation can also be written
! ∞ ! ∞
Rη1 η2 (τ ) = c12 (ω) cos(ωτ )dω + q12 (ω) sin(ωτ )dω (7.9)
0 0

where c12 is the co-spectrum and q12 the quad-spectrum in the cross-spectrum

Sη1 η2 (ω) = c12 (ω) − iq12 (ω) (7.10)

A comparison between eq.(7.8) and eq.(7.9) leads to

! π
c12 (ω) = S(ω, θ) cos(kr12 cos(θ − β12 ))dθ
!−π
π
q12 (ω) = S(ω, θ) sin(kr12 cos(θ − β12 ))dθ
−π

Inserting these expressions for the co- and quad-spectrum into eq. (7.10)
yields
! π
Sη1 η2 (ω) = S(ω, θ) exp(−ikr12 cos(θ − β12 ))dθ
!−π
π
Sη1 η2 (ω)
= H(ω, θ) exp(−ikr12 cos(θ − β12 ))dθ (7.11)
Sηη (ω) −π

Eq. (7.11) relates data, which can be computed from wave elevation mea-
surements, to the unknown directional spreading function. The geometry is
represented by r12 and β12 .

Eq. (7.11) is the basic relation that was in demand. It cannot be solved
analytically. Therefore methods have been proposed in order to make a
reliable estimate.

Common for all methods is that they make assumptions about the shape
of the directional spreading function. That may either be a parameterized

78
analytic expression or a number of discrete values making a step curve. This
initial assumption is very important to the ability of the particular method.
The next step is to fit the unknown parameters to the measured data. The
reliability of the fit depends on the reliability of the measured data but also on
the number of coefficients to fit, i.e. the more coefficients the less reliability.
This is a problem as a high resolution requires many coefficients.

7.1 The Maximum Likelihood Method

In the Maximum Likelihood Method, MLM, the directional spectrum is

calculated from minimizing the errors between measured wave data (cross-
correlations) and fitted directional spectrum.

This is achieved by use of the Maximum Likelihood technique, which has

named the method.

In it most simple form MLM applies for simultaneous wave elevations, but
by use of transfer functions it can readily be extended to apply for various
measurements. The MLM was originally presented by Capon (1969), but has
since been modified by several authors. Especially the papers by Davis and
Regier (1977) and Isobe et al. (1984) are of importance.

As an initial assumption the directional spectrum is expressed as a linear

combination of the cross-spectra. The linear wave theory is assumed valid as
well, i.e.

k = 2π/L
L = L0 tanh(kd)

Hence the estimated directional spectrum, S̃(ω, θ), can be written as


N 
N
S̃(ω, θ) = αmn (ω, θ)Smn (ω) (7.12)
m=1 n=1

Then using

! 2π
Smn (ω) = exp(−ik(xn − xm ))S(ω, θ)dθ (7.13)
0

79
eq. (7.12) leads to


N 
N ! 2π
S̃(ω, θ) = αmn (ω, θ) exp(−ik (xn − xm ))S(ω, θ )dθ(7.14)


m=1 n=1 0

$ %
 k cos θ
where k = .
k sin θ

Eq. (7.14) is expressed as

! 2π
S̃(ω, θ) = w(ω, θ, θ )S(ω, θ )dθ (7.15)
0

where


N 
N

w(ω, θ, θ ) = αmn (ω, θ) exp(−ik (xn − xm )) (7.16)
m=1 n=1

From eq. (7.15) it is seen that the estimated directional spectrum is a convo-
lution of the true directional spectrum and the window function w(ω, θ, θ )
given by eq. (7.16).

It is assumed that the coefficients αmn (ω, θ) can be expressed as

αmn (ω, θ) = γm (ω, θ)γn∗ (ω, θ) (7.17)

80
Inserting this into eq. (7.12) and eq. (7.16) yields


N 
N
S̃(ω, θ) = γm (ω, θ)γn∗ (ω, θ)Smn (ω) (7.18)
m=1 n=1
N  N
w(ω, θ, θ ) = γm (ω, θ)γn∗ (ω, θ) exp(−ik (xn − xm ))
m=1 n=1
N  N
= γm (ω, θ)γn∗ (ω, θ) exp(ik xm ) exp(−ik xn )
m=1 n=1
N 
N
= γm (ω, θ) exp(ik xm ) γn∗ (ω, θ) exp(−ik xn )
m=1 n=1
N N
∗
= γm (ω, θ) exp(ik xm ) (γn (ω, θ) exp(ik xn ))
m=1 n=1

N
= | γm (ω, θ) exp(ik xm )|2 (7.19)
m=1

The window function is normalised by setting w(ω, θ, θ) = 1. Having done

this it is seen from eq. (7.15) that if the window function is Diracs delta func-
tion, then the estimated directional spectrum is equal to the true directional
spectrum. As both the window function and the true directional spectrum
is non-negative the aim is to minimise S̃(ω, θ) as given in eq. (7.18) i.e.


N 
N
minimise γm (ω, θ)γn∗ (ω, θ)Smn (ω)
m=1 n=1

The same problem can be formulated as

  N N ∗
w(ω, θ, θ) n=1 γm (ω, θ)Tmn (ω, θ)γn (ω, θ)
maximise = Nm=1
N (7.20)
S̃(ω, θ) ∗
m=1 n=1 γm (ω, θ)Smn (ω, θ)γn (ω, θ)

where a matrix T has been introduced as

Tmn (ω, θ) = exp(ikxm ) exp(−ikxn ) (7.21)

∗
= γom (ω, θ)γon (ω, θ) (7.22)

The maximum value of eq. (7.20) is equal to the maximum eigenvalue of the
matrix S −1 T . Eq. (7.20) is then the Rayleigh quotient of γ T (T − λS)γ ∗ = 0

81
whereas

γT T γ∗
maximum = λmax (7.23)
γ T Sγ ∗

−1
where λmax is the maximum eigenvalue. Thus by multiplication of γ T and
S −1 and substitution of T

γT T γ∗ = λmax γ T Sγ ∗
T γ∗ = λmax Sγ ∗
S −1 T γ ∗ = λmax γ ∗
S −1 γ ∗o γ To γ ∗ = λmax γ ∗

Further by premultiplication by γ To yields

γ To S −1 γ ∗o γ To γ ∗ = λmax γ To γ ∗ (7.24)

Thus

γ To S −1 γ ∗o = λmax (7.25)

It is given from eq. (7.20) that

S̃(ω, θ) ∝ 1/λmax

Hence introducing a proportionality factor κ the directional spectrum can be

estimated as

& '−1
S̃(ω, θ) = κ γ To S −1 γ ∗o (7.26)

The terms in the RHS of eq. (7.26) are all known when a sample has been
carried out. S is the cross-spectrum matrix and is calculated from the time
series. γo is dependent only on the wavenumber vector and the geometry.
The factor κ is used in order to achieve the correct variance (i.e. the mea-
sured variance). While the MLM (and BDM) provides reasonable results for
directional wave spectra containing incident waves only the methods becomes
less reliable when reflections occur. This problem happen because there is a
corelation of the phases of the incident and reflected waves. This problem

82
arises in most methods.

Several methods has been suggested to handle reflection. Common ap-

proaches is to either provide information about the location and direction
of the reflector i.e. Isobe & Kondo (1984) or by assuming a parametric form
of the spreading functions as suggested by Yokoki, Isobe & Watanabe (1992).
Hiromune et al. (1992) have done this utilising Mitsuyasu’s spreading func-
tion. Further the spreading function has the option to imply a reflected wave
system which in spreading has the same form as the incoming system but
reduced in energy. This method assumes that the waves are reflected along
a structure that has a straight front.

7.2 MLM utilising standard spectra

The purpose of the present section is to present the Maximum Likelihood

method for estimating directional spectra utilising standard spectra. The
presentation is based on Isobe (1990), Yokoki, Isobe & Watanabe (1992)
and Christensen & Sørensen (1994). The directional spectrum is given in a
standard form in terms of some unknown parameters to be estimated from
measured data. In the present section only surface elevations measurements
are treated.

The starting point is M surface elevations, η(x, t), measured at M different

locations x at time t. The total elevation processes ηp (xp , t) , p = 1, 2, ..., M ,
are modelled as stochastic processes. The processes are assumed to be joint
stationary, ergodic and Gaussian distributed. The mean value functions
μηp (t), p = 1, 2, . . . , M , are assumed to be equal to 0. The M time series
can be written as a Fourier-sum as


N
ηp (x, t) = (Ap,l cos ωl t + Bp,l sin ωl t) , p = 1, 2, ..., M (7.27)
l=1

where ωl = l 2π
T
and T is the length of the time series. The coefficients Ap,l

83
and Bp,l are given as the stochastic integrals

! T
2 2
Ap,l = ηp (xp , t) cos ωl t dt l = 1, 2, . . . , N p = 1, 2, . . .(7.28)
,M
T − T2
! T
2 2
Bp,l = ηp (xp , t) sin ωl t dt l = 1, 2, . . . N p = 1, 2, . . . (7.29)
,M
T − T2

In eq. (7.27) the term corresponding to l = 0 has been omitted as it equals

0.

From eq. (7.28) and eq. (7.29) it is seen that all the coefficients Ap,l and Bp,l
are joint Gaussian distributed stochastic variables.

In the following only the l’th components are considered. The primary
aim is to determine the joint distribution of the coefficients Ap,l and Bp,l
p = 1, 2, . . . , M , at the frequency wl .

The reason for limiting the analysis to a single frequency at a time is that the
unknown parameters in S(ω, θ) to be estimated may be frequency dependent.

The coefficients are expressed as a vector

AT = [A1 A2 . . . AM B1 B2 . . . BM ]
= [A1 . . . AM AM +1 . . . A2M ] (7.30)

As the coefficients are joint Gaussian distributed the frequency distribution

function is given in terms of the mean value function
  vector, E[A], and the
cross covariance function matrix κA AT = E A A . The cross covariance
T

function matrix of size 2M × 2M is symmetric and of the form

 
B E
κ A AT = (7.31)
ET D

where the submatrices B, E and D are M × M matrices.

84
From eq. (7.28) and eq. (7.29) the following results are found

! 2
2 T
E [Ap,l ] = E [ηp (xp , t)] cos ωl t dt = 0 , p = 1, 2, . . . , M (7.32)
T − T2

and

! 2
2 T
E [Bp,l ] = E [ηp (xp , t)] sin ωl t dt = 0 , p = 1, 2, . . . , M (7.33)
T − T2

i.e. E [A] = 0.

 
In the following the cross covariance function matrix E A AT will be de-
termined. As an example the submatrix B in eq. (7.31) is considered.

( ! 2 ! 2
)
2 T 2 T
κAm,l An,l (t1 , t2 ) = E ηm (t1 ) cos (ωl t1 ) dt1 · ηn (t2 ) cos (ωl t2 ) dt2
T − T2 T − T2

! 2 ! 2
4 T T
= E [ηm (t1 ) ηn (t2 )] cos (ωl t1 ) cos (ωl t2 ) dt1 dt2
T2 − T2 − T2

! 2 ! 2
4 T T
= κmn (t2 − t1 ) cos (ωl t1 ) cos (ωl t2 ) dt1 dt2 (7.34)
T2 − T2 − T2

It is seen, that κAm n depends on the cross covariance between the elevation
cj Acj
processes ηm (t1 ) and ηn (t2 ). The cross covariance matrix κmn (τ ) is related
to the cross spectral density matrix Smn (ω) in terms of the Wiener-Khinchine
relation
! ∞
Gmn (ω) = 2Smn (ω) = 2 κmn (τ ) e−iωτ dτ
−∞
! ∞ ! ∞
= 2 κmn (τ ) cos ωτ dτ − i · 2 κmn (τ ) sin ωτ dτ
−∞ −∞

= Cmn (ω) − i Qmn (ω) (7.35)

85
where Gmn (ω) is a one-sided spectrum.

! ∞
Cmn (ω) = 2 κmn (τ ) cos ωτ dτ
−∞

and
! ∞
Qmn (ω) = 2 κmn (τ ) sin ωτ dτ
−∞

are the co-spectrum and the quad-spectrum, respectively.

Assuming the period T to be ”very long” the following result is obtained

from eq. (7.34)

1 Δω
κAm,l An,l (t1 , t2 ) = κAm,l An,l (τ ) = Cmn (ωl ) = Cmn (ωl ) = Bmn
(7.36)
T 2π

It is concluded that B in eq. (7.31) equals Δω

2π
multiplied by the co-spectrum
of the elevation processes. Further B is time independent.

Using the same procedure the elements E and D of eq. (7.31) can be found
and after some calculation the following result is obtained
 
Δω C Q Δω
κAAT (ωl ) = = · Ω (ωl ) (7.37)
2π −Q C 2π

i.e. the cross-covariance between the coefficients in eq. (7.27) at a given

frequency ωl is given in terms of the co- and quad-spectra of the elevation
processes.

So far expressions have been established for the mean value vector (eq. (7.32)
and eq. (7.33)) and the cross covariance matrix eq. (7.37) of the vector A.
These relations were established at a given frequency ωl . Besides the fre-
quency the wave pattern is also characterised by a direction of travel. The
dependence of the frequency and the direction of travel is expressed in terms
of the directional spectrum S(ω, θ). The following relation exists between

86
the one-sided autospectral density G(ω) and S(ω, θ)

! 2π
G(ω) = S(ω, θ) dθ (7.38)
0

Furthermore, it can be shown that S(ω, θ) is related to the cross spectral

density matrix, Gmn (ω), of the elevation processes ηm (xm , t) and ηn (xn , t) as
expressed in eq. (7.39)

 2π
Gmn (ω) = 0 S(ω, θ) { exp(ik (xm − xn ) + r(ω, θ) exp(ik (xm − xn T ))
+ r(ω, θ) exp(ik (xm T − xn )) + r2 (ω, θ) exp(ik (xm − xn ) T )}
(7.39)
dθ

where r(ω, θ) is the reflection coefficient. The modelling of reflected waves is

described in Christensen & S>rensen (1994). T is a transformation matrix
equal to
 
−1 0
T = (7.40)
0 1

If the directional spectrum S(ω, θ) was known eq. (7.39) may be used to cal-
culate Gmn (ω), e.g. by numerical integration.

Based on eq. (7.35) C and Q are identified as the real and the imaginary
parts of G. Finally κA AT (ωl ) can be determined from eq. (7.37). However,
S(ω, θ) is generally unknown. In the following a method will be presented
which can be used to determine a directional spectrum expressed in standard
form in terms of some unknown parameters. The method is known as the
Maximum Likelihood (ML) method. As expressed earlier the elements in
A have a joint normal distribution. The general expression for a density
function of a vector A with 2M joint Gaussian elements having a mean
value vector E(A) = 0 is

1 1 T −1
pA (a) = &√ '2M # # 1 exp − 2 a κAAT (τ )a (7.41)
2π #κ T (τ )#
2
AA

# #
where #κAAT # and κ−1 T are the determinant and the inverse matrix of
AA

87
κAAT respectively.

If S(ω, θ) was given eq. (7.41) could be used to calculate the probability of
the observed realisation a of A, where a represents the actual Fourier coef-
ficients obtained from a given time-series. Since S(ω, θ), or some parameters
in S(ω, θ), are unknown, a Likelihood function, L(·), will be formulated.
Expressed in terms of the Likelihood function the unknown parameters in
Sη (ω, θ) are determined as the values corresponding to the maximum value
of L.

In his article Isobe (1990) uses L time-series at each of the M locations

xp , p = 1, 2, . . . , M . Based on each of the L time-series an estimate
a(l) , l = 1, 2, . . . , L, of A(l) is obtained. The probability of a(l) is given
by eq. (7.41). Assuming the L observations to be independent the joint
probability
& (1) ' for &obtaining
' exactly& the' L observed estimates a(l) is given as
pA a · pA a (2)
· . . . · pA a(L) . Therefore Isobe (1990) suggested a
Likelihood function, L(·), as the L’th root of this product, i.e.

& ' * & ' & '+1/L

L a(1) , . . . , a(L) , G = pA a(1) · . . . · pA a(L)

1   −1
2M 2M
1
=  exp − Ω Ω̃lh (7.42)
(Δω)M det (Ω) 2 h=1 l=1 hl

where

2π 
L
Ω̃lh = aml amh (7.43)
Δω · L m=1

where Ω̃ is the measured cross spectral density matrix.

Eq. (7.42) represents the probability of obtaining exactly the estimates a(l) , l =
1, 2, . . . , L. The unknown quantity in eq. (7.42) is G (or S(ω, θ)).

The optimal choice of S(ω, θ) or (in practice) the unknown parameters in

S(ω, θ) are determined in order to maximise L(·), i.e. the optimal param-
eters maximise the probability of obtaining exactly the observed Fourier-
coefficients.

88
7.3 The Bayesian Directional spectrum esti-
mation Method

The Bayesian Directional spectrum estimation Method, BDM, is in principle

similar to MLM. However, BDM makes use of a Bayesian approach in order
to estimate the most likely estimate of the directional spectrum.

BDM has been presented by Hashimoto et al. (1987). It is assumed, that

the directional spreading function H(ω, θ) can be expressed as a piecewise-
constant function, which takes only positive values. The directional spreading
function is discretized into K intervals. Defining xl (ω) as

xl (ω) = ln H(ω, θl ) l = 1, 2, ..., K KΔθ = 2π (7.44)

H(ω, θ) can be approximated to


K $
1 (l − 1)Δθ ≤ θ < kΔθ
H(ω, θ)  exp(xl (ω))Il (θ) Il = (7.45)
0 otherwise
l=1

The relationship between the cross-spectrum and the directional spectrum

has been deducted to
! 2π
Smn (ω) = S(ω, θ) exp(−ikrmn cos(θ − βmn ))dθ (7.46)
0

Inserting the approximation eq. (7.45) and dividing by S(ω) leads to

! 2π 
K
Smn (ω)
 exp(xl (ω))Il (θ) exp(−ikrmn cos(θ − βmn ))dθ
S(ω) 0 l=1

K
 exp(xl (ω)) exp(−ikrmn cos(θl − βmn ))Δθ (7.47)
l=1

If M wave gauges are available, eq. (7.47) can be applied to N = M (M +1)/2

spectra of which M will be autospectra. However, considering the N complex
equations as two separate equations, i.e. the real term and the imaginary
term, 2N real equations are obtained. Expressing these in an arbitrary order

89
eq. (7.47) can be rearranged to


K
Sj (ω) = exp(xl (ω))αjl (ω) + εj (ω) j = 1, 2, . . . , M (7.48)
l=1

where an error εmn (ω) of the cross-spectrum Smn (ω) has been added, and

exp(−ikrmn cos(θl − βmn ))Δθ

αjl (ω) =  (7.49)
Smm (ω)Snn (ω)
S (ω)
Sj (ω) =  mn (7.50)
S(ω) Smm (ω)Snn (ω)
εmn (ω)
εj (ω) =  (7.51)
Smm (ω)Snn (ω)

where it is suggested that 1 ≤ j ≤ N denote real parts and N < j ≤ 2N

denote imaginary parts. However, the number of equations involved can be
altered at will.

It is assumed that εj : N (0; σ 2 ), hence when omitting the argument ω


K
ε j = Sj − exp(xl )αjl : N (0; σ 2 )
l=1

Sj ’s and αjl ’s are given. σ 2 and xl ’s are to be estimated.

The probability density function for εj is

1 ε2j
p(εj ) = √ exp − 2
2πσ 2σ

90
The likelihood function is then

,
2N
1 ε2j
L(ε1 , ε2 , ..., ε2N ; σ) = √ exp − 2
j=1
2πσ 2σ
 2N
1  ε2j
= exp −
(2πσ 2 )N j=1
2σ 2
⎛  2
⎞
1 −1 2N K
= 2 N
exp ⎝ 2 Sj − ⎠
exp(xl )αjl(7.52)
(2πσ ) 2σ j=1 l=1

L(x1 , x2 , ..., xK ; σ) = L(ε1 , ε2 , ..., ε2N ; σ) (7.53)

for the given Sj , j = 1, 2, ..., 2N .

As the estimate of H(ω, θ) becomes smoother, it is assumed, that the differ-

ences of second order of xl − 2xl−1 + xl−2 decrease, i.e. the value


K
(xl − 2xl−1 + xl−2 )2 (7.54)
l=1

where x0 = xK and x−1 = xK−1 , decreases as well.

Maximising the likelihood eq. (7.52) and minimising eq. (7.54) lead to the
estimate which also maximises

u2 
K
ln L(x1 , x2 , ..., xK ; σ) − (xl − 2xl−1 + xl−2 )2
2σ 2 l=1

where u is introduced as a hyperparameter. Further applying the exponential

function leads to the alternative expression


u2 
K
L(x1 , x2 , ..., xK ; σ) exp − 2 (xl − 2xl−1 + 2xl−2 )2 (7.55)
2σ l=1

At this stage the Bayesian approach is introduced utilising p(y|x) ∝ L(y|x)p(y).

When normalised the second term in eq. (7.55) can be regarded as the joint

91
distribution of x = (x1 , x2 , ..., xK )


u2 
K K
u
2
p(x|u , σ ) =2
√ exp − 2 (xl − 2xl−1 + xl−2 )2 (7.56)
2πσ 2σ l=1

This correspond to the prior distribution and is known, when an estimate of

x is given and u and σ have been estimated. Maximising eq. (7.55) leads to
the posterior distribution, when inserting into p(y|x) ∝ L(y|x)p(y), i.e.

ppost (x|u, σ 2 ) ∝ L(x, σ 2 )pprior (x|u, σ 2 ) (7.57)

Hence minimising the following quantity determine the value of x, which

maximises eq. (7.57)

 2 
1 
2N K
u 2 
K
Sj − exp(xl )αjl + 2 (xl − 2xl−1 + xl−2 )2
2σ 2 j=1 l=1
2σ l=1

I.e. σ can be omitted reducing the above quantity to

 2 

2N 
K 
K
Sj − exp(xl )αjl +u 2
(xl − 2xl−1 + xl−2 )2 (7.58)
j=1 l=1 l=1

During the derivation of BDM it has been presumed that

• All values in the spreading functions are larger than zero

• The spreading functions are smooth (the smoothness being dependent

on one parameter)

• The errors on the spectral estimates are outcomes of a Gaussian distri-

bution function

The BDM turns out to be a very useful and reliable method for estimation of
directional wave spectra. But some inaccuracy arises when the wave samples
contain reflections, which in worst cases results in non-converging results.
Similarly to the MLM it is possible to create additional constraints to the
shape of the directional spreading function or make some assumptions about
location and direction of the reflector.

92
7.4 Tests of the presented methods

In principle all methods for estimating reflected waves should be tested upon:

• Ability to process data correctly

• Robustness to possible errors

• Reliability of results from real wave data

In this section the performance of the methods are evaluated with numerical
data in situations where the reflection is known, and no possible noise (errors)
on the data are included.Moreover, the application of the BDM method to
data from physical model tests is demonstrated.

7.4.1 MLM utilising standard spectra on numerical

data

In the numerical tests performed, the generated data were simultaneous re-
alizations of surface elevation time series recorded in a CERC5 wave gauge
array with a radius of 1.0 m positioned at a water depth of 4.0 m, see fig. 7.2.

Two tests were performed: One in which no reflection occurred and one in
which a wall with a reflection coefficient of r=0.5 was positioned at x=0.0
m, see fig. 7.2.

93
 Figure 7.2: Wave gauge arrangement for numerical tests.

In both tests, the incident wave fields were irregular waves corresponding to
a Pierson-Moskowitz spectrum (fp = 0.5Hz, Hs = 0.5m) with a Mitsuyasu-
type spreading function (θ0 = π/3, s = 6). A sample frequency of fs = 4Hz
was applied.

Measured cross covariance matrices were determined from a total of 45 sur-

face elevation subseries each of length T =128 s.

The estimated value of the auto spectral density Sη (ω) and the maximum
likelihood estimates of the parameters θ0 , s and r are given in fig. 7.3 and
fig. 7.4. For comparison, the target values are plotted.

In both tests, the estimated values of the parameters are in good agree-
ment with the target values. However, in the test involving reflection the
method has some problems separating incident and reflected waves at some
frequencies resulting in larger estimated values of the auto spectral density,
see fig. 7.4 compared to the results obtained without reflection, see fig. 7.3.

94
 Autospectral density Peak wave direction
60m 6.3

50m 5.2

40m 4.2
teta [rad]
S [m*m]

30m 3.1

20m 2.1

10m 1.0

0m 0.0
0.3 0.4 0.6 0.7 0.8 1.0 1.1 0.3 0.4 0.6 0.7 0.8 1.0 1.1
frequency [Hz] frequency [Hz]

Degree of directional concentration Reflection coefficient

10.0
9.0 0.90
8.0 0.80
7.0 0.70
6.0 0.60
5.0 0.50
s

4.0 0.40
3.0 0.30
2.0 0.20
1.0 0.10
0.0 0.00
0.3 0.4 0.6 0.7 0.8 1.0 1.1 0.3 0.4 0.6 0.7 0.8 1.0 1.1
frequency [Hz] frequency [Hz]

Figure 7.3: Numerical test results (no reflection).

95
 Autospectral density Peak wave direction
60m 6.3

50m 5.2

40m 4.2

teta [rad]
S [m*m]

30m 3.1

20m 2.1

10m 1.0

0m 0.0
0.3 0.4 0.6 0.7 0.8 1.0 1.1 0.3 0.4 0.6 0.7 0.8 1.0 1.1
frequency [Hz] frequency [Hz]

Degree of directional concentration Reflection coefficient

10.0
9.0 0.90
8.0 0.80
7.0 0.70
6.0 0.60
5.0 0.50
s

4.0 0.40
3.0 0.30
2.0 0.20
1.0 0.10
0.0 0.00
0.3 0.4 0.6 0.7 0.8 1.0 1.1 0.3 0.4 0.6 0.7 0.8 1.0 1.1
frequency [Hz] frequency [Hz]

Figure 7.4: Numerical test results (with reflection).

96
7.4.2 BDM on numerical data

The BDM method has been tested utilising numerical simulations of an uni-
modal wave field and a bi-modal wave field.

Only the bi-directional wave field test is shown, see fig. 7.5. In this case
it is seen that both peaks are estimated well. Though the narrow peak is
estimated most accurate.

0.45
0.4 target
 BDM 
0.35  
0.3

0.25 
H(ω, θ)
 
0.2 
  
0.15  
   
0.1  
   


0.05  
0
0 50 100 150 200 250 300 350
θ

Figure 7.5: Test with simulated bi-directional wave field.

Data. Test 1.
Date of sample 070893
Depth of water 200m
Radius of array 26m
Sample frequency 4Hz
Duration 720sec.
Quantity Elevation
Autospectrum P–M,
Tp = 10s, Hs = 4.0m
Directional spreading function Mits.,
s1 = 5, s2 = 8, θ1 = −90◦ , θ2 = 90◦

97
7.4.3 BDM on Model Test Data

In fig. 7.6 a test from Aalborg University’s 3D wave basin is shown. The
main direction of the incident waves is 290 degrees and the mean direction
of the reflected waves should be 70 degrees.

0.8
generated
0.7 BDM  
 
0.6
 
0.5
 
H(ω, θ) 0.4
 
0.3

0.2     

 
  
0.1    
 
  
0
0 50 100 150 200 250 300 350
θ

Figure 7.6: Test with reflected waves in 3D wave basin. Test 2.

Data. Test 2.
Date of sample 090193
Depth of water 0.6m
Radius of array 0.25m
Sample frequency 10Hz
Duration 720sec.
Quantity Elevation
Autospectrum JONSWAP,
Tp = 0.8s, Hs = 0.1m
Directional spreading function Mits.,
truncated s1 = 3, θ1 = −90◦

The BDM method estimates the incident wave energy to be spread over
a wider range. This is reasonable considering the conditions in the basin,
i.e. limitations due to boundary conditions. The mean direction of the re-
flected waves though deviate about 20 degrees.

98
7.5 Conclusion

Several different methods for estimating reflection have been presented and
tested in the previous chapters.

The methods for estimating reflection are divided into groups:

• Frequency domain methods for 2-dimensional waves.

• Time domain methods for 2-dimensional waves.

• Frequency domain methods for 3-dimensional waves.

The frequency domain methods give the incident wave spectrum and the
reflected wave spectrum. The time domain methods give the incident waves
as function of time.

Some of the methods take into account the 3-dimensionality of the waves. If
waves on location are 3-dimensional it is necessary to use a method for 3-
dimensional waves, though because the methods for estimating reflection in
3-dimensional waves introduce the directional spreading function with many
degrees of freedom, they become less statistically reliable.

99
Chapter 8

Wave Groups

A wave group is generally defined as a sequence of consecutive high waves in

a random wave train.

In sea wave recordings, group formations of high waves occur from time to
time. This phenomenon corresponds to a non-zero correlation between suc-
cessive waves. Information concerning this correlation is of importance when
reproducing waves in the laboratory in order to determine the response of
the modelled structure. Normally, irregular waves are reproduced in accor-
dance with a specific energy spectrum solely defining the distribution of the
variances. The grouping of waves is determined by the distribution of the
phases. Hitherto, independence between successive waves have been applied
and the phases are treated as independent random variables, each with a
uniform probability density on the interval [0;2π] leading to a sea surface
that is Gaussian distributed. However, if the waves during wave propagation
become more non-linear there will be some coupling and thus dependence of
the phases of the component waves at different frequencies, which eventually
will modify the wave grouping.

To illustrate the effect of randomly assigned phases two wave trains are gener-
ated from the same energy spectrum. These two wave conditions are depicted
in Figure 8.1.

Figure 8.1 shows different groupiness characteristics, and clearly it is im-

portant to have information on the wave grouping when coastal structures
respond differently when exposed to the distinctive wave patterns. Espe-
cially, the stability of rubble mound structures appears to be significantly
affected by the wave grouping, but also the slow drift oscillations of moored
vessels is highly dependent on the wave grouping.

Burcharth (1979) and Johnson et al. (1978) found that the wave grouping

101
 5.0

2.5

Elevation (m)
0.0
0 100 200 300 400
30.0
-2.5
25.0
S(f)/m0 (m2*sec/m2)

-5.0
20.0
time (sec.)
15.0

10.0
5.0
5.0
2.5
0.0

Elevation (m)
0.00 0.05 0.10 0.15 0.20 0.25 0.30
frequency (Hz) 0.0
0 100 200 300 400

-2.5

-5.0

time (sec.)

Figure 8.1: Wave energy spectrum and generated grouped and non-grouped
wave trains.

significantly affects the stability of rubble mound breakwaters as well as the

run-up. Johnson et al., (1978) compared the effects of a grouped and a non-
grouped time series generated from the same energy spectrum, thus having
the same statistical properties. Conclusively, the model tests showed that
the breakwater response to the two different wave trains was quite different,
with the grouped wave train causing severe damage and the non-grouped
only causing minor rocking of the armour units. Similar significant influence
on the wave grouping was found in the tests performed by Burcharth (1979).

In irregular seas, model tests by Spangenberg (1980) showed that the wave
grouping has a significant influence on the slow drift motion of moored plat-
forms and vessels. This influence might be explained by the fact that the
period of the slow drift oscillations practically corresponds to the wave group
period where the wave grouping is pronounced.

Both examples illustrate the importance of a correct modelling of natural sea

waves in the laboratory if the structural responses are sensitive to the wave
grouping. A characterization of the wave grouping seems therefore evident.

8.1 Description of Wave Groups

A measure of the wave grouping is obtained by defining the wave envelope

to the time signal. Due to the presence of small waves in the signal the wave
envelope is difficult to determine. However, if the time signal is squared,
the squaring procedure will suppress the relative influence of the small waves
present, and furthermore, a slowly varying part appears which may be inter-

102
preted as the square envelope.

Assuming that the sea surface elevation at a given point is a realization of a

linear stationary Gaussian process defined by its one-sided spectrum Sη (f ),
it can be represented by an ordinary sum of a finite number of waves


N
η(t) = cn cos(ωn t + εn ) (8.1)
n=1

where cn = amplitude, ωn = cyclic frequency, and εn = phase angle. By

squaring the time signal following equation is obtained


N 
N
2
η (t) = cn cm cos(ωn t + εn ) cos(ωm t + εm ) (8.2)
n=1 m=1

N 
N
1
= {cn cm ( cos((ωn + ωm )t + (εn + εm )) +
n=1 m=1
2
1
cos((ωn − ωm )t + (εn − εm )))} (8.3)
2

Equation (8.3) represents a splitting of η 2 (t) into a slowly varying part (rep-
resented by the difference-frequencies) and a more rapid oscillating part (rep-
resented by the summation-frequencies).

By use of symmetry of the double summation, equation (8.3) can be expressed

in terms of four separate contributions

1 2 1 2
N N
2
η (t) = c + c cos(2ωn t + 2εn )
2 n=1 n 2 n=1 n

N 
N
+ cn cm cos((ωn + ωm )t + (εn + εm ))
n=1 m=n+1


N 
N
+ cn cm cos((ωn − ωm )t + (εn − εm ))) (8.4)
n=1 m=n+1

The four terms on the right-hand side of equation (8.4) are identified as fol-
lows: The first term consists of a constant off-set component. The second and
third term constitutes the superharmonic components, i.e. the summation-
frequency terms, and the fourth term constitutes the subharmonic compo-
nents, i.e. the difference-frequency terms. It is the latter that describes the

103
slowly varying part of the squared time signal and the term which may be
interpreted as the square envelope. By means of Bartlett filtering the su-
perharmonic components on the right-hand side of equation (8.4) may be
filtered out after subtraction of the constant off-set as done by Funke and
Mansard (1979).

Funke and Mansard denoted the filtered square of the time signal the SIWEH
(Smoothed Instantaneous Wave Energy History) function as the function
provides a measure of the instantaneous wave energy in the time signal.

The effect of the Bartlett filtering corresponds to a digital low pass filtering
and the efficiency of the SIWEH analysis can best be interpreted by examina-
tion of the energy spectrum of the stochastic process in (8.1) and the energy
spectrum of the squared process.

30.0 30.0

25.0 25.0
S(f)/m0 (m2*sec/m2)

S(f)/m0 (m2*sec/m2)

20.0 20.0

15.0 15.0

10.0 10.0

5.0 5.0

0.0 0.0
0.00 0.10 0.20 0.30 0.40 0.50 0.60 0.00 0.10 0.20 0.30 0.40 0.50 0.60
frequency (Hz) frequency (Hz)

Figure 8.2: a) JONSWAP energy spectrum for a linear stochastic process and
b) energy spectrum for the squared process.

From figure 8.2 it is understood that the SIWEH analysis does not exactly
isolate the slowly varying part; also contributions from the superharmonic
components occur and not the complete amount of energy from the subhar-
monic components is included. Only when the process is narrow-banded does
the SIWEH analysis perform well but as the process becomes more and more
broad-banded the SIWEH function is a poor estimator of the wave envelope,
see Hupspeth and Medina (1988).

Instead of using a Bartlett window to isolate the subharmonic components,

a wave envelope function defined on basis of the time series and its Hilbert
transform isolates exactly the subharmonic components.

104
8.2 Hilbert Transform Technique

From the sea surface elevation η(t) a conjugate signal η̂(t) is uniquely ob-
tained by shifting the phase of each elementary harmonic component of η(t)
by ± π2 . When the phase angles of all components of a given signal are shifted
± π2 , the resulting function η̂(t) is known as the Hilbert transform of the orig-
inal signal η(t). The Hilbert transform is defined by

! ∞
1 η(t)
η̂(t) = dτ (8.5)
π −∞ t−τ

From the definition of the Hilbert transform it is noted that η̂(t) is sim-
ply the convolution of η(t) with a linear filter with the impulse response
1 1
function h(t) = πt . Since a convolution of two functions in the time do-
main are transformed into a multiplication of their Fourier transforms in
the frequency domain 2 a frequency response function H(f ) is related to
the impulse response function. The frequency response function provides an
equally characterization of the linear time-invariant input and output system
in (8.5) and does furthermore visualize the effect of the Hilbert transform op-
eration. Through the Fourier transform the frequency response of the Hilbert
transformer becomes
⎧
1 ⎨ −i f >0
H(f ) = F[ ] = −i sgn(f ) = 0 f =0 (8.6)
πt ⎩
i f <0


The gain of this frequency response function is re2 (H(f )) + im2 (H(f ))
resulting in unity in magnitude, and thus, the amplitudes of the signal does
not change. The phase angle is arctan( im
re
(H(f ))
(Hf ))
) resulting in a phase angle of
− 2 for f > 0 and + 2 for f < 0. Such a system is denoted an ideal 90-degree
π π

phase shifter.

Consequently, applying the Hilbert transform operation to the sea surface

elevation in (8.1) the cosine function simply shifts to the sine function


N
η̂(t) = cn sin(ωn t + εn ) (8.7)
n=1

Associated with the Hilbert transform is the complex analytical signal defined
1
∞
The convolution of the functions g(t) ∗ h(t) is: g(t) ∗ h(t) ≡ −∞
g(τ )h(t − τ )dτ
2
The convolution theorem: g(t) ∗ h(t) ⇔ G(f )H(f )

105
from the original signal η(t) and the Hilbert transform η̂(t)

η̃(t) = | η̃(t) | exp(iψ(t)) = | η̃(t) | cos(ψ(t)) + i | η̃(t) | sin(ψ(t))

= η(t) + iη̂(t) (8.8)


where the envelope or the modulation | η̃(t) | = η 2 (t) + η̂ 2 (t) and the
η̂(t)
associated phase ψ(t) = arctan( η(t) ). The properties of the Hilbert transform
operation entail that the slowly varying difference-frequency terms in the
second order expression η 2 (t) are separated mathematically by the expression

E(t) ≡ re(η̃ ∗ (t)η̃(t)) =| η̃(t) |2 (8.9)

where η̃ ∗ (t) = the complex conjugate and E(t) = the square wave envelope
function.

In order to visualize the effect of the defined envelope function the Hilbert
transform of the sea surface elevation is squared and rewritten by use of
trigonometry and symmetry of the double summation similar to η 2 (t)


N 
N
2
η̂ (t) = cn cm sin(ωn t + εn ) sin(ωm t + εm ) (8.10)
n=1 m=1
N  N
1
= {cn cm ( cos((ωn − ωm )t + (εn − εm )) −
n=1 m=1
2
1
cos((ωn + ωm )t + (εn + εm )))} (8.11)
2
1 2 1 2
N N
= cn − c cos(2ωn t + 2εn )
2 n=1 2 n=1 n

N 
N
− cn cm cos((ωn + ωm )t + (εn + εm ))
n=1 m=n+1


N 
N
+ cn cm cos((ωn − ωm )t + (εn − εm ))) (8.12)
n=1 m=n+1

Remembering that the squared time signal is given by (8.4), the square wave
envelope function, according to (8.9), then becomes


N 
N 
N
E(t) = c2n +2 cn cm cos((ωn − ωm )t + (εn − εm )) (8.13)
n=1 n=1 m=n+1

106
Introducing √12 in the complex analytical signal η̃(t) = √12 (η(t) + iη̂(t)) leads
to the definition of an envelope function which may be interpreted as half
the square envelope.

1
E(t) =| η̃(t) |2 = (η 2 (t) + η̂ 2 (t)) (8.14)
2

This envelope function isolates exactly the slowly varying part of the squared
time signal plus the constant off-set similar to what approximately is achieved
by the SIWEH analysis.

The present method seems to be more convenient than the SIWEH analysis
and it does not require the narrow-band spectrum assumption. The disad-
vantage of this method is however that the sea surface must be described by
a linear model.

8.2.1 Computation of half the square envelope

To compute the Hilbert transform numerically the continuous-time convolu-

tion integral in (8.5) is approximated by a discrete-time Hilbert transforma-
tion. Furthermore, as the Hilbert transformation is non-banded, approxima-
tions limiting the impulse response function are made. A tool to handle the
ideal Hilbert transformation of the sea surface elevation is by using FIR ap-
proximations. In such approximations the 90-degree phase shift is conserved
exactly.

The principle in the FIR approximation is that the convolution integral

in (8.5) is represented by a summation over a finite number of coefficients
where the coefficients are fitted to represent the impulse response function.
Taking an even number of coefficients, easily extended to an odd number,
the non-causal FIR approximation can be written


Nc /2−1

N c −1

η̂j = ck ηj−k = ck ηj+k−Nc /2 (8.15)

k=−Nc /2 k=0

where ck = the k’th coefficient, Nc = number of coefficients or filter length,

η̂j is the Hilbert transform corresponding to the time step j, and ηj+k−Nc /2
are the input elevations to the filter system. The reason why the index on the
filter coefficients remain unchanged is that the coefficients are mirrored in the
Nyquist frequency, i.e. the frequency corresponding to half the filter length.
The coefficients are derived from the frequency response function by FFT

107
to obtain a least-square fit of the coefficients. Opposite the centered format
definition of the Fourier transformation, the FFT is based on a one-sided
format

1  1 
Nc −1 Nc −1
2πjk
ck = Xj exp(iωj kΔt) = Xj exp(i ) (8.16)
Nc j=0 Nc j=0 Nc

where ωj is the cyclic frequency corresponding to the j’th coefficient and Xj

is the desired sampled frequency response of the system. By using the one-
sided format a time delay corresponding to half the filter length is introduced

Nc
τ= Δt (8.17)
2

The corresponding phase delay may then be found as

2πj
ψτ = τ ωj = τ = πj (8.18)
Nc Δt

To compensate for the phase delay the original frequency response function
given by (8.6) only needs to be multiplied by a linear phase shift operator
exp(−iπk) and Xj might be interpreted as

Xj = H(fj ) exp(−iπj) = G(fj ) cos(ψj − πj) + iG(fj ) sin(ψj − πj)

(8.19)

where G(fj ) is the gain of the input amplitude to equal the output amplitude
and ψj − πj is the phase difference between the input and the output signal.

To sample the frequency response function the frequency band is subdivided

into Nc discrete frequencies where fj = j Nfsc and fs is the sample frequency.
Since the phase ψj = − π2 for 0 < fj < fN q and ψj = π2 for fN q < fj < 2fN q
the sampled discrete frequency response function becomes
⎧
⎨ G(fj ) cos(− π2 − πj) + iG(fj ) sin(− π2 − πj) 0 < fj < fN q
H(fj ) = 0 fj = 0, fN q (8.20)
⎩
G(fj ) cos(− 2 − πj) − iG(fj ) sin(− π2 − πj)
π
fN q < fj < 2fN q

Due to the truncation of the Fourier transformation, the filter frequency

response will differ from the desired frequency response. To illustrate the
effect of the least-square fit, both the gain and phase characteristic of a
linear FIR Hilbert filter are plotted in figure 8.3.

108
 1.5 2.0

1.0
1.0
0.5

Phase
Gain

0.0 0.0
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0
-0.5
-1.0
-1.0

-1.5 -2.0
frequency (Hz) frequency (Hz)

Figure 8.3: Gain and phase characteristic of linear FIR Hilbert filter with a
filter length Nc = 64 and fs = 1.0 Hz.

To compare the FIR approximated Hilbert transform with the theoretical

Hilbert transform an irregular time signal is generated from the JONSWAP
spectrum and the two transforms are depicted in figure 8.4. Generally very
good accordance is observed also at the edges where a zone of half the filter
length normally is disturbed.

Original signal
Theoretical Hilbert transform
10.0 Approximated Hilbert transform

7.5

5.0
Signal

2.5

0.0
0 50 100 150 200
-2.5

-5.0

time (sec.)

Figure 8.4: Comparison of theoretical and FIR approximated, Nc = 64,

Hilbert transform. The signal is generated from the JONSWAP spectrum,
fp = 0.1 Hz and γ = 3.3.

To illustrate the envelope function, E(t) is plotted together with half the
squared elevation in figure 8.5 for a time signal generated from the JONSWAP
spectrum.

109
 Original signal
Half the square of signal
10.0 Half the square envelope

7.5

5.0
Signal

2.5

0.0
0 50 100 150 200
-2.5

-5.0

time (sec.)

Figure 8.5: Comparison of half the square envelope E(t) and 12 η 2 (t) for signal
generated from a JONSWAP spectrum, fp = 0.1 Hz, γ = 3.3, and Nc = 64.

8.3 Groupiness Factor

To characterize the actual groupiness of a wave train the energy spectrum

Sη̂ (f ) of half the square envelope function can be evaluated. However, a sim-
pler measure is the groupiness factor that is defined as the standard deviation
of half the square envelope relative to the variance of the original time signal

σ[E(t)]
GF = (8.21)
σ 2 [η(t)]

For a monochromatic (sinusoidal) signal the envelope function E(t) is con-

stant leading to a groupiness factor GF = 0. Taking a completely Gaussian
signal the expected value of the groupiness factor can be shown to be equal
to 1.0 independent of the spectrum shape. The actual values for time signals
generated from a JONSWAP spectrum including approximately 500 periods
are approximately 1.0 in mean with a standard deviation of approximately
σ = 0.13.

Instead of computing one value of the groupiness factor over the complete
length of the time signal, the groupiness factor can be evaluated as instanta-
neous values by computing an average groupiness factor over a time moving
window. The length of the window in time is dependent on the desired degree
of smoothing of the computed groupiness factor function.

In figure 8.6 to figure 8.9 the groupiness factor function is plotted for both
a narrow-banded and a broad-banded JONSWAP spectrum for two different
window sizes.

110
 Original signal
Half the square envelope
10.0 Groupiness factor

7.5

5.0
Signal

2.5

0.0
0 100 200 300 400
-2.5

-5.0

time (sec.)

Figure 8.6: Groupiness factor function GF(t) for signal generated from JON-
SWAP spectrum, fp = 0.1 Hz, fs = 1.0 Hz, γ = 10.0, Nc = 64, and window
size = Tm .

Original signal
Half the square envelope
10.0 Groupiness factor

7.5

5.0
Signal

2.5

0.0
0 100 200 300 400
-2.5

-5.0

time (sec.)

Figure 8.7: Groupiness factor function GF(t) for signal generated from JON-
SWAP spectrum, fp = 0.1 Hz, fs = 1.0 Hz, γ = 10.0, Nc = 64, and window
size = 3 Tm .

Generally, a more smooth groupiness factor function is obtained for a window

size of 3 mean periods and only the largest wave groups are separated as high
and smooth peaks. It should though be noted that the sample frequency is
1.0 Hz and that a higher sample frequency eventually will lead to smoother
groupiness factor function for smaller window sizes.

111
 Original signal
Half the square envelope
10.0 Groupiness factor

7.5

5.0
Signal

2.5

0.0
0 100 200 300 400
-2.5

-5.0

time (sec.)

Figure 8.8: Groupiness factor function GF(t) for signal generated from JON-
SWAP spectrum, fp = 0.1 Hz, fs = 1.0 Hz, γ = 1.0, Nc = 64, and window
size = Tm .

Original signal
Half the square envelope
10.0 Groupiness factor

7.5

5.0
Signal

2.5

0.0
0 100 200 300 400
-2.5

-5.0

time (sec.)

Figure 8.9: Groupiness factor function GF(t) for signal generated from JON-
SWAP spectrum, fp = 0.1 Hz, fs = 1.0 Hz, γ = 1.0, Nc = 64, and window
size = 3 Tm .

8.4 Conclusions and Further Use

Based on a linear assumption a method for calculating the instantaneous

wave energy history and the groupiness factor function has been presented.
The method is based on a temporal Hilbert filter and this approach enables
an exact isolation of the 2nd order subharmonics which describe the slowly
varying part of the time signal. This Hilbert filter approach is thus more
efficient than the SIWEH analysis. The groupiness factor has proven to be
ineffective in describing Gaussian distributed sea surfaces and the groupiness
factor function is defined. Also discussions regarding the implementation of
the Hilbert filter using FIR approximations and choice of window sizes for
computing the groupiness factor functions are made.

112
The method can easily be extended to a three-dimensional motion but a
physical interpretation of the more slowly varying part must then be revised.

The groupiness factor function enables computations of instantaneous groupi-

ness factors in time and hence, the function is suitable for comparing the
correlation between the damage development of e.g. a breakwater and the
wave grouping in the wave train causing the damage.

A further application is the possibility to evaluate the change in wave group-

ing due to shoaling and thus also the change in phase distribution from deep
to shallow water.

113
Chapter 9

References

Aalborg University, 2012.

WaveLab 3 homepage. http://www.hydrosoft.civil.aau.dk/wavelab/

Battjes, J. A. and Groenendijk, H. W.

Wave height distributions on shallow foreshores.
Coastal Engineering, Vol. 40, 2000, pp. 161-182.

Brorsen, M. and Frigaard, P., 1992.

Active Absorption of Irregular Gravity Waves in BEM-models. Boundary El-
ements XIV, Vol.1 Editors: Brebbia, Dominiquez and Paris. Computational
Mechanics Publications, Southampton.

Christensen, M. and N. B. Sørensen, 1994.

”Maximum Likelihood Estimation of Directional Spectrum Expressed in Stan-
dard Form”
Aalborg University.

Davis, Russ E. and Regier, Lloyd A., 1977.

”Methods for estimating directional wave spectra from multi-element arrays”.
Journal of Marine Research.

Frigaard, P. and Brorsen, M., 1993.

”A time-domain method for Separating Incident and Reflected Irregular
Waves.
Coastal Engineering, Vol. 24, Nos. 3 and 4.

115
Frigaard, P., Helm–Petersen, J., Klopman, G., Stansberg, C.T., Benoit, M.,
Briggs, J., Miles, M., Santas, J., Schäffer, H.A. and Hawkes, P.J., 1997.
IAHR List of Sea State Parameters – an update for multidirectional waves.
IAHR Seminar Multidirectional Waves and their Interaction with Structures,
IAHR Congress, San Francisco.

Grønbech, J., Jensen, T. and Andersen, H. (1996).

Reflection Analysis with Separation of Cross Modes.
Proc. 25th Int. Conf. on Coastal Eng., Orlando, Florida, USA, pp. 968-980.

Goda, Y. and Suzuki, Y., 1976.

”Estimation of Incident and Reflected Waves in Random Wave Experiments”
Proceedings, 15th International Conference on Coastal Engineering, Vol. 1,
pp. 828-845, Honolulu, Hawaii.

Goda, Y., 1986.

Effect of wave tilting on zero-crossing wave heights and periods.
Coastal Engineering in Japan, Vol. 29, pp. 79-90.

Goda, Y., 2010.

Random Seas and Design of Maritime Structures, 3rd Edition.
Advanced Series on Ocean Eng., Vol. 38, World scientific publishing.

Hashimoto, N. and Kobune, K., 1988.

”Directional spectrum estimation from a Bayesian approach”.
Proceedings, 21st Conference on Coastal Engineering, Vol. 1, pp 62-76, Costa
del Sol-Malaga, Spain.

Helm–Petersen, J., 1993.

”The Bayesian Directional Spectrum Estimation Method”
Master Thesis, Aalborg University

Isobe, M., 1990.

116
”Estimation of Directional Spectrum Expressed in Standard Form”
Proc. 22nd Int. Conf on Coastal Eng., pp. 647-660.

Isobe, M. and K. Kondo, 1984.

”Method for Estimating Directional Wave Spectrum in Incident and Re-
flected Wave Field”
Proc. 19th Int. Conf on Coastal Eng., pp. 467-483.

Isobe, M., Kondo, K. and Horikawa, K., 1984.

”Extension of MLM for estimating directional wave spectrum”.
University of Tokyo, Japan.

Karl, John H., 1989.

”An Introduction to Digital Signal Processing”
Academic Press, San Diego.

Klopmann, G. and Stive, M.J.F., 1989.

”Extreme waves and wave loading in shallow water.”
E & P Forum, Report No. 3.12/156.

Mansard, E. and Funke, E., 1980.

”The Measurement of Incident and Reflected Spectra Using a Least Squares
Method”
Proceedings, 17th International Conference on Coastal Engineering, Vol. 1,
pp 154-172, Sydney, Australia.

Nelder, J. A. and R. Mead, 1965.

Computer Journal, vol. 7, p. 308.

Neu, H.J.A., 1982.

”11-year deep water wave climate of Canadian Atlantic waters.”
Canadian Tech. Rept. of Hydrography and Ocean Sciences, 13.

Newland, D.E., 1975.

”In introduction to random vibrations and spectral analysis.”
Longman, London, 1975.

117
Ochi, Michel K. 1998
”Ocean Waves - The Stochastic Approach”
Cambridge University Press, Cambridge. ISBN: 0-521-56378

Press, W. H., B. P. Flannery, S. A. Teukolsky and W. T. Vetterling, 1989.

”Numerical Recipes in Pascal”.
Cambridge University Press, Cambridge.

Raunholt, L. & Skjærbæk, P. 1990

Optimering af lodrette perforerede plader som bølgeabsorber (in Danish)
Master Thesis, Aalborg University

Stive, M.J.F., 1986.

”Extreme shallow water conditions.”
Delft Hydraulics, Intern Report H533.

Tucker, M.J., 1991.

”Waves in ocean engineering: measurement, analysis, interpretation”.
Ellis Horwood Ltd.

Yokoki, H., M. Isobe and A. Watanabe, 1992.

”A Method for Estimating Reflection Coefficient in Short Crested Random
Seas”
Proc. 23rd Int. Conf on Coastal Eng., pp. 765-776.

Young, I.R., 1999

”Wind Generated Ocean Waves.”
Elsevier, Kidlington, Oxford 1999. ISBN: 0-08-04331

Zelt, J.A. and Skjelbreia, J., 1992.

”Estimating Incident and Reflected Wave Fields Using an Arbitrary Number
of Wave Gauges”
Proceedings, 23th International Conference on Coastal Engineering, Vol. 1,
pp 777-789, Venice, Italy.

118
ISSN 1901-7286
DCE Lecture Notes No. 33