NUMERICAL ANALYSIS of PROCESSES

NAP3

Black box models and experimental identification of linear systems characteristics:

Transfer characteristics.

Volterra equation and transfer function identification. Convolution and deconvolution.

Fourier analysis (continuous and discrete FT).

Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

NAP3

Black box models based upon transfer characteristics express relationships between inlet x(t) and outlet y(t) of a system, for example

relationship between humidity of air (t) and moisture of dried material X(t)relationship between time course of temperature and concentration of microorganism relationship between inlet and outlet concentration in continuous chemical reactorrelationship between a valve closing and pressure in a pipeline (water hammer) relationship between temperature of fluid and recorded signal of thermocouple

Transfer characteristics

x(t) y(t)

x

y

t [s]

NAP3

E(t)x(t) y(t)

x

y

t [s]

E

At nonlinear systems even a small disturbance of stimulus function can cause instabilities (deterministic chaos). This subject (nonlinearities, chaos, strange attractors…) will be subject of the following lecture. Nothing like this exist at linear systems, where any relatinship b etween stimulus function and response can be expressed by linear integral relationship – by convolution

dxtEty )()()(

Response y(t) is a convolution of stimulus x(t) and transfer function E(t), symbolic form

y=E*x

For explanation, see next slide…

Transfer characteristics

Dirac function can be expressed as the limit of

Gauss function

NAP3

E(t)(t)

E(t)

Let us consider special case of stimulus function x(t) in form of infinitely short impulse with unit surface. This is so called delta function (t) equal zero for t0 and infinity for t=0. Integral from minus to plus infinity is 1, therefore the convolution integral is directly the E(t) function:

( ) ( ) ( ) ( )y t E t d E t

This is why the E(t) is

called impulse responsex=(t)

t [s]

E

Impulse response can be determined experimentally by measuring response to to a short impulse, for example to an instantaneous injection of trace to the inlet. Impulse response E(t) is the measured concentration of tracer at outlet..

2

lim)( nt

nen

t

Transfer characteristics

NAP3

E(t)x(t) y(t)

General stimulus function can be replaced by a sequence of infinitely many short pulses with area x()d. Each impulse is associated with the impulse response shifted by the time delay and multiplied by the area x()d (linear system –response is directly proportional to input). Sum of these elementary responses is the convolution integral E * x

0

( ) ( ) ( ) ( ) ( )t

y t E t x d E t x d

x()d

t [s]

E(t-)x()d

Impulse response and x(t) are zero for negative times

Transfer characteristics

NAP3

There exist two principle problems

Convolution, when characteristic of system E(t) and time course of inlet x(t) are know. Response y(t) is obtained by integration. At each time t it is necessary to calculate a new integral.

Deconvolution, when input x(t) and output y(t) are known and impuls function E(t) should be calculated (system identification). In this case it is necessary to solve integral equation.

A suitable tool for solution of these two problems is Fourier transform

Volterra equation

dxtEty )()()(

Leger

Fourier transformation

Duchamp

NAP3

The functions x(t),y(t),E(t) could have been expressed by ordinary polynomials (see regression) and the convolution as a relationship between coefficients of these polynomials. However, each polynomial goes to infinity for t and regression polynomials of order greater than 7 cannot be identified with sufficient accuracy (this is because the functions 1,t,t2,t3,…,t7,t8 are of a similar shape and therefore are not sufficiently linearly independent. Then the solution of coefficients is distorted by round off errors).

Instead of common polynomials it is better to use orthogonal functions, for example orthogonal polynomials or goniometric functions and to approximate functions x,y,E by Fourier expansions

0

( ) ( )j jj

E t E P t

But what does it mean the orthogonality of functions Pm(t) and Pn(t)?

You should know it, see previous lecture ☺.

For example goniometric functions (sin, cos) are orthogonal in the interval (-1,1).

Orthogonality of Pm(t)=cos 2mt follows from1 1

1 1

11

1cos 2 cos 2 (cos 2 ( ) cos 2 ( ) )

2

1 sin 2 ( ) sin 2 ( )[ ] 1

2 2 ( ) 2 ( )

mt ntdt m n t m n t dt

m n t m n t

m n m n

for m=n, otherwise =0

Orthogonality of sin 2mt can be proved in the same way. Using Euler’s formula follows orthogonality of the complex exponential function

2( ) cos 2 sin 2i mtmP t e mt i mt

Prove!!!

i-imaginární jednotka

NAP3 OrthogonalityFunctions Pm(t) and Pn(t) are orthogonal, if integral of their product (across some interval) is zero

functions cos are even orthonormal, because integral is normalised

(equal one)

integral of product is called scalar product

Any function c(t) can be expressed as a Fourier expansion which is a linear combination of orthogonal functions Pj(t)

0

( ) ( )j jj

c t c P t

NAP3 Fourier expansion

where the coefficients of expansion are defined uniquely by function c(t) (and this association is called discrete Fourier transform )

2

( ) ( )

( )

j

j

j

c t P t dtc

P t dt

jc~

Prove!!!0 1( ) , ,..., nc t c c c

Coefficients of expansion can be evaluated without solving a system of algebraic equations (compare with the regression analysis using common polynomials). Application of orthogonal functions instead of P0(t)=1, P1(t)=t, P2(t)=t2,… is usually advantageous also in regression analysis, where only few initial terms of Fourier series is used. Legendre, Tschebyschev, Hermite or Laguerrovy generalised polynomials are used (differences follow from different definition of scalar product, for example Laguerre

polynomials have integration interval <0,>, Legendre <-1,1>. Selection of suitable orthogonal polynomial is therefore determined by the definition interval of regression function: for<0,>, Laguerre polynomials, for <-,> Hermite and for finite interval Legendre or Tschebyschev polynomials).

In the following text only the Fourier expansion based upon goniometric functions will be discussed 2( ) imt

mP t e (orthogonality on different interval than <-1,1> will be proved in a while)

Sum of infinite number of terms in Fourier series can be substituted by integral. Instead discrete indices j a continuous frequency f will be used. For example

2 2

0

( ) ( ) ( )ijt iftj

j

c t c e c t c f e df

NAP3 Fourier expansion

and this representation is called continuous integral Fourier transformation.

dtetcfc ift2)()(~

2( ) ( ) iftc t c f e df

Back transform

Big advantage is practical identity of forward and back transforms (the same routines are used for forward

and inverse transforms)!

FT is a complex function even if

the c(t) function is real

Fourier series Fourier integral

Continuous forward Fourier transform (time is substituted by frequency)

Proof (simplified)

NAP3 Fourier transform

2 2 2 2

2 ( ) 2 ( )2 ( )

( ) ( ( ) ) ( )( )

( ) ( ) ( ) ( )lim lim2 ( )

sin 2 ( )( ( ) )lim

( )

if ift if ift

a ia t ia tif t

a aa

a a

c t c e d e df c e e df d

e ec e df d c d

i t

a tc d

t

sin sin

( ( )) ( ) ( )lim2

x x xc t dx c t dx c t

a x x

proof that the functions e2if are orthonormal

Substitution 2a(-t)=x

dx=2 ad

PSD Power Spectral Density

22)(~)(~)( fcfcfPc

Energy of low frequencies

Energy of high frequencies (noise)

NAP3 Fourier transform

Graph PSD shows dominant frequencies and is suitable for noise filtering setup

Beauty of FT is in the simplification of difficult operations

1. integrals of products (convolution or correlation) are transformed to product of transformed functions

2. derivative of function is transformed to multiplication by frequency (see next lecture)

3. easy inverse transformation (formard/backward transforms can be realised by the same procedure)

NAP3 Fourier transform

dytxtCxy )()()( ( ) ( ) ( )xyC f x f y f

dytxtRxy )()()( )(~)(~)(~ * fyfxfRxy

Convolution of two functions (change of x(t) passing through the filter y(t))

Correlation of two functions (measure of similarity of mutually shifted functions x,y)

mean time -inlet x(t)=0.5

mean time out y(t)=2.1

mean time of cross-correlation Rxy(t)=1.6

represents time shift between x(t) and y(t)

NAP3 Fourier transform

Calculation of response to a given input signal x(t)

FT of stimulus function and impulse response (forward FT)

FT of response is the product of FT (product of complex numbers)

Response calculated by inverse FT

NAP3 FT aplications

)(~

)( ),(~)( fEtEfxtx )(

~)(~)(~ fEfxfy

( ) ( )y f y t

Identification of system from measured x(t),y(t)

Forward FT of input and response

FT of impulse response is ratio of FT (ratio of complex numbers)

Impulse response calculated by inverse FT

( ) ( ), y( ) ( )x t x f t y f

( ) ( ) / ( )E f y f x f

( ) ( )E f E t

Sampling of signalIn reality we operate only with a point representation of functions, sampled usually with the constant time step .

A function is represented by N-data points and N-frequencies (even number)

t: t0=0, t1=, t2=2, 3,......, tN-1= (N-1),

f: f-N/2=

2

1,f-N/2+1=

1)

