Proveden

1
Experimental methods E181101
EXM9
Time seriesNoise filtrationFourier
transformWaveletsStimulus-response techniques
Rudolf Žitný, Ústav procesní a zpracovatelské
techniky CVUT FS 2010
2
Median (warm up - introductory exercise noise
reduction of a signal y1, y2,.yN)
EXM9
Median N9 5 3 1 7 9 10 1 -4 4
-4 1 1 3 4 5 7 9 10
Sort the numbers according to values
Median is in the middle 4
3
Smoothing (noise reduction of a signal y1,
y2,.yN)
EXM9
Median smoothing The main idea is to run through
the signal replacing each entry with the median
of neighboring entries. The pattern of
neighbors is called the "window", which slides,
entry by entry, over the entire signal. For 1D
signals, the most obvious window is just the
first few preceding and following entries,
whereas for 2D (or higher-dimensional) signals
such as images, more complex window patterns are
possible (such as "box" or "cross" patterns).
Note that if the window has an odd number of
entries, then the median is the middle value
after all the entries in the window are sorted
numerically.
Median smoothing Nb50 total number of points
N1024
4
Smoothing (noise reduction of a signal y1,
y2,.yN)
EXM9
Regression smoothing Savitzky Golay Instead of
median the signal inside the window is
approximated by a polynomial linear regression
(degree of polynomial is obviously limited by the
width of window).
Median smoothing Nb50 total number of points
N1024
Regression smoothing Nb50 total number of
points N1024
5
Smoothing (noise reduction of a signal y1,
y2,.yN)
EXM9
Moving average and median filter Midpoint in the
moving window is replaced by average value. A
moving average may also use unequal weights for
each data value in the window to emphasize
particular values in the subset.
Write a MATLAB program for moving average and
median filter as an excercise
6
Smoothing (noise reduction of a signal y1,
y2,.yN)
EXM9
n1000 for i1n t(i)i/n e(i)exp(-t(i))sin(20
t(i)) end rrand(n,1) erer' plot(t,er)
Moving average filter
w10 for i1n wlmax(1,i-w) wrmin(n,iw) ef(i
)mean(er(wlwr)) end plot(t,ef)
7
Smoothing (median)
EXM9
for i2n-2 wlmax(1,i-w) wrmin(n,iw)
nvalueswr-wl1 medianfloor(nvalues/2)1
for jwlwr ejer(j) ngt0 for
kwlwr if er(k) gt ej
ngtngt1 end end
nltnvalues-ngt if (ngt gt median) (nlt lt
median) ef(i)ej break end
end end
8
System responsestimulus response technique.
EXM9
Balthus
9
System responsestimulus response technique.
EXM9
Some characteristics of identified continuous
system (a reactor, furnace, heat exchanger,
almost anything) can be obtained by monitoring
time response of system y(t) to a stimulus
function x(t). Stimulus can be represented by a
marker injected to the inlet of system and
response its concentration detected at outlet. As
marker (indicator) are frequently used salts
(concentration is measured by conductivity
meter), dyes (detected by colorimeter), bubbles
or fine particles (detected by ultrasound or
laser), or radionuclides (detected by
gamma-radiation detectors).
Dispersion in pipe
10
System responsestimulus response technique.
EXM9
Applications experimental characterization of
continuous systems

• Look at the database using keywords residence
  time distribution, transfer function, Peclet
  number, axial dispersion. Some of the
  following applications can be found on my web
  pages
• Integral characteristics (moments of continuously
  recorded signals at inlet and outlet of
  apparatuses). Example Packed columns evaluation
  of holdup (mass of flowing film as a function of
  flowrate). Holdup calculated from first moments
  tt2-t1, dispersion (Péclet number) from second
  central moments.
• Residence time distribution RTD (flight times of
  particles flowing through an apparatus). Example
  RTD identification of a polydisperse mixture
  (titanium oxide) in a horizontal rotary drum
  using radioisotope tracer Na24 mixed with the
  processed material at inlet (see A in figure).
  Identified RTD enables simulation of drying and
  calcination reactions inside a drum.
  This kind of analysis is typical for chemical
  reactors, combustors, waste water treatment, food
  processing lines.

11
Example rotary dryer- RTD
EXM9
Experimental analysis using radionuclide tracer
(Na-half life 15 hours) mixed with the TiO2 feed.
Concentration of tracer was recorded by 12 gamma
radiation detectors mounted along the rotating
cylinder (from outside, gamma radiation
penetrates steel wall). This is the way how to
track the transported particles and how to
identify their residence times.
Example Residence Time Distribution (rotary
dryer/kiln TiO2 Precheza Prerov)
Radionuclide mixed with processed material
Na
12
Example rotary dryer- RTD
EXM9
Theory
M
It was theoretically derived that the residence
time of a spherical particle (diameter d) moving
in an inclined rotary dryer (diameter of cylinder
D, length L, inclination angle ?, frequency of
rotation N) is in the counter current arrangement
(when gas flows against transported particles)
and in the parallel (cocurrent) flow (MG
and MS is mass flow rate of gas and solid
repectively). It is seen that the residence time
decreases with the increasing frequency of
rotation, diameter and inclination of cylinder.
Influence of particle size depends upon mass
flowrate of gas. Friedman F.,Marshall
A.,Chem.Engng.Progr. 45,482, 1949
L
G
gas
gas


