Fourier Transforms

Properties of Fourier Transform

Match Filter

Sampling

Derivatives Using the Fourier Method

Rayleigh’s theorm

Here the original function, f(t) varies with time and its Fourier transform, F(w), varies with frequency, w = 2?f. The general form of the Fourier transform for a function of many variables (g1, g2, g3, ....) is:

we can also write:

for a function that varies with position, x, and its Fourier transform varies with wavenumber, kx. For a function of time and space, f(t,x), we have

where kx is the spatial frequency and w is the temporal frequency. This representation is very useful in wave propagation problems.

Slow Fourier Transform:

Below is a simple function to calculate a ‘slow’ Fourier transform at one frequency OMEG.

COMPLEX FUNCTION XFOMEGA (X, N, OMEG)

C Calculates the complex Fourier coefficient for one frequency, OMEG.

C The input is a real data vector, X, of length N.

REAL X (*), OMEG

INTEGER N

XFOMEGA = CMPLX (0., 0.)

DO J = 1, N

XFOMEGA = XFOMEGA + X (J) * CEXP (CMPLX (0., (J-1)*OMEG)

END DO

RETURN

END

1) Scalar multiplication:

a f(t) <----------->a F(w)

let

2) Addition of two signals:

h(t) = f(t) + g(t) <-----------> F(w) + G(w) = H(w)

then

e.g., let f(t) = d(t) f(t) <-----------> F(w) = 1

g(t) = .5 d(t – t1) g(t) <-----------> G(w) = .5 e–iwt1

clearly, h(t) = f(t) – g(t) <-----------> F(w) – G(w) = H(w)

3) Shift theorem:

f(t–a) <-----------> e–iaw F(w) .

i.e., if we delay a signal in time, we add a linear phase shift "a" to its Fourier transform.

Consider

let t' = t – a

dt' = dt

and

t = t' + a

then

4) Differentiation:
 

Use òby parts, remember ò udv = uv – òvdu

let u = f(t) and du = e–iwt dt

in òf(t) e–iwt dt

if dv = e–iwt dt, then

So, uv = 0, and we have

So, if F(w) = u(w) + iv(w)

iw F(w) = iwu(w) – wv(w)

= w(–v(w) + iu(w))

For computation, given F(w)

Re -------> IM

–IM------> Re

and amplify by w

So, differentiation in time corresponds to a linear amplifier of frequency, so high frequencies are amplified.

By successive applications of above proof we can show:

e.g.,

etc.

5) Integration:

again use integration by parts

let u = e–iwt

then v(t) = ò f(t) dt

So, ò udv = uv – ò vdu

where

So, integration is the opposite differentiation, i.e damping of high frequencies using a linear attenuation.
 

so, Re -------> IM

–IM------> Re

and divide by w

let

let t' = at

dt' = a dt

then

Note

so,

7) Modulation:

remember

so,

8) Convolution:

let

Interchange integrals

the shift theorem says

so that

which we recognize as the Fourier transform of g(t).

So that

= G(w) F(w)

So, convolution in time becomes multiplication in the frequency domain.

c convolution via fft

parameter (nmax=1024)

real x (nmax+2),y(nmax+2)

complex xc(nmax /2+1),yc(nmax/2+1)

equivalence (x(1),xc(1)),(y(1),yc(1))

c read the data values into x and y

c nx=number of x values

c ny=number of y values

read(*,*)nx, ny

:

:

:

c determine 2**m >nx and ny

do j=2, 16

nftx=j**2

if (ntx.gt.nx) go to 10

end do

10 do j=2, 16

nftx=j**2

if (nfty.gt.ny) go to 15

end do

15 nft=max (nftx,nfty)

c add zero pad

do j=nx+1,nft+2

x(j)=0

end do

do j=ny+1,nft+2

y(j)=0

end do

nfreq=nft/2+1

c get the fft

call four2 (x,nft,1,–1,0)

call four2 (y,nft,1,–1,0)

c convolve via complex multiply

do j=1,nfreq

xc(j) ´ yc(j)

endo

c go back to time

call four2(yc,nft,1,+1,–1)

Let

and

Then
 

or, let

then

so

amplitudes are multiplied and phases add

Correlation is convolution with one signal time reversed. Assume that f(t) and g(t) are real signals:

interchange ò's again

from shift theorem

If we let
 

and

or in polar form

So, amplitudes are multiplied and phases are subtracted

Autocorrelation:

let g(t) = f(t)

then

so that

H(w) = F(w) F(w)*

or in polar form:

but

so the transform is

= A_f²                                        i.e., power spectra

no phase information is in the ACF.

Assume that h(t) = w(t) * r(t),
where w(t) is a wavelet, r(t) is the reflectivity and assume that w(t) is known.
Then, if we correlate with w(t), we can find the 'matches,' i.e. they will occur where the new function is large.

h(t) = w(t) * r(t)

or

H(w) = W(w) R(w)

let

g(t) = [w(t) * r(t)] * w(–t)

G(w) = W(w) R(w) W*(w)

= [W(w) W*(w)] R(w)

In Vibroseis w(t) is a swept frequency, so data can not be recognized as seismic data until it is correlated with the sweep.

? is sweep generated actually what goes into the earth ?

? what kind of problems will this cause ?

Sampling:

The delta function, d(t), is an impulse at time, t, and zero every place else. It can be used to select data from a function at any specified time when the argument of the delta function becomes zero. For example:

f0 = d(t) f(t)

f1 = d(t–1Dt) f(t)

fN = d(t–NDt) f(t)

where Dt is some specified time interval, e.g., 1 s, 1 ms, etc.

When data are collected, it begins as analog (continuous) measurements of time or space during data acquisition. It is implicitly multiplied by equi-spaced delta functions to obtain discrete data. These delta functions are usually evenly spaced at a pre-determined sample interval, e.g., Dt, known as the sample rate.

The Nyquist frequency is defined as 1/ (2 ´ sample interval time).

fNyq =1./ (2Dt)

It is important because frequencies above the Nyquist frequency will be aliased. That is, frequencies greater than Nyquist will be lost from the original signal and appear as lower frequencies, i.e., they will be aliased.

The discrete data set can be represented as a polynomial z transform, where z=exp(iwDt) and w=2pf, where f varies from 0 to Nyquist = 1./ (2Dt).

c bandpass filter program

parameter (nmax=1024)

real x(nmax+2)

complex xc(nmax/2+1)

equivalence (x(1),xc(1))

c read the data

c the number of real data values is ndata

c the sample rate is dt

read(x,x)ndata,dt

:

:

:

c read the low cut and high cut filter

c frequencies

read(*,*)flow,fhigh

c 2**m>ndata

do j=2,16

nft=2**j

if (nft.gt.ndata) go to 10

end do

c zero pad

10 do j=ndata+1,nft+2

x(j)=0.0

end do

nfreq=nft/2+1

fnyq=1./(2*dt)

c low, high frequency cuts

if (flow.gt.0) then

ilow=nfreq* (flow/fnyq)

ilow=max (ilow,1)

ilow=min (ilow,nfreq)

end if

if (fhigh.lt.fnyq) then

ihigh=nfreq* (fhigh /fnyq)

ihigh=max (ihigh ,1)

ihigh=min (ihigh,nfreq)

end if

c fft

call four2 (x,nft,1,–1,0)

c apply filter to xc which is equivalent to x

c bandpass, low end first

if (flow.gt.0) then

do jf=1,ilow

xc(jf)=cmplx (0.,0.)

end do

end if

if (fhigh.lt.fnyq) then

do jf=ihigh, nfreq

xc (jf)=cmplx (0.,0.)

end do

end if

c apply inverse fft

call four2 (xc,nft,1,+1,

Note that the Fourier transform of f(t) is F (w) and vice versa:

If:

then:

and

The derivative of a function, f(t)

has an effect on the original signal, F(w) in the frequency domain. It acts like a high pass filter. The low frequencies in the frequency domain are damped, while the higher frequencies are emphasized.

F(w)

w

The integral of a function, f(t)

also has an effect on the original signal, F(w), in the frequency domain. It acts like a low-pass filter. The high frequencies in the frequency domain are damped, while the lower frequencies are emphasized (see below).

F(w)

w

Integtral of the squared modules of a function is equal to the integral of the squared modules of its spectrum.
 

now consider
 

if w¢ = 0

Now use the Fourier transform to get
 

and since w¢ = 0