1

2

1(

N,.. f-1=

N

1, f0=0,

f1=N

1,.. fN/2-1=

1)

1

2

1(

N, fN/2=

2

1

Nyquist frequency

Nyquist frequency is the highest frequency which can be distinguished in data sampled with the time step

2

NAP3

Discrete FT (Fourier Transform)

DFT has the same properties (concerning convolution and correlation) as the continuous FT.

NAP3

2( ) ( ) iftc f c t e dt

Discrete DFT

1

0

/2)()(~N

k

Niknkn etcfc

where

n=-N/2,-N/2+1,…,0,1,…,N/2

t: t0=0,t1 =,t 2=2 ,3,......,tN-1= (N-1),

f: f-N/2=

2

1,f-N/2+1=

1)

1

2

1(

N,.. f-1=

N

1, f0=0, f1= N

1,.. fN/2-1=

1)

1

2

1(

N, fN/2= 2

1

N

nfn ktk

Remark: DFT coefficient corresponds to unit time step =1, more accurately

nn cfc ~)(~

Discrete FT (Fourier Transform)

N-Fourier coefficients represent the function c(t) defined from - to +, passing through specified points c0, c1,…cN-1 and is periodic with the period N.

NAP3

1 1 1 12 ( /2)/

/20 0 0 0

( ) (cos( ) sin ) cosN N N N

ik N N ikN k k k k

k k k k

c f c e c e c k i k c k

sin(k)=0 for arbitrary k

Fourier coefficients of positive and negative Nyquist frequencies are identical:

1 1 1 12 ( /2)/

/20 0 0 0

( ) (cos sin ) cosN N N N

ik N N ikN k k k k

k k k k

c f c e c e c k i k c k

t0 t1 t2 ….. tN-2 tN-1 tN

N

c0

c1

c2

cN-1

Periodicity of Fourier coefficients follows from NiknNnNki ee /2/)(2 12 ine

Negative frequencies?NAP3

What is relationship between the index of a FT coefficient and frequency?

Example: 8 points sampled with frequency 1 Hz, =1s. Nyquist (max.) frequency is 0.5 Hz. Corresponding FT vector has NINE coefficients (c-4,c-3,…c0,…c3,c4) for frequencies (-4/8,-3/8,-2/8,-1/8,0,1/8,2/8,3/8,4/8). However, FT coefficients of positive and negative Nyquist frequencies are identical (c-4=c4) so that the number of FT coefficients is in fact only eight (as it should be!). These coefficients are usually arranged in a vector starting with ZERO frequency, i.e. (c0, c-1, c-2, c-3,c4=c-4, c3, c2 ,,c1), see MATLAB.

0 1 2 3 4 5 6 7-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

4sin(2 ( ) )

8t

0 1 2 3 4 5 6 7-1

-0.8

-0.6

-0.4

-0.2

0

0.2

0.4

0.6

0.8

1

When transforming a real function it is sufficient to consider only positive (or negative) frequencies because FT coefficients of positive and negative frequencies are not independent (they are complex conjugated). FT is 1:1 transformation, therefore 8 real numbers are transformed to another 8 real numbers, which are real and imaginary parts of only 4 Fourier coefficients.

Cosine components of positive and negative frequencies are identical

Sine components of positive and negative frequencies have opposite signs

4cos(2 )

8t

3cos(2 )

8t )

