Spectral Analysis and Time Series

Andreas Lagg

Part I: fundamentals on time series

classification

prob. density func.

auto-correlation

power spectral density

cross-correlation

applications

pre-processing

sampling

trend removal

Part II: Fourier series

definition

method

properties

convolution

correlations

leakage / windowing

irregular grid

noise removal

Part III: Wavelets

why wavelet transforms?

fundamentals:
FT, STFT and resolution problems

multiresolution analysis: CWT

DWT

Exercises

Basic description of physical data

deterministic: described by explicit mathematical relation

non deterministic

non deterministic: no way to predict an exact value at a future instant of time

Classifications of deterministic data

Periodic

Deterministic

Transient

Almost periodic

Complex Periodic

Sinusoidal

Nonperiodic

Sinusoidal data

time history

frequency spectrogram

Complex periodic data

time history

frequency spectrogram

(T = fundamental period)

Almost periodic data

time history

frequency spectrogram

no highest common divisor -> infinitely long period T

Transient non-periodic data

all non-periodic data other than almost periodic data

Classification of random data

Stationary

Random Data

Special classifications

of nonstationarity

Nonergodic

Ergodic

Nonstationary

stationary / non stationary

collection of sample functions = ensemble

autocorrelation function (joint moment):

mean value (first moment):

stationary:

data can be (hypothetically) described by computing ensemble averages (averaging over multiple measurements / sample functions)

weakly stationary:

ergodic / non ergodic

Ergodic random process:
properties of a stationary random process described by computing averages over only one single sample function in the ensemble

autocorrelation function (joint moment):

mean value of k-th sample function:

ergodic:

Basic descriptive properties of random data

mean square values

propability density function

autocorrelation functions

power spectral density functions

(from now on: assume random data to be stationary and ergodic)

Mean square values

(mean values and variances)

describes general intensity of random data:

rout mean square value:

often convenient:

static component described by mean value:

dynamic component described by variance:

standard deviation:

Probability density functions

describes the probability that the data will assume a value within some defined range at any instant of time

for small :

probability density function

probability distribution function

Illustration: probability density function

sample time histories:

sine wave (a)

sine wave + random noise

narrow-band random noise

wide-band random noise

all 4 cases: mean value

Illustration: probability density function

probability density function

Autocorrelation functions

describes the general dependence of the data values at one time on the values at another time.

(not for special cases like sine waves)

Illustration: ACF

autocorrelation functions (autocorrelogram)

Illustrations

autocorrelation function

autocorrelation function of a rectangular pulse

Power spectral density functions

(also called autospectral density functions)

describe the general frequency composition of the data in terms of the spectral density of its mean square value

mean square value in frequency range :

portion of x(t) in (f,f+∆f)

definition of power spectral density function:

important property: spectral density function is related to the autocorrelation function by a Fourier transform:

Illustration: PSD

sine wave

broadband
noise

narrowband
noise

sine wave
+ random noise

power spectral density functions

Dirac delta function at f=f0

“white” noise:
spectrum is uniform over all frequencies

Joint properties of random data

Joint probability density functions

joint properties in the amplitude domain

Cross-spectral density functions

joint properties in the frequency domain

Cross-correlation functions

joint properties in the time domain

joint probability measurement

until now: described properties of an individual random process

Cross-correlation function

cross-correlation measurement

describes the general dependence of one data set to another

similar to autocorrelation function

functions are uncorrelated

typical cross-correlation plot (cross-correlogram):
sharp peaks indicate the existence of a correlation between x(t) and y(t) for specific time displacements

Applications

Measurement of time delays

2 signals:

different offset

different S/N

time delay 4s

often used:
'discrete' cross correlation coefficient
lag =l, for l >=0:

Applications

Detection and recovery from signals in noise

3 signals:

noise free replica of the signal (e.g. model)

2 noisy signals

cross correlation can be used to determine if theoretical signal is present in data

Pre-processing Operations

sampling considerations

trend removal

filtering methods

sampling

cutoff frequency (=Nyquist frequency or folding frequency)

Trend removal

often desirable before performing a spectral analysis

Least-square method:

time series:

desired fit
(e.g. polynomial):

Lsq-Fit: minimize

set partial
derivatives to 0:

K+1 equations:

Digital filtering