Counter-current
Cocurrent
D
particle
particle
Short residence time in the co-current flow
Long residence time in the counter-current flow
Values t are residence times of identical
particles (monodisperse mixture)
13
Example rotary dryer- RTD
EXM9
Theory
For polydisperse material Rosin Rammler
distribution of particle size must be identified
(dm is the mean diameter of particle, n- is
characteristic of dispersion)
Residence time distribution of particles in
counter-current arrangement
Co-current arrangement
Mean residence time for countercurrent arrangement
Mean residence time for cocurrent arrangement
These functions are distributions of residence
times for polydisperse material (Rosin Rammler
distribution of particle size)
This is only an approximation at co-current flow
(fine particles would have negative residence
time according to this simplified model)
14
Example rotary dryer- RTD
HP9
Experimentally determined impulse response of a
dryer section (by recording response of a gamma
detector to the instantaneous injection of Na
tracer) was used for identification of k1 and k2
parameters of previously derived RTD models.
The values k1 and k2 enable to predict what
happens when the rotational speed N is changed
(and it was the primary goal of analysis to
answer the question how the frequency of rotation
affects lengths of drying and reaction zones in
the TiO2 kiln it is not so easy to change the
frequency, it is not sufficient just only to turn
a knob, such a change needs money and must be
supported by preliminary analysis).

2
1.5
1
0.50
0.75
0.5
1
1.25
1.50
0
0 0.5 1
1.5 2 2.5

15
System responsestimulus response technique.
EXM9
By measuring the mean residence time it is
possible to calculate for example yield of a
chemical reactor (yield depends upon the time
available for reaction)
It tells you something about uniformity of
residence times
16
System models stimulus response technique.
EXM9
Compartment models
Compartment is a basic unit (like a brick in
LEGO) for description complicated systems
structure. Compartment is an ideally mixed tank,
characterfized by mass, temperature,
concentration
17
System models stimulus response technique.
EXM9
Derive response of a perfectly mixed tank of
inner volume V to a short pulse
Initial condition zero concentration c in the
tank at time t0.
?
Constant flowrate at inlet
0 ? ts
Concentration in tank is the same as at outlet
(perfect mikxing, concentration is uniform inside)
Constant flowrate at outlet
18
System models stimulus response technique.
EXM9
Procedure First step describe the system by
differential equation for c(t)
Mass balance of component having concentration c
Mean residence time
19
System models stimulus response technique.
EXM9
Intermediate step calculate response to a unit
step stimulus function
Response valid up to the time ?
20
System models stimulus response technique.
EXM9
Intermediate step response to zero concentration
at inlet for tgt?
0 ? ts
Response valid for t gt ?
21
System models stimulus response technique.
EXM9
Response to puls of width ?2s and ?1 s
22
IMPULSE response E(t)
EXM9
Impulse response E(t) is the response of a
general system to infinitely short pulse (??0)
having unit area (delta function).
Delta function ?(t) is zero for t?0 and
23
System response to stimulus function.
EXM9
Convolution integral
Response y(t) follows from the principle of
superposition of responses to narrow pulses
(representing stimulus x(t))
x
y
E
t s
24
System response to stimulus function.
EXM9
Typical problems Given a stimulus function x(t),
the prediction of system response y(t) is
calculated from the convolution integral Given
measured stimulus and response functions, the
impulse response E(t) is evaluated by
identification from the Volterra equation
(deconvolution)
This is difficult problem because identification
is sensitive to signal noise
These problems are solved using integral
transforms, see
FOURIER transform?
25
FOURIER transform preliminaries
EXM9
Duchamp
26
FOURIER transform preliminaries
EXM9
Goniometric functions (sin, cos) are orthogonal
in interval (-?,?). Orthogonality of Pn(x)cos
nx follows from
for mn, otherwise 0
Proof!!!
In a similar way the orthogonality of sin nx can
be derived. From Eulers formula follows
orthogonality of
i-imaginary unit
Linear combination of Pj(x) is called Fouriers
expansion
The transformation T(x) to Tj for j0,1,2,, is
called Fourier transform and its discrete form is
DFT T(x1), T(x2),. T(xN) to T1,T2,TN . DFT can
be realized by FFT (Fast Fourier Transform).
27
Fourier transforms
EXM9
Continuous Forward Fourier transform from time to
frequency domain
Backward transform
28
Fourier transforms
EXM9
Basic properties of FT
29
Fourier transforms
EXM9
power spectral density
Energy of high frequency (noise)
Energy of long waves (low frequency)
30
Fourier transforms
EXM9
convolution
correlation