8

42sin( t

3sin(2 )

8t

Forward/backward FFT

1 12 2

0 0

/2 /22 2

/2 /2

1 1

kn k nN Ni iN N

n k kk k

kn k nN Ni iN N

k n nn N n N

c c e c e

c c e c eN N

NAP3 Discrete FFT (Fast FT)

There is no difference between DFT and FFT. The FFT is very fast

algorithm which can be applied as soon as the number of point is a

power of two, e.g. N= 512,1024,2048,…

Note, that FFT is transformation of N-numbers ck to N-coefficients, corresponding to sampling interval =1s. You should take it into account when calculating convolution and correlation.

( ) ( ) /n n k kc f c c t c

fn

tk

dt

dtetcfc ift2)()(~

2( ) ( ) iftc t c f e df

N

df1

This screen ios an attempt to explain relationship between continuous and discrete FT. Notice the role of time step .

Parseval theorem

1

0

21

0

2 |~|1

||N

kk

N

kk c

Nc

NAP3 Discrete FFT (Fast FT)

dffcdttc 22 )(~)(

N

df1

dt kk cfc ~)(~( )k kc t c

In words: power of a signal calculated from the time course is the same as the power calculated from frequency domain. Why power? At electric signal is power at a constant resistance proportional to square of voltage.

NAP3 MATLAB examples of applications FFT

Braque

FFT functions in MATLAB

12 ( 1)

1 1

12 ( 1)

1 1

1 1(cos(2 ( 1) ) sin(2 ( 1) ))

1 1 1 1(cos(2 ( 1) ) sin(2 ( 1) ))

nN Ni kN

n k kk k

nN Ni kN

k n nn n

n nc c e c k i k

N N

n nc c e c k i k

N N N N

Forward transform of vector c

cf=fft(c,N)

Inverse transform

c=ifft(cf,N)

Vector of data c(1),c(2),….,c(N)

calculated Fourier coefficients cf(1),…,cf(N)

NAP3

NAP3

The following screen compares previously defin itiopn of FT (with positive and negative frequencies) with the definition used by MATLAB.

FFT functions in MATLAB

NAP3

( 1)2

1

k nN iN

n kk

c c e

( 1)( 1)

2

1

k nN iN

n kk

c c e

MATLAB (n=1 až N)FFT (n=-N/2 up to N/2 index expresses frequency)

It seems that MATLAB defines FFT in a different way and operates with frequencies which are even higher than the Nyquist frequency. Nevertheless the both expressions are correct and the calculated coefficients are identical (only sorted in a different way). This is a consequence of periodicity of goniometric functions (exp(2i(k-1))=1), and thus for example

1 2 3 /2 /2 1 /2 2

0 1 2 /2 1 /2 /2 /2 1

.... ....

.... = .... N N N N

N N N N

c c c c c c c

c c c c c c c

1

1 11 2 1 12 20 .... .... 2

c

N N

N N N N N

MATLAB

frequency=

f

index of coef. in MATLAB

Nyquist frequency

)12