end of part I ...

Spectral Analysis and Time Series

Andreas Lagg

Part I: fundamentals on time series

classification

prob. density func.

auto-correlation

power spectral density

cross-correlation

applications

pre-processing

sampling

trend removal

Part II: Fourier series

definition

method

properties

convolution

correlations

leakage / windowing

irregular grid

noise removal

Part III: Wavelets

why wavelet transforms?

fundamentals:
FT, STFT and resolution problems

multiresolution analysis: CWT

DWT

Exercises

Fourier Series and Fast Fourier Transforms

Standard Fourier series procedure:

if a transformed sample record x(t) is periodic with a period Tp (fundamental frequency f1=1/Tp), then x(t) can be represented by the Fourier series:

where

Fourier series procedure - method

sample record of finite length, equally spaced sampled:

Fourier series passing through these N data values:

Fill in particular points:

coefficients Aq and Bq:

inefficient & slow => Fast Fourier Trafos developed

Fourier Transforms - Properties

Linearity

Symmetry

Circular time shift

Using FFT for Convolution

Convolution Theorem:

Fourier transform of the convolution is product of the individual Fourier transforms

(note how the response function for negative times is wrapped around and stored at the extreme right end of the array)

discrete case:

Convolution Theorem:

constraints:

duration of r and s are not the same

signal is not periodic

original data

response function

convolved data

Treatment of end effects by zero padding

constraint 1:

simply expand response function to length N by padding it with zeros

constraint 2:

extend data at one end with a number of zeros equal to the max. positive / negative duration of r (whichever is larger)

FFT for Convolution

zero-pad data

zero-pad response function
(-> data and response function have N elements)

calculate FFT of data and response function

multiply FFT of data with FFT of response function

calculate inverse FFT for this product

Deconvolution

-> undo smearing caused by a response function

use steps (1-3), and then:

divide FFT of convolved data with FFT of response function

calculate inverse FFT for this product

Correlation / Autocorrelation with FFT

definition of correlation / autocorrelation see first lecture

Correlation Theorem:

Auto-Correlation:

discrete correlation theorem:

Fourier Transform - problems

spectral leakage

reducing leakage by windowing (1)

Applying windowing (apodizing) function to data record:

(original data record x windowing function)

No Window:

Bartlett Window:

Hamming Window:

Blackman Window:

Hann Window:

Gaussian Window:

reducing leakage by windowing (2)

without windowing

with Gaussian windowing

No constant sampling frequency

Fourier transformation requires constant sampling (data points at equal distances)
-> not the case for most physical data

Solution: Interpolation

linear:
inear interpolation between yk and yk+1
IDL> idata=interpol(data,t,t_reg)

quadratic:
quadratic interpolation using yk-1, yk and yk+1
IDL> idata=interpol(data,t,t_reg,/quadratic)

least-square quadratic
least-square quadratic fit using yk-1, yk, yk+1 and yk+2
IDL> idata=interpol(data,t,t_reg,/lsq)

spline
IDL> idata=interpol(data,t,t_reg,/spline)
IDL> idata=spline(t,data,t_reg[,tension])

important:
interpolation changes sampling rate!
-> careful choice of new (regular) time grid necessary!

Fourier Transorm on irregular gridded data - Interpolation

original data: sine wave + noise

FT of original data

irregular sampling of data (measurement)

interpolation: linear, lsq, spline, quadratic

're-sampling'

FT of interpolated data

Noise removal

Frequency threshold (lowpass)

make FT of data

set high frequencies to 0

transorm back to time domain

Noise removal

signal threshold for
weak frequencies (dB-threshold)

make FT of data

set frequencies with amplitudes below a given threshold to 0

transorm back to time domain

Optimal Filtering with FFT

normal situation with measured data:

underlying, uncorrupted signal u(t)
+	response function of measurement r(t)
= smeared signal s(t)

+	noise n(t)
=	smeared, noisy signal c(t)

estimate true signal u(t) with:

= optimal filter
(Wiener filter)

Calculation of optimal filter

reconstructed signal and uncorrupted signal should be close in least-square sense:
-> minimize

additional information:
power spectral density can often be used to disentangle noise function N(f) from smeared signal S(f)

Using FFT for Power Spectrum Estimation

Periodogram

