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

1. Experimental methods E181101 EXM9 Time seriesNoise filtrationFourier transformWaveletsStimulus-response techniques Rudolf Žitný, Ústav procesní a zpracovatelské techniky ČVUT FS 2010

2. Median (warm up - introductory exercise noise reduction of a signal y1, y2,….yN) EXM9 Median N=9 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 Nb=50 total number of points N=1024

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). Regression smoothing Nb=50 total number of points N=1024 Median smoothing Nb=50 total number of points N=1024

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 n=1000; for i=1:n t(i)=i/n; e(i)=exp(-t(i))*sin(20*t(i)); end r=rand(n,1); er=e+r'; plot(t,er); Moving average filter w=10; for i=1:n wl=max(1,i-w); wr=min(n,i+w); ef(i)=mean(er(wl:wr)); end plot(t,ef)

7. Smoothing (median) EXM9 for i=2:n-2 wl=max(1,i-w); wr=min(n,i+w); nvalues=wr-wl+1; median=floor(nvalues/2)+1; for j=wl:wr ej=er(j); ngt=0; for k=wl:wr if er(k) >= ej ngt=ngt+1; end end nlt=nvalues-ngt; if (ngt >= median) & (nlt <= median) ef(i)=ej; break end end end

8. System responsestimulus response technique. EXM9 Balthus

9. cin cout E t [s] 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 t=t2-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 just info 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 Přerov) Radionuclide mixed with processed material Na

12. Example: rotary dryer- RTD just info 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 Cocurrent D Counter-current 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 M L G gas gas just info Cocurrent D Counter-current particle particle 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 d n k k t [mm] [min/m] [min] m 1 2 2.24 1.274 1.692 0.0037 47.2 5.70 3.750 1.830 0.0150 48.3 A * just info Dehydration D3 t / t 0 DN 0 . 9 / D N 0 . 9 0 0 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. E t [s] 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. 0  t[s] 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 t=0. ? Constant flowrate at inlet 0  t[s] Concentration in tank is the same as at outlet (perfect mikxing, concentration is uniform inside) Constant flowrate at outlet

18. 0  t[s] 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. 0  t[s] 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 t> 0  t[s] Response valid for t > 

21. 0  t[s] System models stimulus response technique. EXM9 Response to puls of width =2s and =1 s

22. x(t) x(t) 1/ E(t-) E(t)   t [s] t [s] 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. x(t) y(t) E(t) x y E t [s] 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 m=n, otherwise 0 Proof!!! In a similar way the orthogonality of sin nx can be derived. From Euler’s formula follows orthogonality of i-imaginary unit Linear combination of Pj(x) is called Fourier’s expansion The transformation T(x) to Tj for j=0,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. Heater Fourier transforms EXM9 Example calculated by MATLAB Random signal shifted by 100 time steps

33. 2 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 cf=fft(c,N) Vector of calculated Fourier coefficients) Vector of samples c(1),c(2),….,c(N) Inverse transformation c=ifft(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 = 0:0.001:0.6; x = sin(2*pi*50*t)+sin(2*pi*120*t); y = x + 2*randn(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*(0:256)/512; Normal distribution (mean=0, variance=1) 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 y=ifft(Y,512)

40. Fourier transform convolution/correlat. EXM9 for i=1:512 c1(i)=t(i)*exp(-2*t(i)); c2(i)=t(i)^4*exp(-t(i)*3); end f1=fft(c1,512); f2=fft(c2,512); Fourier coef. of convolution and correlation for i=1:512 c12(i)=f1(i)*f2(i); r12(i)=f1(i)*conj(f2(i)); end cc=ifft(c12,512); rr=ifft(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áček, 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áček, 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) 98–124 R.C. Guido / Applied Mathematics Letters 24 (2011) 1257–1259 EXM9 Deconvolution Wavelet blurred Tichonov regularization Wiener filter

47. Wavelet (papers) EXM9