(1

2)1(2)12

(1

2)12

(1

2

22/

~N

N

ki

k

kiN

N

ki

k

N

N

ki

kMATLAB ecececcN

FFT

previous definition FT

Definition in MATLAB

FFT functions in MATLAB

FFT MATLAB example PSDA typical application of FFT is identification of dominating frequencies, masked by a random noise. Let us con sider sampl;ed with frequency 1000 Hz. Let us generate signal with components at frequencies 50 Hz and 120 Hz superposed with random noise of zero mean value:

t = 0:0.001:0.6;x = sin(2*pi*50*t)+sin(2*pi*120*t);

y = x + 2*randn(size(t));

It is probably difficult to distinguish basic frequencies in time record of signal. Let us transform to frequency space using N=512 points and FFT:

Y = fft(y,512);

Power spectrum

Pyy = Y.* conj(Y) / 512;

It is sufficient to evaluate only one half of points (due to symmetry: y(t) is a real function and the Fourier coefficients of negative frequencies are conjugated with the coefficients of positive frequencies)

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7-6

-4

-2

0

2

4

6

Normal distribution of random numbers (mean value=0,

variance=1)

DFT using 512 points (power of two)

PSD Dominant frequencies

NAP3

FFT MATLAB example PSD

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

Dominant frequencies

NAP3

12 ( 1)

1

12

1

1( ( 1) )

n

nN i kN

n kk

ki t

Nn k

c c eN

c t n c eN

Horizontal axis in PSD graph are only indices of Fourier coefficients and it is seen that dominant frequencies correspond to the indices range 25 to 30 (first peak) and 60 to 65 (second peak). Corresponding frequencies follow from the definition of inverse transformation

N

kfk

1

k index in Fourier coefficients vector

starting from 1

For N=512 and =0.001s is therefore

f25=47 Hz, f30=57 Hz,

f60=115 Hz, f65=125 Hz

Nyquist frequencycase when only the k-th coefficients is

nonzero

f = 1000*(0:256)/512;

PSD noise filterThe simplest method of noise filtering is suppression of Fourier coefficients at high frequencies, for example components Y(65),Y(66),…Y(449) (assigned zero). Resulting PSD

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7-6

-4

-2

0

2

4

6

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5-4

-3

-2

-1

0

1

2

3

4

0 100 200 300 400 500 6000

10

20

30

40

50

60

70

80

Original signal for comparison

Vynulované fourierovy koeficienty

NAP3

n 514-nplease take into account that the last coefficient Y(512) does not correspond to zero but to lowest frequency.

Reconstructed signal y=ifft(Y,512)

FFT convolution/correlation

0 1 2 3 4 5 60

1

2

3

4

5

6

7

8x 10

-3

0 1 2 3 4 5 60

1

2

3

4

5

6

7

8x 10

-3

for i=1:512c1(i)=t(i)*exp(-2*t(i));c2(i)=t(i)^4*exp(-t(i)*3);endf1=fft(c1,512);f2=fft(c2,512);

for i=1:512c12(i)=f1(i)*f2(i); r12(i)=f1(i)*conj(f2(i));end

cc=ifft(c12,512)*0.001; rr=ifft(r12,512)*0.001;

Fourier coef. for convolution and correlation

effect of mirroring

(periodicity of FT)

NAP3

Mirroring effect is a consequence of DFT periodicity. The effect can be mitigated by using sufficiently long sampling interval.

at convulution and correlation is necessary multiply results by , at deconvoution the result

must be divided by

Inverse Fourier transform

FFT identification E(t)NAP3

The following screen illustrates basic techniques using radionuclide tracers for identification of impulse response of continuous systems.

Many specific applications in industry (fluidised combustion, furnace, reactors, abspní spalování, pece, reaktory, absorbers, columns, filters, aeration tanks, mills, containers) are discussed in monography Thýn J.,Žitný R.:Analysis and diagnostics of industrial processes by radiotracers. ČVUT Praha 2000

FFT identification E(t)NAP3

dt=0.01;t=0:dt:10.23;n=2;tm=1;x=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm);n=6;tm=3;y=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm);x=x+0.01*randn(size(t));y=y+0.005*randn(size(t));