discrete Fourier transform of c(t)
-> Fourier coefficients:

-> periodogram estimate of power spectrum:

end of FT

Spectral Analysis and Time Series

Andreas Lagg

Part I: fundamentals on time series

classification

prob. density func.

auto-correlation

power spectral density

cross-correlation

applications

pre-processing

sampling

trend removal

Part II: Fourier series

definition

method

properties

convolution

correlations

leakage / windowing

irregular grid

noise removal

Part III: Wavelets

why wavelet transforms?

fundamentals:
FT, STFT and resolution problems

multiresolution analysis: CWT

DWT

Exercises

Introduction to Wavelets

why wavelet transforms?

fundamentals: FT, short term FT and resolution problems

multiresolution analysis:
continous wavelet transform

multiresolution analysis:
discrete wavelet transform

Fourier: lost time information

6 Hz, 4 Hz, 2 Hz, 1 Hz

6 Hz + 4 Hz + 2 Hz + 1 Hz

Solution: Short Time Fourier Transform

(STFT)

perform FT on 'windowed' function:

example: rectangular window

move window in small steps over data

perform FT for every time step

Short Time Fourier Transform

STFT

STFT-spectrogram shows both time and frequency information!

Gauss-functions as windows

Short Time Fourier Transform: Problem

narrow window function -> good time resolution
wide window function -> good frequency resolution

Solution: Wavelet Transformation

time vs. frequency resolution is intrinsic problem (Heisenberg Uncertainty Principle)
approach: analyze the signal at different frequencies with different resolutions

-> multiresolution analysis (MRA)

Continuous Wavelet Transform

similar to STFT:

signal is multiplied with a function (the wavelet)

transform is calculated separately for different segments of the time domain

but:

the FT of the windowed signals are not taken
(no negative frequencies)

The width of the window is changed as the transform is computed for every single spectral component

Continuous Wavelet Transform

mother wavelet:

finite length (compactly supported) ->'let'

oscillatory ->'wave'

functions for different regions are derived from this function -> 'mother'

scale parameter s replaces frequency in STFT

The Scale

similar to scales used in maps:

high scale = non detailed global view (of the signal)

low scale = detailed view

scaling dilates (stretches out) or compresses a signal:
s > 1 -> dilation
s < 1 -> compression

in practical applications:

low scales (= high frequencies) appear usually as short bursts or spikes

high scales (= low frequencies) last for entire signal

Computation of the CWT

signal to be analyzed: x(t), mother wavelet: Morlet or Mexican Hat

start with scale s=1 (lowest scale, highest frequency)
-> most compressed wavelet

shift wavelet in time from t0 to t1

increase s by small value

shift dilated wavelet from t0 to t1

repeat steps for all scales

CWT - Example

signal x(t)

axes of CWT: translation and scale (not time and frequency)

translation -> time
scale -> 1/frequency

Time and Frequency Resolution

every box corresponds to a value of the wavelet transform in the time frequency plane

all boxes have constant area
Δf Δt = const.

low frequencies: high resolution in f, low time resolution

high frequencies: good time resolution

STFT: time and frequency resolution is constant (all boxes are the same)

Wavelets: Mathematical Approach

Mexican Hat wavelet:

Morlet wavelet:

WL-transform:

inverse WL-transform:

admissibility condition:

Discretization of CWT: Wavelet Series

-> sampling the time – frequency (or scale) plane

advantage:

sampling high for high frequencies (low scales)
scale s1 and rate N1

sampling rate can be decreased for low frequencies (high scales)
scale s2 and rate N2

discrete wavelet

continuous wavelet

WL-transformation

reconstruction of signal

Discrete Wavelet Transform

discretized continuous wavelet transform is only a sampled version of the CWT

The discrete wavelet transform (DWT) has significant advantages for implementation in computers.

(DWT)

excellent tutorial:
http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html

IDL-Wavelet Tools:
IDL> wv_applet

Matlab-Wavelet Toolbox Documentation

Wavelet expert at MPS:
Rajat Thomas

end of Wavelets

Exercises

Part I: Fourier Analysis

(Andreas Lagg)

Part II: Wavelets

(Rajat Thomas)

Instructions:

http://www.linmpi.mpg.de/~lagg

Seminar room
Time: 14:00

