- About this Journal
- Abstracting and Indexing
- Aims and Scope
- Article Processing Charges
- Articles in Press
- Author Guidelines
- Bibliographic Information
- Citations to this Journal
- Contact Information
- Editorial Board
- Editorial Workflow
- Free eTOC Alerts
- Publication Ethics
- Reviewers Acknowledgment
- Submit a Manuscript
- Subscription Information
- Table of Contents
Modelling and Simulation in Engineering
Volume 2012 (2012), Article ID 362849, 9 pages
Basic Characteristics of IEC Flickermeter Processing
Poznań University of Technology, Poznań, Poland
Received 5 April 2012; Revised 15 October 2012; Accepted 30 November 2012
Academic Editor: MuDer Jeng
Copyright © 2012 Jarosław Majchrzak and Grzegorz Wiczyński. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Flickermeter is a common name for a system that measures the obnoxiousness of flicker caused by voltage fluctuations. The output of flickermeter is a value of short-term flicker severity indicator, . This paper presents the results of the numerical simulations that reconstruct the processing of flickermeter in frequency domain. With the use of standard test signals, the characteristics of flickermeter were determined for the case of amplitude modulation of input signal, frequency modulation of input signal, and for input signal with interharmonic component. For the needs of simulative research, elements of standard IEC flickermeter signal chain as well as test signal source and tools for acquisition, archiving, and presentation of the obtained results were modeled. The results were presented with a set of charts, and the specific fragments of the charts were pointed out and commented on. Some examples of the influence of input signal’s bandwidth limitation on the flickermeter measurement result were presented for the case of AM and FM modulation. In addition, the diagrams that enable the evaluation of flickermeter’s linearity were also presented.
Flickermeter is a common name for a system that measures the obnoxiousness of flicker caused by voltage fluctuations. According to , flickermeter is an instrument designed to measure any quantity representative of flicker. Flicker is an impression of unsteadiness of visual sensation induced by a light stimulus whose luminance or spectral distribution fluctuates with time . According to , flickermeter is a flicker measuring apparatus to indicate the correct flicker perception level for all practical voltage fluctuation waveforms. The processing performed by this system’s signal chain is complicated  to such an extent that it is not easy to obtain output values only analytically. Bibliography concerning this subject contains numerous works that describe the operation and features of flickermeter [5–10]. In most of them, authors reconstruct flickermeter’s characteristics for the case of amplitude modulation (AM) [3, 11, 12], input voltage with interharmonic component [12–21], or step changes of input voltage phase [12, 22–25]. The bibliography lacks publication on reconstructing the characteristics of flickermeter processing in a systematic and complex manner. Many discussions and comments contain opinions that vaguely describe the way in which flickermeter actually operates. For example, one of the most common mistakes is the erroneous limitation of frequency range for which the characteristics of IEC flickermeter are constructed. This appears to come from wrong conclusion that modulation with frequencies higher than 40 Hz does not cause obnoxious flicker because the only flicker obnoxious for human is limited to the frequency of 40 Hz. Nevertheless, because of the way incandescent light sources function (which is mapped in an IEC flickermeter block diagram with operation of voltage squaring), flicker could be sensed for frequencies of modulating signal higher than 40 Hz.
Herein, the results of numerical simulations that reconstructed the processing of signal chain of IEC flickermeter will be presented. The research results presented here are concerned with the processing of IEC flickermeter; in other words, the flickermeter is built accordingly to IEC 61000-4-15, and thereby does not describe the flicker phenomenon in general, for instance, for light source other than the one defined in the standard. Therefore any modification of the test signal path is abandoned. The novelty is a comprehensive study of flickermeter for various voltage variation.
Owing to space constraint, this paper contains only selected simulation results that are important from the analytical point of view or have some practical significance .
2. Test Signals for Flickermeter Signal Chain
The dynamic properties of measuring devices are determined in frequency domain. The results of the research on frequency domain are amplitude characteristics and, if applicable, phase characteristics. The above-mentioned means of determination of dynamic properties form the classic approach to the subject. Therefore, an attempt has been made to use them to test the flickermeter signal chain while taking into account the occurrence of carrier component in the input voltage. At the same time, an assumption has been made that the tests will be carried out in a steady state of the signal chain, that is, after the transient component has completely faded out. This indicates that, for example, step change of voltage will appear after the signal (with some preset amplitude ) has been applied previously for a time needed to fade out transient component completely. One of the benefits of carrying out the analysis of signal chain in time domain, besides enriching the set of reference signals, is the possibility to obtain the values of output and inner signals during the transient state (which is helpful when trying to determine if any saturation occurs) and to determine how long it takes for the transient component to fade out.
An ideal input signal in steady state could be expressed as: where is the voltage amplitude, is the carrier wave frequency (corresponding to the time period ), is the pulsatance of carrier wave, and is time.
The test in the frequency domain will utilize the following input signals : (i)AM modulated with not suppressed carrier wave where is the modulating signal that satisfies condition and is the modulation depth; (ii)with single interharmonic component where is the interharmonic component amplitude and is the interharmonic component pulsatance; (iii)FM modulated  where is the carrier frequency modulation depth scaling coefficient and is the modulating signal, frequency deviation
3. IEC Flickermeter Simulator Structure
The system is composed of two main parts (Figure 1): flickermeter signal chain modeled according to the standard specification and simulation support block. The task of simulation support block is to generate test signals with the chosen time plots, amplitudes, and frequencies to determine the rms value () of signal, to determine the value of signal, as well as to record and visualize selected signals. The structure of the test signals source enables generation of signal as a result of AM modulation, FM modulation, or as a sum of carrier and interharmonic components. Low-pass Butterworth-type filter of 6-order with a cut-off frequency limits the bandwidth of signal.
The AGC block replaces an input transformer with branches taps and input voltage conditioning circuitry. The parameters of RMS/LPF filter were chosen so as to achieve raising/falling time of signal equal to 1 min for the step change of rms value. By assuming Hz as a frequency of voltage in power network, filters from Figure 1 could be specified as follows: Hz, Hz, Hz and Hz (the detailed specification in ). IFL is a signal that appears on output 5 in standard flickermeter signal chain schema . The value of signal corresponds to the maximal value of IFL signal. Statistical analysis block calculates the value of short-term flicker severity indicator, , on the basis of statistical distribution of IFL signal.
All instantaneous values of , , and other signals are recorded by the recording and visualization block. The process model of flickermeter and simulation researches were carried out by using the software application Matlab and Matlab-Simulink [29, 30], using the solver ODE45 with variable step (Dormand-Prince).
4. Basic Characteristics of Flickermeter Processing
4.1. Specification of Flickermeter Basic Characteristics
The following groups of characteristics were chosen to describe the processing of flickermeter signal chain: (1)time plots of internal signals and IFL during of fading of intial transient component, (2)grouping of output signal values for no-modulation state (signal in accordance with (1)), (3)dependence of output signal and internal signal on frequency , modulating signal type, and modulation depth (amplitude modulation—signal in accordance with (2)), (4)dependence of the modulation depth on frequency and modulating signal type for the preset values of output signal or internal signal (amplitude modulation—signal in accordance with (2)), (5)dependence of output signal on interharmonic component frequency, (signal with interharmonic component—signal in accordance with (3)), (6)dependence of interharmonic component amplitude on frequency for preset values of output signal (signal with interharmonic component—signal in accordance with (3)), (7)dependence of output signal on frequency and frequency deviation (frequency modulation—signal in accordance with (4)), and(8)dependence of frequency deviation on sinusoidal modulating signal for the preset values of output signal (FM modulation—signal in accordance with (4)).
Characteristic groups, 1–3, 5, and 7, were obtained by setting appropriate properties of the input signal while recording the values of the respective signals. The characteristic group denoted as 4 differs from the other groups in the way it was obtained. In that case, an iterative algorithm for determination of modulation depth was used, which is presented in Figure 2. In the standard 61000-4-15 , the action of the flickermeter for rectangular and sine signal was defined, which modulates the amplitude with a frequency of up to 33 Hz. However, the flicker phenomenon also occurs for higher frequencies. Therefore, this study has presented the characteristics of the full range of frequencies that is likely to cause noticeable flicker. This applies to the modulation characteristics of both AM and FM. As the initial-phase modulating signal did not affect the rate of , all simulations were performed for the initial phase equal to zero. It is worth to point out the time-consuming aspect of the simulation process. Observation time was set to to satisfy the standard requirements. To eliminate the influence of the transient state of signal chain on the simulation result, of extra time was added before each start of the standard observation period. With the computer used for calculations, determination of the simulations results took significantly more time. It took several dozens of hours to obtain multipoint characteristic, especially when iterative algorithm was utilized. Characteristic groups 6 and 8 were build with the use of algorithm presented in Figure 2, while substituting the modulation depth with interharmonic component amplitude and frequency deviation , respectively.
4.2. Measurement Result for Sinusoidal Input Signal (without Modulation)
When a steady sinusoidal signal without modulation (in accordance with (1)) is applied on flickermeter input, the output signal should remain at zero. However, in case of standard flickermeter signal chain , the value for such an input signal is greater than zero. Table 1 summarizes the values of signal obtained with simulations for different orders of filter.
The test for the absence of voltage variation indicates that the real flickermeter shall not be less than approx. 0.01 (if the meter shows , it can be inferred that in the processing the “trick” was used).
4.3. Fading of Initial Transient Component
By assuming zero initial conditions, when the signal, in accordance with (1) is applied, some transient component appears. The fading time of that component is determined using the features of flickermeter signal chain. Figure 3 presents the fading out of the transient component at the two distinctive points of the signal chain: signal at the output of AGC block and IFL signal at the output 5 (see diagrams in Figure 3).
Analysis of Figure 3 leads to the following conclusions: (i)the maximum value of signal is about times greater than the value at the steady state, (ii)the maximum value that the IFL signal takes on is about , while the value at the steady state is zero, and(iii)the fading out time of signal transient component is about 180 s, and the fading out time of IFL signal transient component is about 120 s.
One of the effects of such a large short-lasting value of IFL signal may be the occurrence of saturation of flickermeter signal chain. The occurrence of such state may have an effect on additional error of indicator measurement, which could be hard to estimate. Such erroneous state occurs when the fluctuation of voltage is sufficiently strong and repeats over time. An example of such conditions is the fluctuation of voltage in power circuit of arc furnace.
4.4. Flickermeter Processing Characteristics for AM Modulation of Input Signal
The characteristics of flickermeter processing for AM-modulated input signals, described with (2), can be divided into two groups. The first group contains characteristics for a preset modulation depth, . The second group contains characteristics for and . Figure 4 presents the graph of dependence for sinusoidal, triangular, and rectangular modulating signal with constant value. The modulation depth was set to the value that guarantees unitary indicator value at frequency Hz, that is, for maximum sensing of flicker. Flickermeter signal chain interrelates the output value with the parameters of input signal: modulation depth , frequency , and the shape of the signal. To determine characteristics, the iterative procedure was used, in which the shape of the modulating signal was set along with and values. Figures 5 and 6 present the sets of characteristics for three modulating signals: rectangular, sinusoidal, and triangular with and , respectively.
Verification of comparison of the selected characteristics with points specified in the standard IEC61000-4-15  was carried out. Normative characteristics of the flickermeter’s processing specification are subject only to amplitude modulation by using signals with frequency of up to 33 Hz. The credibility of the simulation for higher frequencies has already been evidenced by the observation of flickering lights and by comparing with the results of model tests. Comparison of the determined characteristics leads to the conclusion that the most “obnoxious” modulation signal is the rectangular one, and the least “obnoxious” is the triangular one. Similar to the plot presented in Figure 4, three local minima are observed to exist, with the distinction that the minimal value of the modulation depth occurs for the frequency Hz. The other extrema occur for Hz and Hz. Relation for the case of modulation with rectangular signal cab be distinguished by peculiar non-monotonicity in a frequency range of 25 Hz–40 Hz.
The two groups of characteristics, and are complemented with characteristic, which was determined for modulation with rectangular signal, as presented in Figure 7. This enables verification of flickermeter linearity, while assuming the modulation depth as an input quantity and as an output. Accordingly, it is possible to state that, in general, flickermeter is not a linear system, but for inputs that correspond to (which means sensing of flicker), this system is nearly linear.
Figure 8 presents the plot of for AM modulation with rectangular signal. It gives the information on how the bandwidth of the input signal reflects in the flickermeter measurement result. The influence of bandwidth limitation becomes visible for Hz.
4.5. Flickermeter Processing Characteristics for Input Signal with Single Interharmonic Component
An input signal with single interharmonic component used to determine the characteristics of flickermeter processing is defined using (3). Figure 9 presents the characteristic constructed with the use of the algorithm presented in Figure 2 while exchanging a modulation depth with amplitude of interharmonic . For the sake of comparison, the characteristic for AM modulation with sinusoidal signal is also presented. The main difference between the two is the value of frequency, for which the local minimum occurs. For the input signal with interharmonic component, the maximum flicker sensing occurs at Hz and Hz. It is worth noting the difference between the signal with single interharmonic component and the AM-modulated signal, defined with (2). Amplitude-modulated signal contains at least two interharmonics: in the case of modulation with sinusoidal signal, it contains two interharmonics, and in the case of modulation with deformed signals, it may contain, theoretically, infinite number of interharmonics.
4.6. Flickermeter Processing Characteristics for FM Modulation of Input Signal
The characteristics of flickermeter processing for input signals obtained as a result of FM modulation were obtained with signal defined as (4). Frequency modulation is a nonlinear operation. This fact complicates the reproduction of flickermeter processing characteristic, because the input signal must be specified in a way that takes into account the working point of the signal chain. Figure 10 presents a dependency for frequency deviation of 0.05 Hz, 0.5 Hz, and 2 Hz, and modulation with rectangular signals. The maximum of characteristics occur at Hz. For a case of modulation with rectangular signal, the dependency is highly non-monotonic in a Hz frequency range.
The result of FM modulation is usually a broadband signal. According to , “the pass bandwidth of input stage … should not introduce an extensive suppression at least up to 700 Hz.” Figure 11 presents the characteristic obtained for the changing values of frequency. The bandwidth of the input signal was limited with a low-pass LPF filter with a cut-off frequency adjusted in the range of 50–800 Hz. Surprisingly, to some extent, the resulting characteristic shows that the limitation of the bandwidth leads to increased value of output signal for almost all of the preset values of frequency .
The influence of bandwidth limitation on measurement result depends on frequency value. Figure 12 presents the characteristic for the case of FM with rectangular signal. The increase in the value for Hz can be clearly observed.
Figure 13 presents the characteristic for modulation with rectangular signal. The evaluation of characteristic linearity is complex for both the cases. For some values of frequency (i.e., 109 Hz, 91 Hz, and 78 Hz), it could be treated as linear, for other values, it is nonlinear, while for the smallest values, it is non-monotonic.
Figure 14 presents the characteristic for modulation with sinusoidal and rectangular signals. This characteristic was reconstructed using the algorithm presented in Figure 2, where the modulation depth is replaced by frequency deviation .
5. Discussion of Results
Based on the presented figures, the following conclusions could be derived.
(1) Fading of Initial Transient Component
the fading out time of transient component of IFL signal for zero initial conditions (Figure 3) is about 120 s (even though the simulated signal chain includes a block of raising/falling time equal to 60 s, and the fading out time of the transient component of the block equals to 180 s).
(2) For a Case of No Modulation
for a case of no modulation, the value of indicator is greater than zero (see Table 1) and depends on the order of low-pass F1B filter. For order 6 of this filter, which is recommended in , the indicator value is about 0.01. This means that the measurement result of a real flickermeter cannot be lower than 0.01.
(3) For AM Modulation
(i)if an input signal is a result of AM modulation (2), then the processing characteristic (Figure 4) covers frequencies up to 155 Hz, and three local maxima occur for frequencies and Hz. The global maximum does not depend on the shape of the modulating signal and occurs at Hz. These signify that there are three maxima with respect to sensing the obnoxiousness of the flicker for incandescent lamp,(ii)characteristics (Figure 4), (Figure 5), and (Figure 6) show that for the threshold value of indicator (), the smallest modulation depth occurs for rectangular signal and the greatest modulation depth occurs for triangular signal. Ipso facto, the most obnoxious modulation is the modulation with rectangular signal, followed by the one with sinusoidal signal, and the least obnoxious is the modulation with triangular signal,(iii)with regard to characteristics (Figure 4), (Figure 5), and (Figure 6), one can observe a non-monotonicity in a frequency range of 28–37 Hz, and hence, this fragment of the characteristics is very useful during the tests of the performance of real flickermeters,(iv)by taking characteristic (Figure 5) as a reference criterion when evaluating flickermeter linearity, we can state that the flickermeter signal chain is, in general, nonlinear, but for inputs that correspond to (which means sensing of flicker), this system is nearly linear,(v)in a case of AM modulation, the influence of the input signal bandwidth limitation (Figure 8) is visible for Hz.
(4) For Input Signal with Interharmonic Component
comparative combination of and characteristics (Figure 9) makes a good basis to conclude on the difference between obnoxiousness of the flicker caused by amplitude modulation and occurrence of single interharmonic in a voltage that supplies incandescent lamp,
(5) For FM Modulation
(i)on the basis of characteristics for FM modulation (Figure 10), it is difficult to estimate the range of frequency in which a changeability of indicator value occurs; for modulation with sinusoidal signal, the changeability of indicator fades for Hz,(ii) characteristic for FM modulation with rectangular signal (Figure 10) is strongly non-monotonic,(iii)on the basis of plot (Figure 13), one can state that the indicator value could be greater than that for the frequency deviation greater than 0.25 Hz. Thus, FM modulation of the input voltage with Hz should not lead to obnoxious flicker (Figure 14),(iv)in the case of FM modulation, the measurement result of indicator strongly depends on the flickermeter bandwidth; limitation of input signal bandwidth (i.e., decreasing ) surprisingly leads to increased value of indicator (Figures 11 and 12),(v)taking characteristic (Figure 13) as a reference criterion when evaluating flickermeter linearity, we can state that the flickermeter signal chain in the case of FM modulation is, in general, nonlinear.
The results give a comprehensive overview of the signal chain and supplement the standard specification for the case of AM modulation of the input signal. The results also complement the specification for the case of FM modulation of the input signal and for the input signal with single interharmonic component. The presented results thoroughly describe the performance of IEC flickermeter in a full frequency range that influences the result of indicator measurement (as opposed to other results given in the literature that describe the flickermeter only in limited frequency range). The results of the simulations make it easier to understand the operation of the IEC flickermeter. They describe the influence of the input voltage parameters on indicator measurement result. Furthermore, the reaction of the IEC flickermeter to different types of input signals is also demonstrated. The analysis of the presented characteristics helps to determine the requirements with regard to flickermeter signal chain and suggests the potential source of measurement error. Any peculiar fragments of the characteristics define the optimal condition for checking the accuracy of the performance of the IEC flickermeter and, at the same time, help to shorten the time of flickermeter testing. It can be presented in the future.
- Flickermeter, International Electrotechnical Vocabulary, IEC, number 604-01-28.
- Flicker, International Electrotechnical Vocabulary, IEC, number 161-08-13.
- IEC 61000-4-15, Flickermeter—Functional and Design Specifications, 2010.
- H. De Lange Dzn, “Eye’s response at flicker fusion to square-wave modulation of a test field surrounded by a large steady field of equal mean luminance,” Journal of the Optical Society of America, vol. 51, no. 4, pp. 415–421, 1961.
- X. Yang and M. Kratz, “Power system flicker analysis and numeric flicker meter emulation,” in Proceedings of the IEEE Lausanne POWERTECH, pp. 1534–1539, July 2007.
- R. Cai, J. F. G. Cobben, J. M. A. Myrzik, J. H. Blom, and W. L. Kling, “Flicker responses of different lamp types,” IET Generation, Transmission and Distribution, vol. 3, no. 9, pp. 816–824, 2009.
- L. W. White and S. Bhattacharya, “A discrete matlab-simulink flickermeter model for power quality studies,” IEEE Transactions on Instrumentation and Measurement, vol. 59, no. 3, pp. 527–533, 2010.
- P. Clarkson and P. S. Wright, “Sensitivity analysis of flickermeter implementations to waveforms for testing to the requirements of IEC, 61000-4-15,” IET Science, Measurement and Technology, vol. 4, no. 3, pp. 125–135, 2010.
- J. Slezingr and J. Drapela, “Verification of Flickermeters under new edition of IEC, 61000-4-15,” in Proceedings of the IEEE Trondheim PowerTech, p. 6, June 2011.
- I. Sadinezhad and V. G. Agelidis, “Frequency adaptive least-squares-kalman technique for real-time voltage envelope and flicker estimation,” IEEE Transactions on Industrial Electronics, vol. 59, no. 8, pp. 3330–3341, 2012.
- D. Gallo, C. Landi, R. Langella, and A. Testa, “Implementation of a test system for advanced calibration and performance analysis of flickermeters,” IEEE Transactions on Instrumentation and Measurement, vol. 53, no. 4, pp. 1078–1085, 2004.
- E. W. Gunther, “A proposed flicker meter test protocol,” in Proceedings of the Quality and Security of Electric Power Delivery Systems Symposium, pp. 235–240, October 2003.
- L. Peretto, E. Pivello, R. Tinarelli, and A. E. Emanuel, “Theoretical analysis of the physiologic mechanism of luminous variation in eye-brain system,” in Proceedings of the IEEE Instrumentation and Measurement Technology Conference (IMTC '05), pp. 128–133, 2005.
- M. De Koster, E. De Jaeger, and W. Vancoetsem, “Light flicker caused by interharmonics,” in IEEE Power Engineering Society Transmission and Distribution Committee General Systems Subcommittee, Harmonics Working Group, 2001.
- S. M. Halpin and V. Singhvi, “Limits for interharmonics in the 1-100-Hz range based on lamp flicker considerations,” IEEE Transactions on Power Delivery, vol. 22, no. 1, pp. 270–276, 2007.
- T. Kim, E. J. Powers, W. M. Grady, and A. Arapostathis, “Detection of flicker caused by interharmonics,” IEEE Transactions on Power Delivery, vol. 58, no. 1, pp. 152–160, 2009.
- G. Wiczyński, “Analysis of flickermeter’s signal chain for input signal with two sub/interharmonics,” Electrical Review, vol. 86, pp. 328–335, 2010.
- D. Gallo, R. Langella, and A. Testa, “Toward a new flickermeter based on voltage spectral analysis,” in Proceedings of the International Symposium on Industrial Electronics (ISIE '02), vol. 2, pp. 573–578, 2002.
- D. Gallo, R. Langela, and A. Testa, “Light flicker prediction based on voltage spectral analysis,” in Proceedings of the IEEE Porto Power Tech Conference (PPT '01), vol. 1, p. 6, 2001.
- W. Mombauer, “Flicker caused by interharmonics,” EtzArchiv, vol. 12, no. 12, pp. 391–396, 1990.
- A. Testa and R. Langella, “Power system subharmonics,” in Proceedings of the IEEE Power Engineering Society General Meeting, vol. 3, pp. 2237–2242, June 2005.
- M. Rogóz, A. Bień, and Z. Hanzelka, “The influence of a phase change in the measured voltage on flickermeter response,” in Proceedings of the 11th International Conference on Harmonics and Quality of Power (ICHQP '04), pp. 333–337, September 2004.
- J. J. Gutierrez, L. A. Leturiondo, J. Ruiz, A. Lazkano, P. Saiz, and I. Azkarate, “Effect of the sampling rate on the assessment of flicker severity due to phase jumps,” IEEE Transactions on Power Delivery, vol. 26, no. 4, pp. 2215–2222, 2011.
- J. Ruiz, A. Lazkano, J. J. Gutierrez, et al., “Influence of the carrier phase on flicker measurement for rectangular voltage fluctuations,” IEEE Transactions on Instrumentation and Measurement, vol. 61, no. 3, pp. 629–635, 2012.
- W. Mombauer, “Flicker caused by phase jumps,” European Transactions on Electrical Power, vol. 16, no. 6, pp. 545–567, 2006.
- G. Wiczyński, “Simple model of flickermeter signal chain for deformed modulating signals,” IEEE Transactions on Power Delivery, vol. 23, no. 4, pp. 1743–1748, 2008.
- S. Haykin, Communication Systems, John Wiley & Sons, 1994.
- G. Wiczyński, “A model of the flickermeter for frequency modulation of the input voltage,” IEEE Transactions on Instrumentation and Measurement, vol. 58, no. 7, pp. 2139–2144, 2009.
- A. Bertola, G. C. Lazaroiu, M. Roscia, and D. Zaninelli, “A Matlab-simulink flickermeter model for power quality studies,” in Proceedings of the 11th International Conference on Harmonics and Quality of Power (ICHQP '04), pp. 734–738, September 2004.
- D. Hanselman and B. Littlefield, Mastering Matlab 6: A Comprehensive Tutorial and Reference, Prentice Hall, New Jersey, NJ, USA, 2001.