(t)y(t)

x(t)

Let us consider a system consisting in 4 serially connected mixed tanks. The impulse response E(t) of this series can be expressed in an analytical form (for example by using Laplace transformation but this is not important now)

mtNtNm

NN

etN

tNtE /

1

)!1()(

N-number of mixed tanks, tm is mean residence time (volume/flowrate)

At inlet the signal x(t) and at outlet the signal y(t) (for example concentration of a tracer) are known, for example measured. These signals will be created artificially: x(t) as a response of system consisting in 2 mixers (N=2). Then the output y(t) must be impulse response of 6 identical mixers (N=2+4) with mean time corresponding to the sum of mean times tm of input and identified system. The generated signals x(t) and y(t) will be distorted by a random noise

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

x(t)

y(t)

the same can written in MATLABu shortly as

e=@(t,tm,n) n^n.*t.^(n-1)./(factorial(n-1)*tm^n).*exp(-n.*t/tm)t=linspace(0,10.23,1024);x=e(t,1,2)+0.01*randn(size(t));y=e(t,3,6)+0.005*randn(size(t));

anonymous function @ with 3 parameters. Linspace is a

standard function.

FFT identification E(t)NAP3

xf=fft(x);yf=fft(y);ef=yf./xf;e=ifft(ef)/dt;plot(t,e)

Impulse response of identified system can be obtained by deconvolution, by application of FFT to input and output (vectors of Fourier coefficients xf(1:1024), yf(1:1024)), Fourier coefficients of function E(t) are calculateds as ration of coefficients yf(i)/xf(i) and impulse response is obtained by inverse FFT

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

xf=fft(x);yf=fft(y);ef=yf./xf;ef(17:1024-16)=0;e=ifft(ef)/dt;

0 2 4 6 8 10 12-15

-10

-5

0

5

10

15

The result is a wildly oscillating function having little in common with the theoretical impulse response (blue curve n=4; tm=2; et=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm); ). This is because deconvolution is a badly conditioned problem when even a small deviation of the output signal y(t) will cause a large deviation of solution. Noise is possible to suppress by resetting higher frequencies, eg