Exercise: Galileo magnetic field

data set from Galileo magnetometer (synthesized)

file: gll_data.sav, contains:

total magnetic field

radial distance

time in seconds

Tips:

restore,'gll_data.sav'

use IDL-FFT

remember basic plasma physics formula for the ion cyclotron wave:

your tasks:

Which ions are present?

Is the time resolution of the magnetometer sufficient to detect electrons or protons?

Background:

If the density of ions is high enough they will excite ion cyclotron waves during gyration around the magnetic field lines. This gyration frequency only depends on mass per charge and on the magnitude of the magnetic field.

In a low-beta plasma the magnetic field dominates over plasma effects. The magnetic field shows only very little influence from the plasma and can be considered as a magnetic dipole.

http://www.sciencemag.org/cgi/content/full/274/5286/396

Literature

Random Data: Analysis and Measurement Procedures
Bendat and Piersol, Wiley Interscience, 1971

The Analysis of Time Series: An Introduction
Chris Chatfield, Chapman and Hall / CRC, 2004

Time Series Analysis and Its Applications
Shumway and Stoffer, Springer, 2000

The Wavelet Tutorial
Robi Polikar, 2001
http://users.rowan.edu/~polikar/WAVELETS/WTtutorial.html

Numerical Recipies in C
Cambridge University Press, 1988-1992
http://www.nr.com/

whole lecture: 1-D functions of time, space, frequency or wavelength

incomplete

length untested

most physical data: function f(t) or f(x) or f(lambda)...

special interest for physicists: find a mathematical law behind these functions

rigid body on a linear, massless spring

deterministic: satellite orbit, potential across a condenser during discharge through a resistor, temperature of water if heat is applied

non deterministic: height of waves in confused sea, acoustic pressure generated by air rushing through a pipe, electrical output of noise generator

no physical data are 'truly deterministic': unforeseen event in the future might influence phenomenon

no process is 'truly non-deterministic': mathematical formulas allow to describe almost every process

waveform exactly repeats itself at regular intervals

fundamental frequency f1

with a few exceptions expandable to Fourier series
(see x(t) = a0/2 ...)

with a_n, b_n Fourier coefficients

fundamental component may be absent: eg frequencies 60, 75, 100 Hz -> fundamental frequency is 5 Hz (highest common divisor)

occur in nature when two or more unrelated physical phenomena are mixed

eg vibration response of a multiple engine propeller airplane when engines are out of synchronization

discrete frequency spectrum possible, only difference to periodic data: components are not related by rational numbers.

example 1: could represent temperature of water after heating is turned off

example 2: damped oscillation

example 3: output of wave generator

difference to periodic data: discrete spectral representation is not possible!

however: continuous spectral representation can be obtained in most cases from a Fourier integral (complex number)

right hand side shows real part of Fourier integral

cannot be described by mathematical relationship – each observation of the phenomenon will be unique.

Any observation will only present one of many possible results.

“Sample record” = single time history over a finite interval representing random phenomenon

collection of all possible sample functions == random or stochastic process.

Properties of phenomenon can hypothetically be described at any instant of time by computing average values over the collection of sample functions describing the random process.

mean value (first moment): instantaneous value of each sample at time t1, summing values and divide by number of sample functions

autocorrelation function: taking ensemble average of the product of instantaneous values at two times, t1 and t1+tau.

non-stationary: mu_x and R_x vary with t1

stationary (weakly): mu_x is independent of t1 (mu_x(t)=mu_x) and Rx is dependent only on the displacement tau.

fully stationary: all moments and joint moments are time invariant

Previous slide: properties of random process described by computing ensemble averages at specific instants of time.

In most cases it is possible to describe the properties of a random process by computing time averages over specific sample functions in the ensemble.

if the random process {x(t)} is stationary and mu and R_x do not differ when computed over different sample functions, then the random process is ergodic.

The time-averaged mean value and the autocorrelation function are equal to the corresponding ensemble averaged values.

Only stationary processes can be ergodic.

Ergodic random processes are important: all properties can be determined by performing time averages over a single sample function.

Fortunately, random data representing stationary physical phenomena are generally ergodic.

-> properties of phenomenon can be measured from a single observed time history record.

From now on we make the assumption that our time series are ergodic, so that the properties of the data can be determined from time averages of individual sample records.

