Introduction to Sampling Theory and Data Analysis (Part 1)

Principles of Oceanographic Instrument Systems: Sensors and Measurements Spring 2004 INTRODUCTION TO SAMPLING THEORY AND DATA ANALYSIS These notes are meant to introduce the ocean scientist and engineer to the concepts associated with the sampling and analysis of oceanographic time series data, and the effects that the sensor, recorder, sampling plan and analysis can have on the results. In order to plan the optimum sampling and analysis plan, one needs to understand what information and analysis are required, and how all these factors will affect the final result. To get the most from these lecture notes, the student should do supplemental readings from the references listed below. Exercises utilizing the MATLAB software package will be assigned at the appropriate place in the lectures. An outline of this section is given below, and covered in handouts. 1. Time Series and Analysis: •Properties of a random, stochastic processes •Statistical description: mean, variance, correlation/covariance, spectra •Fourier transforms, frequency domain/time domain description of a process •Digital filtering and filters: Convolution product, filters, and filter response 2. Sampling Theory: •Sampling process, sampling theorem, and sampling effects on statistics •Aliasing and the Nyquist frequency •Power density spectra, coherence, degrees of freedom, confidence limits 3. Environmental Sampling in the real world: •Calibrations: static, dynamic •Digitizing effects, prefiltering •Sensor frequency response effects •Sensor noise limitations Suggested references and readings: Jenkins, G.M. and D.G. Watts, Spectral Analysis and its Applications, Holden-Day, San Francisco, 1968. Koopmans, L.H., The Spectral Analysis of Time Series, Academic Press, New York, 1974. Bendat, J.S. and A.G. Piersol, Random Data: Analysis and Measurement Procedures, Wileyinterscience, New York, Second Edition, 1986. Daley, R., Atmospheric Data Analysis, Cambridge University Press, New York, 1991. Cochran, W.T., et al, “What is the Fast Fourier Transform,” IEEE Trans. Audio and Electroacoustics, AU-15(2), 45-55, 1967. 1 Glossary of Terms, from Blackman, R.B. and J.W. Tukey, The Measurement of Power Spectra, Dover, 1958. Bingham, C., M.D. Godfrey and J.W. Tukey, “Modern Techniques of Power Spectrum Estimation,” IEEE Trans. Audio and Electroacoustics, AU-15(2), 56-66, 1967. Welch, P.D., “The use of Fast Fourier Transform for the Estimation of Power Spectra: A Method Based on Time Averaging of Short, Modified Periodograms,” IEEE Trans. Audio and Electroacoustics, AU-15(2), 70-73, 1967. Carter G.C., C.H. Knapp and A.H. Nutall, “Estimation of the Magnitude-Squared Coherence Function Via Overlapped Fast Fourier Transform Processing,” IEEE Trans. Audio and Electroacoustics, AU-21(4), 337-344, 1967. Background Everyone has some idea vague of what is involved in making measurements of the environment. However, few people have the background to really know how to do it properly. This is an introduction to how to measure the environment and analyze the results to obtain information for scientific studies and management decisions. As an ocean scientist or engineer you desire to make and analyze observations that will give you certain statistics describing the environment. To get these statistics, you need to design an experiment, place sensors in the field, digitize and record the results, analyze them on a computer, and finally present them in a meaningful manner. In reality, all these processes that you must go through can be thought of as a filter, or that you are looking at the ocean through “colored” glasses. In order to know what your glasses are doing to your view of the ocean, you need to know how to design an experiment to get the data that you want, select the sensors which will properly measure the environment, use recorders that will satisfactorily record the data, and utilize analysis techniques which will give the desired results. What follows is a simple introduction to the background that you will need to know in order to sample the environment properly. To simplify the discussions, much of the statistical complexity has been removed, so in order to become really professionally involved in data analysis, further course work is required to fill in this statistical information. Properties of Random variables We make the assumption that the environmental data of interest is a stationary, random, stochastic process. If this is so, then the environmental process that we wish to study can be fully described by its statistics. Random Variables - A deterministic variable is one whose value may be determined or estimated exactly. An example of a variable which can be predicted is the result from an explicit mathematical relationship, e.g. y(x) = a + bx or y(x,t) = cos(kx - ωt + θ). A random variable is one in which perfect prediction of succeeding values is impossible. Examples of random variables are the time until the next alpha particles is emitted from a radioactive source, the next direction taken by a particle in Brownian motion, or the elevation of the sea surface at a specific latitude, longitude and time. A set of observations of a random variable represents only one of many possible realizations. 2 Stochastic Process - A stochastic process is a collection of random variables. One observes a stochastic process when he examines a process developing in time in a manner controlled by probabilistic laws. A single set of observations is called a "sample function" or “sample record.” A random stochastic process is described by all its possible sample functions. Repeated observations will result in sample functions that are different, or are not the same function of time, but have the same statistics. Some examples of stochastic processes are the number of particles emitted from a radioactive source, the path of a particle in Brownian motion, or the sea surface elevation variations due to surface wind waves. One can not predict exactly any succeeding values, but one can describe succeeding values statistically. Stationary processes - The assumption that a random process is stationary is the most important assumption made in time series analysis. Perhaps this assumption is bad, at best it is only approximately true. A process is stationary when its "statistics" remain constant with time. Examples would be the output of a white noise generator, or the path of a particle in Brownian motion. An economic time series of Gross National Product (GNP) of the U.S. tends to increase with time, so is non-stationary. Hence, a stationary process is in statistical equilibrium and contains no trends or ramps. In actual practice we see three kinds of processes, 1) stationary, such as the output of the white noise generator, 2) quasi- stationary over a short period, such as atmospheric turbulence over a few minutes, or ocean waves over a few tens of minutes, and 3) non-stationary, such as the GNP, where the properties (statistics) are obviously changing with time. An oceanic example of a non-stationarity in a time series would be the surface wave field at the WHOI dock measured for 3 hours. During that time the tide would change the mean sea level by an amount which would "look" like a trend, but is really just a large, low frequency signal which is not resolved by our short record length. The wave field also might be growing in response to wind forcing, so the wave height and wavelength are changing with time. Another example would be the short-term temperature fluctuations observed for a month during the spring which are superimposed on the yearly warming and cooling. Since most geophysical spectra are “red,” or have more energy at low frequencies than at high frequencies, this may be a significant problem. Therefore, “beware” of unresolved low frequencies. Time - Frequency Domain Description of a Process Time Series - Everyone is familiar with a time series representation of some phenomena, that is a series of values at successive points in time. For example the tides or surface waves seen at the coast appear to vary somewhat like a sine wave as a function of time. Observations of oceanographic temperature or currents at a particular position as a function of time is a time series. Time is a continuous variable, and is represented as a function of time by x(t) as shown in Figure 1 below. Although geophysical processes are continuous in time, practical considerations require that we sample environmental process at discrete intervals in time, δt, for a finite length of time, T. That is, at regular intervals such as once a second, one minute, two hours, etc., an observation is made of the continuous process for a finite length of points. In order to simplify the analysis, sampling should always be done at equally spaced intervals in time. (It should be noted that the monthly averages often tabulated are not evenly spaced in time, so should not be analyzed as such.) We will consider a process, x, which is sampled at equally spaced time intervals, δt. 3 Time Series Example 1.28 1.26 Geophysical Value - x 1.24 1.22 X(t) 1.2 1.18 1.16 1.14 1.12 1.1 0 20 40 60 80 100 Time Figure 1. X is a continuous function of time, and is plotted with its value on the ordinate axis with value increasing upward and time increasing to the right on the abscissa. If the number of samples or terms in the time series is n, then the length of the series (which was T) will be nδt. The relationship between the discrete sample set, xt, (what we have to work with) and the original continuous function, x(t), will be covered below, but it is obvious that one needs to sample fast enough to be able to “see” the higher frequency fluctuations in the signal. Discrete Series 1.28 1.26 Geophysical Value - x 1.24 1.22 1.2 ->|dt|<1.18 1.16 1.14 1.12 1.1 |<------------------------------T---------------------------------->| 0 20 40 60 80 100 Time Figure 2. Discrete time series, xt, of observations shown in Figure 1, but sampled every 4 units in time, or δ t = 4. 4 Statistical descriptions - The second assumption (after stationarity) that we make is that the process may be adequately described by the lower moments of its probability density distributions. These include the mean, the variance, and the covariance, with its transform, the power spectrum. In our discussions below, we shall make a simplification by eliminating the probability density distribution function from our expected values. This would show up as multiplying the probability that x has a certain value in the integrals below. This makes the concepts involved in time series analysis a bit easier to follow in this introduction by eliminating some of the statistical considerations. However, in order to fully understand and use the power involved in proper data analysis procedures, one must go back and cover the concepts discussed below with a full statistical approach. Mean. The mean is just the average or expected value for a process. Consider x(t) as a continuous realization of a sample function of a stochastic processes. Then the mean is just the "expected value" or most probable value, ⌠∞ ⎮ x(t)dt ⌡−∞ 1 µ = ----------- = lim --⌠∞ T→∞ 2T ⎮ dt ⌡ −∞ ⌠T ⎮ x(t)dt ⌡ -T Eq 1 This is just the "static" component of the process. In practice we take the discrete time series, xt, as the sample series which extends from t = 1 to m, at internals of time, δt. We can look at our statistics as being estimates of the actual statistics of the continuous process. As the number of terms in the series becomes large, then the estimated statistics converge to the actual statistics. The error in estimating the statistics is made up of several parts: 1. the sampling accuracy related to the length of the series, m, and sample interval, δt - that is how closely the sampled series, xt, represents the true series, x(t), also of importance is the 2. digitizing accuracy used to create the actual numbers, xt, and finally the 3. computer accuracy in doing the actual analysis calculations. On modern computers (workstations and PC's) this is generally not a problem. The discrete representation of the mean is just the sum over all the terms, normalized by the number of terms m m µ ≈ 1/m ∑ xi or µ ≈ µ’ = 1/m ∑ xi i=1 i=1 Eq 2 as m, the number of terms, goes to infinity, this converges to the real mean. Mean Square and Root Mean Square. The "intensity" of the series is given by the mean square value 5 ⌠ ∞ ⎮ x(t)²dt ⌡ −∞ 1 Ψ² = --------- = lim --⌠∞ T→∞ 2T ⎮ dt ⌡ −∞ ⌠ T ⎮ x(t)²dt ⌡ -T Eq 3 Its positive square root is called the rms or "root mean square" value. Variance. The "dynamic component" of the series is given by the variance. This is the mean of the square of the differences from the mean, ⌠∞ ⎮ [x(t)-µ]²dt ⌡ −∞ 1 ⌠ T σ² = --------------- = lim ---⎮ [x(t)-µ]²dt Eq 4 ⌠∞ T→∞ 2T ⌡ -T ⎮ dt ⌡ −∞ It is easy to expand the above equations to show that Ψ² = σ² + µ² The discrete representation of the variance is just m σ² ≈ 1/m ∑ [xi-µ]² i=1 Eq 5 The positive square root of the variance is called the standard deviation, σ. Gaussian or Normal Distribution: One often hears about a normal or gaussian distribution, where the process is random, but the observations are grouped about the mean with a greater probability that they are nearer the mean, than farther away. This distribution is often found in geophysical processes and is often the statistics assumed by a processes when we calculate the statistics. When one samples the real environment and calculates the statistics of a process, the resultant is a good approximation of a Gaussian distribution. The probability density distribution function for a process x(t) with a mean µx and standard deviation σx is given by p(x) = (σx √2π)−1 e-(x-µx)²/2σx² The normal or Gaussian distribution is plotted below in Figure 3 for σx = 1 and µx = 0. 6 N o rm a l o r G a us s ia n D is trib utio n 0 .4 0 .3 5 Probability Density 0 .3 0 .2 5 0 .2 0 .1 5 0 .1 0 .0 5 0 -5 0 x 5 Figure 3. A Gaussian distribution of quantity x, giving the probability of the occurrence of a deviation in x from its mean value of 0. *** Assignment #1 *** A. Work with MATLAB until you are comfortable. Create a 2048-point sine wave with an amplitude of 2.0 and a period of 64. First you have to create an array of angles as the argument of the sine wave. Create the numbers from 1 to 2048 by i=1:2048 and multiply them by 2π and divide by 64 to get the argument. Take the sine of this and multiply by the amplitude, 2.0, to get your final series. You can change the initial phase of the sine wave by adding a constant to the angle series. Plot the resultant series and label the axes. B. Calculate the statistics of the series (maximum, minimum, mean, mean square, root mean square (rms), variance and standard deviation). Some of these computations can be done with standard MATLAB functions, and others you will have to create. Do these results make sense with what we have discussed in class? Write your own MATLAB “stats” function for calculating these statistics and outputting the results. C. The series you generated is deterministic, i.e. not a random series. Generate a 2048-point random series using MATLAB’s random number generator. If you use the RAND function, the series that has amplitude of 0 to 1.0 with uniform probability of any value between 0 and 1. If you use the RANDN function, the series is normally distributed about 0 – any value from –½ to ½. If you use the first, you will want to offset the series by 0.5 to bring the mean to zero. Then you can multiply the result by 0.1 to make a random series of amplitude 0.1, and add this to your sine wave. Calculate the statistics of this random "noise" series, and the sine wave with the random series added (this series is now random since you can not exactly predict the next value in time). Does this agree with what we have discussed in class? D. Plot the normal and uniform distribution series and discuss of this simple exercise in terms of the statistics and your understanding of time series so far. Note: Discussion is important! 7 Covariance. The auto-covariance function describes the dependence of values of the sample function at one time, on those at another time. The auto-covariance for the continuous series is a function of time lag, τ. ⌠∞ ⎮ [x(t)-µ][x(t-τ)-µ]dt ⌡ −∞ R(τ) = -----------------------⌠∞ ⎮ dt ⌡ −∞ Eq 6 1 R(τ) = lim --T→∞ 2T ⌠ T ⎮ [x(t)-µ][x(t-τ)-µ]dt ⌡ -T Note that R(τ) is an even function (symmetric about τ = 0, or the same at τ = -τ) with a maximum at τ = 0. The auto-covariance is useful in detecting a deterministic signal in the presence of random background noise. It is obvious that if τ = 0, then the autocovariance function, R(0) = σ², the variance. Correlation. A normalized form of the covariance function is often used. The autocorrelation function or auto-correlation is R(τ) ρ(τ) = -----R(0) Eq 7 where ρ(τ) is a dimensionless number with |ρ(τ)| ≤ 1. If the sample function has a predominate periodicity, then at lag τ corresponding to that period, R(τ) will have a relative maximum or minimum since the series will be shifted over one period and line up. Example. Consider the arbitrary function of time = A Cos(kx - ωt + φ) where A is an arbitrary amplitude of the sinusoid, k is the wavenumber (2π/wavelength), x is the horizontal coordinate in the direction of the sinusoid, ω is the frequency (2π/period), t is time and φ the initial phase. Then, since R(τ) is symmetric about t=0, we can take the integral over only positive times by doubling the right hand side, and holding x constant, 1 R(τ) = lim T→∞ T ⌠ T ⎮ A Cos(kx-ωt+φ) A Cos(kx-ω(t-τ)+φ)dt ⌡ 0 8 i. take t = 0, and changing the range of the integration from t=0 into blocks with t going from 0 to π/ω, and then multiplying by n where we let n go to infinity in the limit. Then T in the denominator becomes nπ/ω and nω ⌠ π/ω R(0) = lim --- ⎮ A² Cos² (kx - ωt + φ)dt n→∞ nπ ⌡ 0 and since Cos²(arg) where arg goes from 0 to π is ½, R(0) = ½A² and ρ(0) = 1.0 ii. take t = 2π/ω and again change the range of the integration as above 2π nω ⌠ π/ω R(--) = lim --- ⎮ ACos(kx-ωt+φ) ACos(kx-ωt+φ2π)dt ω n→∞ nπ ⌡ 0 ⌠ π/ω = ω/π ⎮ A² Cos²(kx-ωt+φ)dt ⌡ 0 2π R(--) = ½ A² and ω 2π ρ(--) = 1.0 ω iii. take τ = π/ω π ⌠ π/ω R(-) = ω/π ⎮ A Cos(kx-ωt+φ) A Cos(kx-ωt+φ-π)dt ω ⌡ 0 ⌠ π/ω = - ω/π ⎮ A² Cos²(kx-ωt+φ)dt ⌡ 0 π R(-) = -½ A² and ω π ρ(-) = -1.0 ω iv. Finally take τ = π/2ω π ⌠ π/ω R(--) = ω/π ⎮ A Cos(kx-ωt+φ) ACos(kx-ωt+φ-π/2)dt 9 2ω ⌡ 0 ⌠ π/ω = ω/π ⎮ A² Cos(kx-ωt+φ) Sin(kx-ωt+φ)dt ⌡ 0 π R(--) = 0 2ω and π ρ(--) = 0 2ω Again, it is clear that R(0) = σ², the variance. Power Spectrum. If x(t) is a time series made up of a sum of “m” number of cosines each with its own amplitude, Ai, frequency, ωi and phase θi, i.e. m x(y,t) = ∑ Ai Cos(kiy - ωit + φi) i=1 Eq 8 where we have also included a wavenumber, ki and distance, y. Then the variance is twice the integral from 0 to T in the limit as T →∞, 2 ⌠ T m σ² = R(0) = lim - ⎮ ∑ Ai²Cos²(kiy - ωit + φi)dt T→∞ 2T ⌡ 0 1 Now we can break the length T up in to n pieces, each π in length, and let the number of pieces go to infinity. n ⌠ π m σ² = R(0) = lim - ⎮ ∑ Ai² Cos²(kiy - ωit + φi)dt n→ ∞ nπ ⌡ 0 1 1 ⌠ π m σ² = R(0) = - ⎮ ∑ Ai² Cos²( kiy - ωit + φi)dt π ⌡ 0 1 m = 1/π ∑ Ai² π/2 i=1 m σ² = R(0) = ∑ ½Ai² i=1 10 Eq 9 So if x(t) can be regarded as being made up of a sum of sinusoids, its variance can be decomposed into components of average power, ½ Ai², at the various frequencies, ωi. Assuming a continuous distribution of frequencies, we obtain (without proof), σ² = R(0) = ⌠ ∞ ⎮ S(f)df ⌡ −∞ Eq 10 where S(f) is called the power spectrum (or variance spectrum). Thus S(f) df is the measure of the average power or variance in the frequency band f - ½df to f + ½df – which really says that S is how the variance is distributed with frequency. It can further be shown that S(f) = ⌠ ∞ ⎮ R(τ) e-2πifτ dτ ⌡ −∞ Eq 11 This can be recognized as the Fourier transform of the covariance function (see below for definition and discussion of the Fourier transform). Then we must have R(τ) = ⌠ ∞ ⎮ S(f) e2πifτ df ⌡ −∞ Eq 12 with τ = 0, we again obtain equation 10. Cross-Covariance and Cross-Correlation. Given two different sample functions, x and y, with means of µx and µy, the cross correlation and cross covariance function can be taken as discussed above. We had from equation 6, ⌠∞ ⎮ [x(t)-µ][x(t-τ)-µ]dt ⌡ −∞ Rx(τ) = -----------------------⌠∞ ⎮ dt ⌡ −∞ for one series. For two series this just becomes ⌠∞ ⎮ [x(t)-µx][y(t-τ)-µy]dt ⌡ −∞ Rxy(τ) = ----------------------⌠∞ 11 Eq 13 ⎮ dt ⌡ −∞ Again, this can be normalized to give the cross-correlation, Rxy(τ) ρxy(t) = ------------√[Rx(0) Ry(0)] Eq 14 and again |ρxy| ≤ 1.0 The cross-correlation of two sets of data describes the dependence of the values of one set of data on those of the other set as a function of lag, τ. Note that now Rxy(τ) is not an even function and the maximum does not necessarily occur at τ = 0. i.e. consider x = cosine and y = sine. When τ = π/(2ω), they line up so you get a peak in ρ, so the maximum occurs at τ=π/(2ω) and not τ=0. The discrete covariance functions are, m Rxτ = 1/m ∑ [xi - µx][xi-τ - µx] i=1 Eq 15 m Rxyτ = 1/m ∑ [xi - µx][yi-τ - µy] i=1 Eq 16 The correlation is again the normalized covariance ρx = Rx/σ² Eq 17 ρxy = Rxy/√[Rx(0) Ry(0)] Eq 18 *** Assignment #2 *** A. Create two sine waves of the same frequency and but slightly different initial phase. Add in two different random noises (uniform and normal distribution), and calculate the Autocovariance and Auto-correlation functions for one of these series. Hint: the MATLAB CONV (convolution) function does the multiplication and summing as required for the covariance and correlation. The MATLAB COV function returns the variance of the vector, not the covariance. If you use MATLAB’s covariance function you will get a different result. Therefore, before using any of the MATLAB function, look at it and understand what it does. This is the power of MATLAB - the documentation of a function is always available as an “m” file. If you make your own, describe the results you get in terms of finite record length. 12 B. Calculate the Cross-covariance and Cross-correlation functions from these two series. Again describe the results. Do you understand what the results are telling you and how you might use this tool in your research? C. Plot the original time series, the auto- and cross-correlation results. Are the results what you expected from our discussion in class? Fourier Transforms Time-frequency space - In the case of a sine wave, everyone recognizes that it can be expressed by an amplitude and phase at a specific frequency, and that this representation is more representative of the geophysical process than expressing the sine wave as a function of time. i.e. instead of x(t) now we have a separate representation of the process which we express as X(f) where f is the frequency. X(f) has an amplitude, A, and phase, θ, associated with each frequency, f. However, it should be obvious that expressing a sine wave as a function of time or as a function of frequency are really just two ways of expressing the same thing. Furthermore, it makes more physical sense to look at the sine wave as an amplitude and phase at a frequency than as a function of time. We can extend this to say that our continuous process x(t) can be represented by a sum of sinusoids with a given amplitude and phase at each frequency. This concept is now quite accepted, but in 1807 when Fourier first suggested that one could reconstruct exactly any function of time by an infinite sum of sinusoids, he astounded many contemporaries. Fourier's theorem states that any function can be reconstructed from a sum of sinusoids and that the Fourier transform is exactly that sum. Therefore, this says that we can express an observed process by a function of time, or by a Fourier transform of this as a function of frequency. For our discrete series, we have a time series in “time space” or a Fourier series in “frequency space.” Fourier Transform For our function x(t) we define its Fourier transform Z(f) as Z(f) = ⌠ ∞ ⎮ x(t) e-2πift dt ⌡ −∞ Eq 19 Where Z is generally a complex number. It is made up of real and imaginary parts which are in reality the coefficients of the cosine and sine functions as Z(F) = A(f) + i B(f) x(t) = A(f) Cos(2πft) + B(f) Sin(2πft) Therefore, a series x which is made up of a sum of cosine waves only will have only a real part or real coefficients, (only A coefficients) and be an even function symmetric around 0. Similarly a series x which is made of a sum of sine waves only, will have a Fourier transform which is only imaginary (only B coefficient) and be an odd function. Depending on the normalization of the discrete Fourier transform, the coefficients A and B are the amplitudes of the cosine and sine waves which go to make up the time series x. The inverse Fourier transform of Z(f) is then defined by 13 x(t) = ⌠ ∞ ⎮ Z(f) e2πift df ⌡ −∞ Eq 20 These two expressions (equations 19 and 20) define a Fourier transform pair. It is obvious that (1) Z(f) is a function of frequency, (2) Z(f) is an exact mathematical representation of x(t) (which was our original series as a function of time), and (3) Z(f) is defined by the Fourier transform of x(t). Similarly x(t) is a representation of Z(f) as given by the inverse Fourier transform. The Fourier transform is the way of getting from a function of time to a function of frequency. Mathematically the Fourier transform of a function x(t) exists if ⌠∞ ⎮ |x(t)|dt ⌡−∞ exists. Some useful functions do not have transforms. x(t) = a non-zero constant x(t) = A Sin(2πft) x(t) = 1 for t > 0 = 0 for t < 0 However, the problem exists because these series are infinite in length. A time series resulting from a geophysical process which is sampled for a finite period of time (or what we call “gated”) always has a Fourier transform. Mathematically, we are saying that any gated function (a function which has a beginning and end and no infinite values) has a Fourier transform. Since our observations start and stop in a finite time interval, they have Fourier transforms. Distance-wavenumber representation. Usually one represents processes as a function of time and uses the Fourier transform to obtain a representation as a function of frequency. Similarly one can define a Fourier transform pair which are functions of distance and wavenumber. So variations as a function of space (as for example vertical CTD profiles of temperature and salinity in the Moonakis river) can be represented as a space series, or a wavenumber series. The Fourier transform is the mechanism relating these two representations of the process. For example, if we have an observed series such as a surface wave field of the form Y(x,t) = A Cos(kx - ωt) then we have two Fourier transforms, first as a wavenumber spectrum at one point in time (constant t) Z(k) = ⌠ ∞ ⎮ Y(x) e-2πikx dx ⌡ −∞ 14 or as a frequency spectrum at one point in space( constant x) Z(f) = ⌠ ∞ ⎮ Y(t) e-2πift dt ⌡ −∞ Fourier transform relationships - Let x(t) and Z(f) be a Fourier transform pair, and be complex functions of real variables. We will use the shorthand x(t) ⊃ Z(f) and also Eq 21A ____ Z(f) = x(t) Eq 21B then, x(-t) ⊃ Z(-f) Eq 22 x*(-t) ⊃ Z*(-f) Eq 23 where * denotes the complex conjugate. Note that if we have a complex number defined as follows, Z = R + i I, where then i = √−1 Z* = R - i I So from Eq 19 and Eq 20, we can then say if x(t) is even, x(t) = x(-t) ⊃ Z(f) = Z(-f) which is even x(t) is odd, x(t) = -x(-t) ⊃ Z(f) = -Z(-f) which is odd X(t) is real, x(t) = x*(-t) ⊃ Z(f) = Z*(-f) which is Hermitian X(t) is imag. x(t) = -x*(-t) ⊃ Z(f) = -Z*(-f) which is anti-hermitian Define X(f) and Y(f) as the Fourier transforms of x(t) and y(t) respectively, and let “a” be a scalar. Some useful theorems concerning Fourier transform pairs are, x(t) + y(t) ⊃ X(f) + Y(f) additive a x(t) ⊃ a X(f) x(t-a) ⊃ e-2πiaf X(f) 15 shift x(at) ⊃ 1/|a| X(f/a) scale ∂x(t)/∂t ⊃ 2πif X(f) ⌠ ∞ ⎮ x(t)y*(t)dt = ⌡ −∞ ⌠ ∞ ⎮ |x(t)|²dt = ⌡ −∞ derivative ⌠ ∞ ⎮ X(f)Y*(f)df ⌡ −∞ power ⌠ ∞ ⎮ |X(f)|²df ⌡ −∞ We shall return to these relationships between transform pairs when we consider filtering. *** Assignment #3 *** A. Create a sinusoidal time series of amplitude 2.5 to which you have added a random component of amplitude 0.2. Note that the length of you series should be a power of two. B. Calculate and plot the Fourier transform (FFT) of your original series and your noisy series. Does it agree with what you expect? Note that MATLAB normalizes the transform by the length of series. I like to normalize it so the FFT returns the amplitude of the sine or cosine. C. Create your own MATLAB function (which we will eventually evolve into the power spectrum and coherence function) with proper normalization. Convolution Product and Filters Convolution Product: To discuss filters we need to define the convolution product which we will use to represent the filtering process. The convolution product of our series x(t) with some arbitrary set of weights w(t) is defined by the integral Y(t) = ⌠ ∞ ⎮ w(τ) x(t-τ)dτ. ⌡ −∞ Eq 24 We will abbreviate the convolution product by the symbol, *, and so Y(t) = w(t) * x(t) = x(t) * w(t) Eq 25 Then we can show that this product is g * h = h * g communitive f * (g * h) = (f * g) * h associative f * (g + h) = f * g + f * h 16 distributive Eq 26 Eq 27 Eq 28 Convolution theorem If we denote the Fourier transform of x(t) by X(f) or x ⊃ X, as defined in Eq 21, then the convolution theorem states that (remember f ⊃ F and g ⊃ G) f * g ⊃ F • G f * g = f • g or Eq 29 f * g ⊃ f • g. This says that the Fourier transform of the convolution product of f and g is equal to the product of the Fourier transform of f and the Fourier transform of g. Conversely, f • g = f * g Eq 30 and we also have f * g * h = f • g • h f * (f • h) = f • (g * h) Eq 31 Eq 32 A time series can be expressed either in the time domain or by its Fourier transform in the frequency domain. The convolution theorem tells us that we can substitute multiplication in one domain for convolution in the other. It is easier to visualize the product of two functions than the convolution, and it is often easier (quicker) to compute. The convolution theorem tells us that if we express the time series x as a Fourier transform, and multiply this Fourier transform by the Fourier transform of w, then inverse transform, we will obtain the same result as doing the convolution of the filter w and the series x. This is written w * x = w • x Eq 33 With the FFT (fast Fourier transform) algorithms available on digital computers, it is often faster (and hence cheaper) to transform, multiply, and inverse transform, rather than calculate the convolution. (For information on the introduction of the FFT, see the IEEE Trans. Audio and Electroacoustics, June 1967, Vol AU15, vol2 Pg 43-93.) The convolution process is the filtering process where a filter, w, is applied to the series, x. If the time series x(t) is a discrete series (see section following on the sampling process for details on how the discrete series is created), xt, t=0,1,2,…,n, and the filter, wt, t=0,1,2,3,…,m, is defined over a different number of terms, m, but with the same sample interval, then the convolution product of the discrete series can be represented as the finite sum Yt m = ∑ wτ xt-τ for t = 0,1,2,…n+m τ=0 17 This will only have meaning if we define xt = 0 ; t=n+1,n+2,..,n+m wt = 0 ; t=m+1,m+2,..,n+m For example, consider the filter to be a triangular filter represented by a set of 5 weights, wt = {0, 0.25, 0.5, 0.25, 0}, and let the series be represented by a set of numbers, xt = {1.2, 2.2, 1.6, 1.8, 1.2, 1.6, 1.4}. Note that the filter series weights sum to one. This is so that the values of the series will have the same value after the filter is applied. Then n=6 and m=4, and Y0 = 0.0 + 0.55 + 0.8 + 0.45 + 0.0 = 1.8 Y1 = 0.0 + 0.4 + 0.9 + 0.3 + 0.0 = 1.6 Y2 = 0.0 +.0.45 + 0.6 + 0.4 + 0.0 = 1.45 Note that the filtered series is shorter than the original series by m-1=3 terms so the filtered series is now n-(m-1) = 3 terms. It is obvious that x * w = w * x since the sum is symmetric. Filter Example 3 2.5 X or Wor Y 2 1.5 1 0.5 0 0 1 2 3 4 Time 5 6 7 8 Figure 3. The example given in the text plotted. Filters and filtering: Filters are selective devices which are used to discriminate and reduce time series. In a broad sense, the world, our experiment, the sensors, recorder, and data analysis are really filters which shift and reduce the data. We begin by considering some simple filters such as would be used in data reduction. Consider a filter as a box with an input x(t) and an output y(t) x(t) → F → y(t) 18 We shall write the filtering operation as x(t) => y(t). Time Invariant Linear Filters. We want well behaved filters which will enhance some frequencies in a time series and suppress others. This type of filter is invariant in time, that is, if x(t) => y(t) Eq 34 then x(t+τ) => y(t+τ) for all τ. and is linear, that is if and x1(t) => y1(t) x2(t) => y2(t) Eq 35 then x1(t)+x2(t) => y1(t)+y2(t) For any constant, a, it follows that a x(t) => a y(t) Eq 36 A time invariant, linear filter has the important property that it does not confuse frequencies (no non-linearities). A sinusoidal input produces a sinusoidal output of the same frequency, that is A Cos(kx-ωt+φ) => B Cos(kx-ωt+θ) where A, B, ω and θ are constant in time. The filter has altered the amplitude from A to B (a gain of B/A) and altered the phase by θ - φ, but the frequency of the sinusoid is unchanged. In order to describe the response of a filter, let us observe what it does to an ideal input. We want to consider a complex filter so, consider two parallel filters A1 Cos(ωt) => y1(t) A2 Sin(ωt) => y2(t) and regard the first as the real and the second as the imaginary part of this complex filter, then we can rewrite the two inputs as one and linearity gives A1 Cos(ωt) + i A2 Sin(ωt) => y1(t) + i y2(t) This can be written, A eiω(t+τ) => y(t+τ) Eq 37 where the phase is written as a frequency, ω, times a constant, τ. This can be broken up into the sinusoidal part and a phase part as 19 A eiω(t+τ) = A eiωt eiωτ => y(τ) eiωt where exp(iωt) is regarded as a multiplier. For τ = 0, or zero phase shift, the above relationship gives y(t) = y(0) eiωt Eq 38 Where we now have the output as an initial amplitude which is expressed as y(0) and the sinusoidal time varying part expressed as the exponential, so that A ei(ωt) => y(o) ei(ωt+θ) Eq 39 The real part of Eq 39 is A Cos(ωt+θ) => |y(0)| Cos(ωt+θ+φ) Eq 40 B Cos(ωt+ψ) where B = |y(0)|, ψ = θ + φ, y(0) = |y(0)|exp(iφ). We characterize the response of a filter (as above) by its transfer function, L, by use of Eq 38 L(ω) = y(0)/A Eq 41 which is the ratio of the output to the input at a given frequency ω . L(ω) is generally complex, it has a real and imaginary part, which can also be written as an amplitude gain and a phase shift. In terms of the real input Eq 40 L(ω) = B/A eiφ = B(ω)/A(ω) eiφ(ω) so that L(ω) is a complex number whose modulus is the amplitude gain (B/A) and whose argument is the phase shift (φ) produced by the real input at each frequency ω. A property of L(ω) is that L(ω) = L(-ω)* Eq 42 This transfer function, L(ω), characterizes how a filter responds to a sharp spike in frequency (a sinusoid). This is done in the laboratory, but applying a sine wave and observing the output voltage amplitude and phase change. An alternate description can be made by specifying its response to a spike or impulse in time. Let δ(t) represent a unit pulse in time. This pulse is narrow compared with the resolving power of our sensors. For example take the gate function Π = 1 Π = 0 for for |t| ≤ ½ |t| > ½ 20 Eq 43 This represents an impulse when written as 1/τ Π(t/τ) = δτ(t) as τ goes to 0, δτ(t) becomes arbitrary sharp (centered near t=0). We have ⌠ ∞ ⎮ δτdt = lim ⌡ −∞ τ→0 ⌠ ∞ ⎮ (1/τ)Π (t/τ)dt = 1 ⌡ −∞ Eq 44 and as Π is non zero between -½τ and ½τ, ⌠ ½τ = lim ⎮ (1/τ)Π(t/τ)dt = 1 τ→0 ⌡ −½τ Eq 44 It can be shown by the limiting process that ⌠ ∞ ⎮ δ(t) F(t) dt = F(0) ⌡ −∞ and by the shifting property of δ(t) ⌠ ∞ ⎮ δ(t-a) F(t) dt = F(a). ⌡ −∞ If δ(t) => l(t), we say that the impulse δ(t) produces an impulse response, l(t). An input which can be expressed as a linear combination of impulses, has an output which is described by the impulse response function. m m ∑ Ciδ(t-ti) => ∑ Cil(t-ti) i=1 i=1 Assuming that our impulses are infinitely close, we can go to the integral representation, and make the input exactly any input function of time. Hence, we can express a filter by either its response to a sum of sinusoids (in the frequency domain, L(ω)) or by its response to a sum of impulse (in the time domain, l(t)). This can be expressed, if eiωt => L(ω)eiωt 21 then Eq 45 ⌠ ∞ x(t) = ⎮ g(ω)eiωtdω ⌡ −∞ > ⌠∞ ⎮ g(ω)L(ω)eiωtdω = y(t) ⌡ −∞ where g(ω) is the Fourier transform of x(t), so that x(t) is the inverse Fourier transform. This results is an integral representation of the transform, which is the amplitude and phase of the sinusoids at each frequency which add up to the initial function of time, which are now shifted in amplitude and phase by the complex response function, L(ω). Now we also can write, δ(t) => l(t) Eq 46 then x(t) = ⌠ ∞ ⎮ x(t)δ(t-τ)dτ => ⌡ −∞ ⌠∞ ⎮ x(τ)l(t-τ)dτ = y(t) ⌡ −∞ Therefore there are two ways of describing the response of a filter, as its frequency response function, L(ω) or as its impulse response function, l(t). It should be obvious that these are just Fourier transforms of each other, L(ω) ⊃ l(t) and l(t) ⊃ L(ω). Examples of filters 1. The do-nothing filter, x(t) => x(t) and L(ω) = 1 (e.g. amplitude multiplied by 1 & phase shifted by 0º) l(t) = δ(t) (e.g. you get out what you put in) 2. The lag filter, x(t) => x(t-τ). The output is simply a delayed version of the input. L(ω) = eiωτ l(t) = δ(t-τ) 3. Consider x(t) => y(t) = x²(t). An input consisting of two different frequencies produces sums, and difference frequencies of the inputs. These are called distortion (harmonic, intermodulation) and are the result of squaring a linear system. This is a bad filter to have to describe based on the output. 4. A time variant filter. If x(t) => g(t)x(t) where g(t) is any non-constant function of time. Such filters are commonly used on time series. Consider the gate function or “Boxcar,” 22 g(t) = Π(t). Since observations must be started and stopped, multiplication by Π(t) is unavoidable. We generally scale so the filter is Π(t/T) where the record length = T, so the resulting frequency confusion does not negate the experiment. Let us more fully explore the gate function, II(t), as a filter. It has a Fourier transform Π(t) ⊃ ⌠ ∞ ⎮ Π(t) Cos(2πft)dt ⌡ −∞ where we need only to use the Cosine part of the transform since the gate function is an even function. Since the gate function is zero outside the interval |1/2|, the limits of the integral become plus and minus 1/2. We can further simplify this by noting that an even function is twice its value from 0 to 1/2. Therefore, Π(t) Π(t) ⊃ ⊃ 2 |½ --- Sin(2πft) | 2πf |0 Sin(πf) ------- ≡ Sinc(f) πf Applying this filter by multiplication in the time domain (our gating process) is the same as convolving the Fourier transform of your observations with the sinc function. X(t) • Π(t) ⊃ X(f)* Sinc(f) This often produces problems and must be dealt with properly in the design and sampling of the experiment. Therefore, filters are actually applied to the data in the process of sampling, and their effects must be known. To study these effects we must examine the sampling process. *** Assignment #4 *** A. Create a sine wave with noise as we did in assignment 3 and calculate the Fourier coefficients for your sine wave. Is the energy at the frequency it is supposed to be? Is your amplitude and phase properly represented in the Fourier coefficients? C. Create a 5-point triangular filter of weights 0.0, 0.25, 0.50, 0.25, 0.0. To be normalized as a low pass filter, the weights must sum to one. Apply your triangular filter to your sine wave with noise. D. Plot the initial and filtered series. Does the filter remove some of the random noise? E. Using FFT, transform the filter weights (cosine transform or real part since the filter weights are an even function) to get the filter gain (filter frequency response function). Make a log-log plot of your filter's frequency response function. Does it look reasonable to you? Discuss. 23