theoretical impulse response

( ) ( ) / ( )E f y f x f

FFT identification E(t) WienerNAP3

In the previous case, we filter the noise from the resulting impulse response. Maybe a little more correct is to filter noise just from the output function y(t) by suppressing the higher frequency coefficients in the Fourier transform yf using Wiener filtration (its theoretical description can be found eg in Teukolski et al.: Numerical Recipes). Again, procedure is based on the graph PSD of the noisy response, from which the frequency of "useful" signal and the frequency which is to be suppressed are estimated.

0 10 20 30 40 50 60 70 80 90 100-15

-10

-5

0

5

10

2

2

|)(~|

|)(~|)(

~fy

fyf c

PSD „noised“ signal

estimated PSD of a „net“ signal

by this correction function is multiplied Fourier transform of

measured output

log PSD only first 100 frequencies PSD are shown. The red dotted line is a noise extrapolation separating useful signal and noise.

FFT signal time delayNAP3

Transformations of “time shifted" signal is realized only by multiplying FT by exponential. This facilitates the modeling of systems in which, for example trasport delay occurs (which causes problems for the numerical solution of differential equations).

0 2 4 6 8 10 12-0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

dt=0.01;t=0:dt:10.23;for i=1:1024 f(i)=(i-1)/(1024*dt);end% vstup tm=1, N=2n=2;tm=1;x=n^n.*t.^(n-1)/(factorial(n-1)*tm^n).*exp(-n.*t/tm);xf=fft(x);tau=-2;xfp=xf.*exp(2*pi*tau*complex(0,1).*f);yf=xfp.*xf;xp=ifft(xfp);y=ifft(yf)*dt;

(t)y(t)x(t) xp(t)

)(~)()()(~ 2222 fcedtetcedtetcfc ififtififtp

FT of function delayed by

Response of 2 mixers

Piston flow (delayí )

x(t)xp(t)

y(t)

Vector of frequencies

Shift by tau

Correlation of stimulated or random signal from two locations (technically for example correlation of two temperatures - thermocouples)

dttTtTR )()()( 2112

T1 T2

Heater (source of random pulses)

NAP3 FFT mutual correlation

A time shift of signals and from it for example velocity, flowrate,…can be evaluated from correlation function. But also necessary deformation of the surface load components captured using a digital camera (image processing)

dttTtTR )()()( 2112

Example simulated in MATLAB

Heater

Random signal shifted by 100 time steps

NAP3 FFT cross correlation

location of this peak determines transit time and therefore flowrate

EXAMNAP 3

Fourier transform

NAP3 Exam - remember

dxtEty )()()(What is it convolution, deconvolution, transfer function and correlation

dtetcfc ift2)()(~ 2( ) ( ) iftc t c f e df

Fourier transform

Fourier transform convolution, correlation ( ) ( ) ( )xyC f x f y f )(~)(~)(~ * fyfxfRxy

Discrete Fourier transform1 /22 2

0 /2

1

kn knN Ni iN N

n k k nk n N

c c e c c eN

Nyquist frequency 1/2

( ) ( ) ( )xyR t x t y d