power spectral density function and autocorrelation function furnish similar information in the time domain and frequency domain respectively. For stationary data the power spectral density function does not supply new information over the autocorrelation function since the two are Fourier pairs.

well known relations, mentioned for completeness

mean square value = average of squared values of the time history

often desirable: physical data as a combination of static (=time invariant) component and a dynamic (=fluctuating) component.

static component: described by a mean value (=average of all values)

dynamic component: described by a variance (=mean square value about the mean value)

probability that x(t) assumes a value in the interval x and (x+ dx) is obtained by taking the ratio Tx / T, where Tx is the amount of time that x(t) falls into (x, x+dx) during an observation time T.

Exact probability description: T -> infinity

proability density function: describes the probability tha a function assumes a certain value

P(x) is (cumulative) probability distribution function

P(x) is bounded by 0 and 1 (probability of x(t) being less than -inf is 0 and probability of x(t) being greater than inf is 1)

mean value and mean square value can be expressed in terms of probability density function p(x)

note: sine wave is deterministic

b) noise in amplitude

c+d) noise in frequency

probability density plot of random data is ideally Gauss shaped

Application:

the principal application is to establish probabilistic descriptions for instantaneous values of the data.

Often: non-linear physical processes exhibit their presence in the probability density function.

Estimate of autocorrelation function between values of x(t) at times t and t+tau can be obtained by taking the product of the two values x(t) and x(t+tau) and averaging over over observation time T.

formula: resulting average product will approach an exact autocorrelation function as T approaches infinity.

R_x(tau) is always real valued, even and has a maximum at t=0

mean value is given by square root of R_x(inf)

mean square value is R_x(0)

sine function: maximum correlation is when tau is equal to the period of the sine wave. phase information is lost.

sine wave + random noise: the underlying sine wave determines the maxima for the autocorrelation function. It is the sum of the autocorrelogram of a sine wave and of random noise.

White random noise (=random data with energy distributed uniformly over all frequencies) would show up as a delta peak at tau=0 in an autocorrelogram

narrow band noise: autocorrelogram diminishes to 0 (in the case of a non-zero mean-value mu_x it will approach mu_x)

Application: establish the influence of values at any time over values at a future time.

Tool for detecting deterministic data which might be masked in a random background.

better tool: Fourier transform

power spectral density function is also called auto spectral density function

Applications:

measurement of physical data to establish the frequency composition of the data.

bears important relationships to the basic characteristics of the physical system involved.

until now: we described the properties of an individual random process.

Often desirable: describe certain common or joint properties from two or more random processes.

the 3 functions are effective extensions and describe the joint properties in the amplitude domain, time domain and frequency domain respectively.

Joint probability density function: only mention here.

describes the probability that both sample records will simultaneously assume values within some defined pair of ranges at any instant of time.

similar definition to the probability density function and the probability distribution function, but in 2-D

signals: eg. crossing of a boundary layer by two different satellites.

'discrete' cross correlation: directly applicable to data sets.

no complicated signal fitting is required to obtain exact time delay between two data sets.

signal shape known or anticipated.

every digital measurement requires sampling.

sampling interval has to be small enough to recover signal shape

if sampling frequency is too small -> data are not reproduced correctly

oversampling: sampling frequency higher than the width of the measurement pulse.

Easy expressions for K=0 (remove constant background) and K=1 (remove linear trend)

important for spectral analysis (e.g. with Fourier): trend in data can be interpreted as very low frequency signal. Will destroy the power spectrum of the signal

highpass / lowpass

easy to realize with Fourier transform

Fourier series and Fourier transforms of data differ in their theoretical properties, but not, for most practical purposes, in their digital computational details.

Reason: only a finite range Fourier series or transform can actually be computed with digitized data, and this finite range can always be considered as the period of an associated Fourier series.

Assume sample record x(t) is of finite length T_p, the fundamental period of the data.

Further assume that record is sampled at an even number of N equally spaced points a distance h apart (h>cutoff frequency)

proceed now to calculate Fourier series which passes through these N data values. For any point in (0,T_p) the result is x(t) = ...

This represents an easy method to program your own FT. But it is very slow compared to more sophisticated computations like FFT.

response function with finite duration N