Mean time of x(t)0.5
Mean time of y(t)2.1

Mean time of cross correlation Rxy(t)1.6
represents a time shift between x(t) and y(t)
31
Fourier transforms
EXM9
Cross-correlation of stimulated or random signal
detected at two locations (technically it can be
a heater and thermocouples)
Heater
T1
T2
32
Fourier transforms
EXM9
Example calculated by MATLAB
Random signal shifted by 100 time steps
33
Discrete Fourier transform
EXM9
Sampling of data at constant time step ? N-data
points (even number)
Nyquist frequency
Nyquist frequency is the maximum frequency which
can be described by data with the sampling
interval ?
34
Discrete Fourier transform
EXM9
only the sum is called DFT
Discrete FFT has the same properties
(convolution, correlation) as the continuous FT.
35
Discrete Fourier transform
EXM9
Discrete FT



Parseval theorem
36
Wiener filtering
EXM9
s(t)
r(t)-impulse response
u(t)
corrupted c(t)
noise n(t)
37
Discrete Fourier transform
EXM9
MATLAB programming Forward transformation cffft
(c,N)
Vector of calculated Fourier coefficients)
Vector of samples c(1),c(2),.,c(N)
Inverse transformation cifft(cf,N)
38
Fourier transform PSD example
EXM9
A common use of Fourier transforms is to find the
frequency components of a signal buried in a
noisy time domain signal. Consider data sampled
at 1000 Hz. Form a signal containing 50 Hz and
120 Hz and corrupt it with some zero-mean random
noise t 00.0010.6 x sin(2pi50t)sin(2
pi120t) y x 2randn(size(t)) It is
difficult to identify the frequency components by
looking at the original signal. Converting to the
frequency domain, the discrete Fourier transform
of the noisy signal y is found by taking the
512-point fast Fourier transform (FFT) Y
fft(y,512) The power spectrum, a measurement of
the power at various frequencies, is Pyy Y.
conj(Y) / 512 Graph the first 257 points (the
other 255 points are redundant) f
1000(0256)/512
Normal distribution (mean0, variance1) of
random numbers (result is of the same size as t)
Discrete Fourier Transform of 512 points
Dominant frequencies
Power spectral density
the upper half is not important because x is real
and not a complex vector
39
Fourier transform PSD filter
EXM9
The simplest way how to filter a noise is to
suppress high frequencies, for example all
frequencies corresponding to fourier components
Y(65),Y(66),Y(451) (assign zeroes to these
entries). Resulting PSD is
Annulated fourier coef.
Original signal for comparison
Filtered signal is reconstructed by inverse
FT yifft(Y,512)
40
Fourier transform convolution/correlat.
EXM9
for i1512 c1(i)t(i)exp(-2t(i)) c2(i)t(i)4
exp(-t(i)3) end f1fft(c1,512) f2fft(c2,512)
Fourier coef. of convolution and correlation
for i1512 c12(i)f1(i)f2(i)
r12(i)f1(i)conj(f2(i)) end ccifft(c12,512)
rrifft(r12,512)
effect of mirroring (periodicity of FT)
Inverse Fourier transformation
41
Wavelets
EXM9
Dalí
42
Wavelets
EXM9
Integral transform with scale k and time shift t
scale
translation
Example of mother wavelet g(z) is Mexican hat
scale k
Time shift t
43
EEG and Visual information processing
EXM9
Information direction from eyes to brain cortex
http//people.deas.harvard.edu/gstanley/research_
topics_vision.html
44
EEG and Visual information processing
EXM9
M. Hulan, J. Kremlácek, R. Žitný New methods for
automatic detection of evoked potentials of a
human primary visual cortex. HBM 2006
Significant VEP peak manually marked
Reversing checker stimulus
U mV
Top of peak - cca 100 130 ms after stimulation
nasion
6 electrodes system
Signal averaging from 40 realizations
Referential 7th ear electrode
t 1 sample 2ms
inion

• Recorded electric activity of primary visual
  cortex related to the visual stimulus
• -Signal is averaged from many realizations to
  increase signal-to-noise ratio
• We can measure differences in the el. activity of
  the primary visual cortex related to the
  stimulus.
• -For clinical evaluation of neural diseases is
  used record from area near the visual cortex

45
EEG and Visual information processing
EXM9
M. Hulan, J. Kremlácek, R. Žitný New methods for
automatic detection of evoked potentials of a
human primary visual cortex. HBM 2006
VEP-Visually Evoked Potential CWT-Continuous
Wavelet Transformation
46
Wavelet (papers)
A.Z. Averbuch et al. / Appl. Comput. Harmon.
Anal. 31 (2011) 98124 R.C. Guido / Applied
Mathematics Letters 24 (2011) 12571259
EXM9
Deconvolution
Wavelet
blurred
Tichonov regularization
Wiener filter
47
Wavelet (papers)
EXM9