special response function: r0=1 and all other rk's are zero. -> identity filter, convolution gives identically the original signal

usually: all measurements have response functions:

- filter width when recording a spectrum

- ...

end effects: repeat response function at end of signal

deconvolution almost impossible to do in time domain

FT makes it trivial in the frequency domain

often apparently spurious components will appear in output of fft computation.

Spectral leakage occurs when frequency components in the input function f(t) are not harmonics of the fundamental frequency k0 = 2pi/(N dT), defined by the length of data record N dT. Under such conditions the lines in the DFT do not accurately reflect the amplitude of the component, and spurious components appear adjacent to the component.

The reason for the appearance of leakage components in a DFT is the convolution of the data function and the truncation window. Each sinusoidal component in f(t) has a spectrum F(j omega) that is a pair of impulses in the frequency domain: multiplication by the truncating function causes a spread in the width of the component.

Leakage may be reduced (but not eliminated) by using another function to truncate the series instead of the implicit rectangular function. Such a function is known as an apodizing, or windowing, function and is chosen to taper the data record to zero at the extremities while minimizing the spectral spreading of each component.

we have: C(f)

we know how to deconvolve the response function: divide C(f) by R(f)

now: search filter function PHI(f) which, when applied to C(f)/R(f) gives signal close to U(f)

underlying noise spectrum may be flat, tilted, smoothly varying. This is not important, as long as we can guess a reasonable hypothesis as to what it is.

left time series: 4 frequencies in for different time intervals

right: all 4 frequencies over the whole time interval

FFT right: peak at all 4 frequencies contained in the spectrum. Noise due to spectral leakage

ALL FREQUENCIES AT ALL TIMES

FFT left: peak at all 4 frequencies – very similar to right figure

the FT gives what frequency components (spectral components) exist in the signal. Nothing more, nothing less.

FT: no resolution problem in frequency domain (we know exactly which frequencies exist)!

time domain: no time resolution problem, we know the value of the signal at any instant of time

-> time resolution in FT is 0

-> frequency resolution in time domain is 0

Perfect frequency resolution in FT is the fact that the window used in the FT (exp{jwt} function) lasts at all times from minus infinity to plus infinity.

Now, in STFT, our window is of finite length, thus it covers only a portion of the signal, which causes the frequency resolution to get poorer.

We no longer know the exact frequency components that exist in the signal, but we only know a band of frequencies that exist.

If we use a window of infinite length, we get the FT, which gives perfect frequency resolution, but no time information.

Narrow window ===>good time resolution, poor frequency resolution.

Wide window ===>good frequency resolution, poor time resolution.

MRA is designed to give good time resolution and poor frequency resolution at high frequencies and good frequency resolution and poor time resolution at low frequencies.

-> Useful when the signal at hand has high frequency components for short durations and low frequency components for long durations. Fortunately, the signals that are encountered in practical applications are often of this type.

CWT was developed as an alternative approach to the STFT to overcome the resolution problem.

The wavelet analysis is done in a similar way to the STFT analysis, in the sense that the signal is multiplied with a function, {the wavelet}, similar to the window function in the STFT, and the transform is computed separately for different segments of the time-domain signal. Difference:

1. The Fourier transforms of the windowed signals are not taken, and therefore single peak will be seen corresponding to a sinusoid, i.e., negative frequencies are not computed.

2. The width of the window is changed as the transform is computed for every single spectral component, which is probably the most significant characteristic of the wavelet transform.

The parameter scale in the wavelet analysis is similar to the scale used in maps. As in the case of maps, high scales correspond to a non-detailed global view (of the signal), and low scales correspond to a detailed view. Similarly, in terms of frequency, low frequencies (high scales) correspond to a global information of a signal (that usually spans the entire signal), whereas high frequencies (low scales) correspond to a detailed information of a hidden pattern in the signal (that usually lasts a relatively short time).

Cosine signals corresponding to various scales are given as examples in the following figure .

Unlike the STFT which has a constant resolution at all times and frequencies, the WT has a good time and poor frequency resolution at high frequencies, and good frequency and poor time resolution at low frequencies.

Mexican hat: second derivative of Gaussian function

which wavelet to use: experience

admissibility constant: if this condition is fulfilled, a successful reconstruction is possible.

2002-05-09T12:58:24