- Email: [email protected]
Estimator-analyzer of power quality: Part I – Methods and algorithms
Measurement 44 (2011) 238–247
Contents lists available at ScienceDirect
Measurement journal homepage: www.elsevier.com/locate/measurement
Estimator-analyzer of power quality: Part I – Methods and algorithms Tomasz Tarasiuk ⇑ Gdynia Maritime University, Department of Ship Electrical Power Engineering, Morska 83, 81-225 Gdynia, Poland
a r t i c l e
i n f o
Article history: Available online 3 November 2010 Keywords: Power quality Signal processing Algorithms Digital signal processors Spectral analysis
a b s t r a c t The paper presents original methods and algorithms of measurement of power quality parameters which were implemented in a new measurement device called an estimatoranalyzer of power quality. The main idea was to overcome the shortcomings of current power quality analyzers, namely sensitivity to fast fluctuations of input signal parameters. This was attained with regard to the algorithm’s overall computational complexity. Complementary application of discrete wavelet transform (DWT) with Fourier transforms and chirp z-transform (CZT) was proposed for the aim. Next, a new concept of estimator mode was devised. It consists in multistage signal processing and applied methods depending on input signal character – its frequency and existing waveform distortions – as well as the chosen set of reference parameters. It further improves device performance in terms of time of computation. Description of the estimator-analyzer hardware and results of experimental research will be presented in the second part of the paper. Ó 2010 Published by Elsevier Ltd.
1. Introduction Nowadays, power quality and devices for its monitoring are hotly discussed topic in the field of power systems engineering. The term ‘‘power quality” (PQ) has various meaning but usually it is considered as a combination of voltage and current quality [1,2] or sometimes only the voltage quality [3]. In this paper the term power quality will be used in the latter meaning. It worth adding that the similar digital signal processing techniques can be used for determining both: the voltage and the current quality. The indispensable part of power quality assessment is the measurement of different voltage parameters. The measurement of PQ parameters is usually carried out by dedicated devices called analyzers of power quality. These devices are capable of measurement of many various parameters but arguably the most investigated part of their algorithm is devoted to waveform assessment of the input signal. It consumes a great deal of resources of the measurement device (computational power and RAM memory). Numerous studies have been carried out and ⇑ Tel.: +48 58 690 14 14; fax: +48 58 620 67 01. E-mail address: [email protected] 0263-2241/$ - see front matter Ó 2010 Published by Elsevier Ltd. doi:10.1016/j.measurement.2010.09.049
various methods of estimation of waveform parameters have been proposed. The studies concern methods of evaluation of low frequency disturbances, transient detection and evaluation or data compression. The methods for assessment of low frequency phenomena up to 9 kHz are based on various principles, e.g. methods based on wavelet packet decomposition (WPD) [4], neural networks [5], nonlinear adaptive technique [6], digital filters [7,8] or mathematical morphology [9]. But by far the most often investigated method is still discrete Fourier transform (DFT), especially its fast algorithm (FFT). It is the method recommended in current IEC Standard 61000-4-7 [10]. The FFT algorithm has been implemented in PQ analyzers for many years. In the 1980s, the method was implemented in a device equipped with a Motorola 16-bit 68,000 microprocessor clocked at 4.9152 MHz [11]. The applied constant sampling frequency was 7680 Hz. Unfortunately, FFT has some commonly known downsides. In particular, it is hard achieving both: desired window width and appropriate number of input signal samples. So, synchronization of sampling frequency with actual input signal frequency has been investigated [12,13] and implemented. Additionally different interpolation techniques have been developed [14]. The two techniques are
T. Tarasiuk / Measurement 44 (2011) 238–247
commonly applied in commercial power quality analyzers. Further, signal windowing is still considered [15] and commonly added as an option in digital oscilloscopes. The author’s research has demonstrated that all the methods based on FFT lack versatility and they can lead to dubious results for some real signals registered in electric power systems [17,18]. Simply put, observed maximum errors were above the permissible limit set in the previously mentioned IEC Standard [10]. However, microprocessor technique and especially digital signal processor (DSP) technique has developed rapidly over the last few years. Present DSPs can be clocked at hundreds of MHz (e.g. 32 bit ADSP 21364 maximum clock frequency 333 MHz [16]) and their Harvard architecture and one clock cycle multiplication enable processing of large volumes of data in real time. So, should we still use the above mentioned and recommended in relative standard methods of signal processing based on FFT? The problem of proper waveform assessment increases in difficulty if high frequency phenomena such as transients have to be analyzed. This entails an increase in sampling frequencies and the introduction of additional tools of signal processing, e.g. application of DWT [19–21], mathematical morphology [22] or subtracting the preevent steady-state waveform with sliding window and varying thresholds [23]. In the paper and its second part the original methods and prototype of the new instrument for PQ monitoring is described. It was designed for monitoring of PQ in relatively broad frequency bandwidth and the FFT as basic tool of waveform assessment was discarded. The basic tools of harmonic and interharmonic subgroups measurement are DFT and chirp z-transform (CZT – devised as a fast method of estimation of Fourier coefficients for an arbitrary number of samples [24]), applied for an exact integer number of input signal cycles, independently of its momentary frequency. Therefore, no sampling frequency synchronization or signal windowing is needed. The instrument was called estimator-analyzer of PQ since it has two main modes of operation: analyzer and estimator. Both modes are described in this paper but especially the idea of estimator mode has to be underlined. It is based on the concept of multistage signal processing with varying algorithms, depending on the input signal characteristic and user defined set of reference parameters. In this mode the conformity of PQ parameters with their permissible limits is monitored, with as limited time of computation as possible. The features and merits of methods applied in the original instrument were discussed in the author’s previous papers [17,18,25–30] separately. Therefore, in this paper only the main conclusions from these papers are summarized and related previous papers are recalled. Eventually, the methods were combined as one complex algorithm, appropriate software was developed and the estimator-analyzer was constructed. So, in the paper and its second part the details of the previously unpublished algorithms, hardware and research results of the new instrumentation are presented.
239
2. Standard framework for measuring spectral components in low frequency range IEC Standard 61000-4-7:2002 recommends methods of measurement of harmonics, interharmonics and components above the harmonic frequency range up to 9 kHz [10]. Although other principles of analysis are not excluded, DFT is specified as the reference instrument. According to the standard, new instruments are likely to use DFT, normally using the FFT algorithm. Further, synchronization of measurement time window with actual input signal periods is recommended. The window duration should be 10 periods for 50 Hz systems or 12 periods for 60 Hz systems with rectangular weighting. A Hanning window can be applied only in the case of loss of synchronization. Maximum permissible error of sampling frequency synchronization with actual duration of recommended window width should not exceed ±0.03% [10]. Further, the concept of harmonic subgroup was introduced in the discussed standard [10] in order to account for signal magnitude fluctuations and improve waveform assessment accuracy. The idea consists in calculating the square root of the sum of squares of harmonic component amplitude and amplitudes of two spectral components (interharmonics) directly adjacent to it. The result is considered as the rms value of the harmonic subgroup. The interharmonic centred subgroup is calculated in a similar way on the basis of remaining spectral bins between two harmonic subgroups [10]. However, it was concluded in Refs. [17,18] that the methods recommended in IEC Standard 61000-4-7 [10] could lead to significantly flawed results in some rare but real cases. In the discussed standard method of assessment of the components above the harmonic frequency range up to 9 kHz was proposed. They are to be calculated on the basis of a 100 ms window corresponding to approximately five or six fundamental periods. The FFT with a rectangular window can be used for the measurement, since sampling frequency need not be synchronized with the actual one. But attenuation of the fundamental component is required and it should exceed 55 dB [10]. Finally, the components calculated by FFT are to be grouped in bands of 200 Hz, according to formula [10]:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u bþ100 u X Gb ¼ t C 2f
ð1Þ
f ¼b90ðHzÞ
where Gb – group of bth frequency band, b – centre frequency of considered frequency band, Cf – rms value of spectral line. The above concepts of harmonic subgroups, interharmonic centred subgroups and groups of components in the frequency range from the 50th harmonic up to 9 kHz are utilized in the presented estimator-analyzer of PQ. 3. Estimator-analyzer of power quality 3.1. Hybrid wavelet-Fourier method of spectrum analysis The hybrid wavelet-Fourier method is the basis of all methods and algorithms behind the software of the
240
T. Tarasiuk / Measurement 44 (2011) 238–247
presented estimator-analyzer of PQ. It consists in complementary application of discrete wavelet transform (DWT) and various methods of computation of Fourier series coefficients [17,25,30]. Basically, DWT is applied for transient and notching detection and evaluation. Transients usually contain high frequency components from a few 100 Hz up to a few MHz but their duration is up to 50 ms (in real cases even less) [31]. Unfortunately, the high sampling frequency required for transients analysis would be an enormous burden for the measurement device in the case of low frequency phenomena measurement. Analysis of components in the frequency range up to 9 kHz requires lower sampling frequencies but longer time of observation (window width equals approximately 200 ms for components in the harmonic frequency range and 100 ms above this range up to 9 kHz) [10]. Fortunately, wavelet coefficients after signal decomposition contain all required information about low-frequency components [4,32,33]. Moreover, the number of wavelet coefficients is decimated after each decomposition level [33] and only coefficients representing lower frequencies should be preserved for the whole measurement window for measurement of low frequency parameters. Therefore, if one calculates coefficients of Fourier series on the basis of these wavelet coefficients (instead of original samples of input signal) a huge decrease in the required mathematical operation would be achieved (in practical terms by a factor of 24) [17,25]. Additionally, RAM memory of the used DSP would be saved as well [27]. One more feature of DWT was utilized. Wavelet approximation coefficients which represent low-frequency components with the fundamental one can be treated as a smoothed version of the original input signal. Simply put, the higher frequency components are filtered out. Then, measurement window width can be calculated by analysis of these coefficients’ zero-crossings, improved by simple interpolation. That way, the influence of multiple zero-crossing or higher order harmonics would be minimized [28]. This method of window width estimation was implemented in the described measurement device. 3.2. Estimator-analyzer of power quality – functions overview The basic measurement functions of the prototype estimator-analyzer are as follows [26]: – – – – – – – – – –
–
voltage rms value, frequency, unbalance coefficient, newly introduced distortion band factors (DBF) [27], voltage dip (value and duration), harmonic subgroups up to the 50th harmonic, interharmonic centred subgroups in the range up to the 50th harmonic, group of sub-harmonics, DC component, groups of high frequency components with bandwidths approximately 200 Hz in the range from the 50th harmonic up to 9 kHz, transient and notching detection and assessment,
– power quality temperature factor for assessment of additional temperature rise of windings of induction motors [34]. In all cases, a waveform analysis is carried out for exactly 10 periods of input signal (the prototype was constructed for 50 Hz systems) and only a rectangular window is adopted, independently of input signal character. Voltage dips are estimated on the basis of voltage rms values calculated for one cycle of input signal refreshed each half-cycle Urms(1/2) by the method described in IEC Standard 61000-4-30 [35]. Transients and notching are detected and evaluated on the basis of wavelet coefficients after filtering out low-frequency components by DWT. 3.3. Analyzer mode The discussed measurement device operates in two basic modes: analyzer and estimator. But a significant part of the developed software is common to both modes. The device operates with constant sampling frequency (148917.8 Hz) and first DWT is applied for signal samples in both modes. Daubechies filters with 10 filter coefficients and eight levels of wavelet decomposition are applied. The differences appear after DWT calculation. The estimator mode consists in multistage signal processing with an algorithm that depends on actual characteristics of the input signal. It will be described in the following subsection. In the analyzer mode the device operates on the basis of a constant algorithm. This means that all required PQ parameters are always determined independently of the actual characteristics of the input signal (frequency and waveform distortions). After DWT the following main tasks are carried out: 1. The measurement window width and frequency are calculated on the basis of zero-crossing detection of wavelet approximation coefficients of the seventh decomposition level. These coefficients represents frequency components up to 581.7 Hz (cut-off frequency of respective filters of seventh decomposition level). 2. Signal rms value is calculated on the basis of wavelet coefficients. Details are set out in the following subsection. 3. Partial synthesis of wavelet detail coefficients of third, fourth, fifth, sixth and seventh decomposition levels is carried out. FFT is applied to results of the synthesis in order to estimate frequency components above harmonic frequency range up to 9 kHz. The method enables attenuation of the fundamental component as required by the IEC Standard [10]. Details of the method were provided in Ref. [29]. 4. CZT is applied for wavelet approximation coefficients of the fourth decomposition level and subsequently harmonic subgroups and interharmonic centred subgroups are calculated. The wavelet approximation coefficients of the fourth decomposition level represent frequency components below 4653.7 Hz (cut-off
241
T. Tarasiuk / Measurement 44 (2011) 238–247
frequency of respective filters of fourth decomposition level), so it contains components of harmonic frequency range. The important part of the algorithm is CZT. It is used for calculation of components in the range up to the 50th harmonic, since it enables relatively fast calculation of DFT for an arbitrary number of samples [24]. The block diagram of CZT for the considered aim is depicted in Fig. 1 [24,29,36,37]. It consists of two FFTs and one inverse fast Fourier transform (IFFT). The W constant can be defined as follows [30]:
j 2p W ¼ exp n
ð2Þ
n value represents a theoretical (non-integer) number of wavelet coefficients of the fourth decomposition level and it is calculated by means of the formula:
n¼
Tw fs
ð3Þ
24
where Tw – duration of considered measurement window calculated in the above-described step 1, fs – sampling frequency. Since only components’ magnitudes are estimated, the last multiplication in the CZT algorithm was omitted. The input data sets y1(k) and y2(k) can be calculated by formulae (4) and (5) [29,30].
k2 2 y1ðkÞ ¼ C 4;k W 0 k2 2 W y2ðkÞ ¼ 0 ðkMÞ2 W 2
k ¼ 0; . . . ; N 1 k ¼ N; . . . ; M 1 k ¼ 0; . . . ; N 1 k ¼ N; . . . ; M N 1
ð4Þ
ð5Þ
k ¼ M N; . . . ; M 1
where c4,k – wavelet approximation coefficients of fourth decomposition level, N – number of coefficients calculated for 10 cycles (in this case n from formula (3) means theoretical (non-integer) number of wavelet approximation coefficients of considered decomposition level), M/2 – the nearest power of 2 greater than N. Further, part of the CZT algorithm is used for measurement of frequency components above the harmonic frequency range up to 9 kHz. Since these components may be calculated on the basis of 5 or 6 cycles of the input signal [10] (a window two times shorter than the measurement window for harmonics and interharmonics calculation) the software for parallel FFTs (see Fig. 1) is used for measurement of these components. Data calculated by means of wavelet synthesis during the above-mentioned step 3 for the first half of the basic 10-cycle measurement window replace array y1(k) and the respective data calculated for the second half replace array y2(k). This enables optimization of the two FFTs’ algorithm and improves performance of the algorithm for measurement of the higher frequency components in terms of computational complexity [30]. Finally, a block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels and
Fig. 1. Block diagram of DFT calculation (X(k)-bins of Fourier spectrum) by means of CZT.
main subroutines in analyzer mode is presented in Fig. 2 (excluding transients and dips detection and evaluation). The main difference between the considered device operating in analyzer mode and commercial analyzers consists in the DFT processor, as defined in IEC Standard 61000-4-7 [10]. The measurements of components in the harmonic frequency range are carried out on the basis of the combined DWT–CZT algorithm and FFT is discarded altogether, so synchronization of sampling frequency is not needed. This means that coefficients of Fourier series are always calculated on the basis of an integer number of input signal cycles. So, the basic idea of waveform assessment set out in IEC Standard 61000-4-7 [10] is still preserved. 3.4. Estimator mode The estimator mode is a completely new concept. It is based on the assumption that in many cases only conformity of PQ parameters with their permissible limits specified in related standards or defined by user are necessary. So, it would be useless to determine all PQ parameters and find out that most of them or all are well below their maximum permissible limits. This realization is the basis of the estimator mode algorithm of the discussed measurement instrument. The purpose of measurements in this mode is detection of breaching of the user defined limits, with as reduced number of processor’s instructions as possible. This aim is achieved by multistage signal processing. The algorithm applied at each stage depends on the results of the previous one. After DWT four stages of signal processing are implemented. The algorithm of multistage waveform assessment is graphically presented in Fig. 3 and a block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels and main subroutines in the estimator mode is presented in Fig. 4. 3.4.1. First stage In this stage basic PQ parameters are calculated, namely input signal rms value, its frequency, unbalance coefficient as well as THD coefficient defined as the ratio of rms value remaining after subtracting the fundamental component to the fundamental component [38].
THD ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U 2rms U 21 U1
100
ð6Þ
where Urms – input voltage rms value, U1 – fundamental component rms value.
242
T. Tarasiuk / Measurement 44 (2011) 238–247
Fig. 2. Block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels (from third to seventh levels of signal wavelet decomposition) and main subroutines in analyzer mode for measurements of low frequency PQ parameters. Dashed line marks flow of information among different subroutines and processor’s Data Memory bank DM required to carry out calculation; h – low-pass filter bank, g – high-pass filter bank, c – wavelet approximation coefficients, d – wavelet detail coefficients.
The value of U1 is calculated by means of ‘‘ordinary” DFT applied for wavelet approximation coefficients of the seventh decomposition level as the input signal. The average number of these coefficients calculated for the duration of the basic measurement window is 27 lower than the number of input signal samples taken by the A/D converter during the same time interval. Obviously, it has a positive impact on the algorithm’s computational complexity. The computation of DFT, CZT or FFT requires appropriate sine and cosine values. It should be mentioned that generally these values are taken from the sine table (4096 values) located in the processor’s data memory bank DM. But THD assessment is an exception. In this case, sine and cosine values are calculated by a min–
max polynomial approximation algorithm [39]. The former method is much faster, but it can lead to unacceptable error in case of low value THDs (below 2% in practical terms due to used formula (6)) and some singular input signal frequencies. This phenomenon was observed during the first tests of the presented device in a ship network. Because increase in volume of the sine table was not possible, other solutions were investigated. A method of calculation of sine and cosine values based on the sine table and interpolation was tested, but the results were mixed at best. Finally, the above-mentioned method of min–max polynomial approximation algorithm led to satisfactory results and it was used for calculation of only U1 in estimator mode.
T. Tarasiuk / Measurement 44 (2011) 238–247
243
Fig. 3. Algorithm of multistage waveform distortion assessment in low frequency range up to 9 kHz in estimator mode, PV – permissible value for chosen set of parameters – 0.5% in the case of EN Standard 50160 [40] or 3% in the case of IEC Standard 60092-101 [38], ncDFT – relative number of mathematical operations for DFT and actual input signal waveform distortion, ncCZT – relative number of mathematical operations for CZT and actual input signal waveform distortion.
If the calculated THD value is lower than the maximum permissible value (PV) of any required higher frequency
component the whole analysis for the considered measurement window is finished and further analysis stages are
244
T. Tarasiuk / Measurement 44 (2011) 238–247
Fig. 4. Block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels (from third to seventh) and main subroutines in estimator mode for measurements of low frequency PQ parameters (third to seventh levels of signal wavelet decomposition); dashed line marks flow of information among different subroutines and processor’s Data Memory bank DM required to carry out calculation; h – low-pass filter bank, g – high-pass filter bank, c – wavelet approximation coefficients, d – wavelet detail coefficients.
omitted. PV are chosen by the user and it depends on the chosen set of permissible values of PQ parameters. There are two predefined sets of these values included in the RAM memory of the used DSP. The first set is based on EN Standard 50160 [40]. The second set is based on IEC Standard 60092-101 [38]. In the case of the former standard [40] PV was assumed to be equal to 0.5%. In the case of the latter standard [38] PV was assumed to be equal to 3.0%.
3.4.2. Second stage This stage consists in rough waveform assessment by means of distortion band factors (DBF). The concept of these new factors was described in previous papers by the author, most importantly in Ref. [27]. In shorthand, DBFs are calculated by means of the formula:
DBF miþ2 ¼
U rmsðfmi fmði1Þ Þ U1
100½%
ð7Þ
245
T. Tarasiuk / Measurement 44 (2011) 238–247
where U rmsðfmi fmði1Þ Þ – rms value of voltage components of (fmi – fm(i1)) frequency band, U1 – rms value of the fundamental component. The rms value of signal components for each frequency band can be calculated by means of wavelet detail coefficients for the respective decomposition level as follows [27]:
U rmsðfmi fmði1Þ Þ
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u i 1 u 1 NX ¼t i ðdmi;k Þ2 2 ni k¼0
ð8Þ
where dmi,k – wavelet detail coefficients of respective decomposition level, Ni – number of wavelet coefficients of considered decomposition level calculated for a time interval of basic measurement window, n – theoretical (non-integer) number of wavelet coefficients in an analysed window calculated by means of appropriately modified formula (3), i – calculation order of the respective decomposition level (i 2 1. . . m), m – number of decomposition levels (m = 8 for the described device). The last DBF for the lowest frequency band is calculated differently, by the formula [27]:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u m X u DBF 1 ¼ tTHD2 DBF 2miþ2
ð9Þ
i¼1
DBF1 comprises components up to a cut-off frequency of low-pass filter of the last decomposition level without the fundamental component. DBFs can be calculated without an additional requirement for computational power of the measurement device. The considered signal rms value was calculated concurrently with these DBFs by the formula [27,41]:
U rms
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # u N m1 m1 i 1 u X NX X1 1 2 t 2 ¼ cm7;k þ ðdmi;k Þ 2m1 nm1 k¼0 i¼1 k¼0 ð10Þ
where cm7,k – wavelet approximation coefficients of seventh decomposition level, m = 8. One should take into account that the applied number of wavelet coefficients N in formula (8) is decimated for each consecutive layer of wavelet decomposition. This means that calculation of DBFs together with the considered signal rms value by means of wavelet coefficients requires a roughly similar number of mathematical operations as the calculation of only the signal rms value on the basis of the original samples of the signal. So, implementation of the method does not increase the computational complexity of the algorithm for PQ estimation on the one hand. But it increases the scope of information about the signal behaviour on the other hand [27]. (In fact, DBFs are calculated independently of THD value. If THD < PV further analysis based on DBFs is skipped.) DBFs can be treated as a rough measure of waveform distortions in the whole frequency band that is determined by the cut-off frequency of the anti-aliasing filter. More information about possibilities of DBFs utilization can be found in Ref. [27]. Since each DBF represents information about waveform distortion in the respective frequency band, it can be used
in a similar way as the THD factor. Namely, the content of the respective frequency band should be analysed pffiffifficarefully only if the product of corresponding DBF and 2 exceeds the permissible limit of any higher frequency component from the analyzed frequencyp band. The necessity of multiffiffiffi plication value of DBF with 2 results from the shape of wavelet filters response [27]. The permissible limit depends on the chosen set of reference values and pffiffiffi analysed frequency band. If all products of all DBFs and 2 are lower than the respective permissible limits, the whole analysis is finished and the next steps are omitted. Otherwise a further analysis is required. In the case of chosen standards all permissible limits equal PV defined for THD assessment in the first stage. Firstly, if values of DBF6 or DBF5 corresponding with wavelet detail coefficients of the fourth and fifth decomposition levels pffiffiffi (frequencies 2326.8–9307.4 Hz) are greater than PV/ 2, a subroutine for calculation of components above the harmonic frequency range up to 9 kHz is called. (The same subroutine is always called in analyzer mode.) Secondly, an affected frequency band in the harmonic frequency range is determined by analysis of respective DBFs associated with the harmonic frequency range. These DBFs, corresponding frequency bands and number of decomposition levels are presented in Table 1 (assuming sampling frequency 148917.8 Hz). By DBFs analysis a set of wavelet coefficients for detailed analysis of harmonics or interharmonic subgroups is defined. It should be mentioned that in many cases only lower order harmonics are presented in real power system signals. It diminishes algorithm computational complexity due to decimation of the wavelet coefficient after each decomposition level. For example, DBF2 is calculated on the basis of wavelet detail coefficients of the eighth pffiffiffi decomposition level. If DBF2 has a value above PV/ 2, a computation of Fourier series coefficients would be carried out on the basis of wavelet approximation coefficients of the sixth decomposition level. This way a possible magnitude damping of Fourier coefficients (due to the shape of wavelet filters’ frequency response) would be less than 0.01 and aliasing of components above 581.7 Hz would be minimized as well. Fortunately, the number of wavelet approximation coefficients of the sixth decomposition level is decimated by a factor of 26 in comparison to the number of input signal samples. As a result a decrease in computational complexity would be achieved. 3.4.3. Third stage In this stage FFT is performed. Input data set for the procedure consists of wavelet coefficients defined in the Table 1 DBFs of harmonic frequency range and related frequency bands. DBF
Frequency band (Hz)
Number of wavelet decomposition level
DBF5 DBF4 DBF3 DBF2 DBF1
2326.8–4653.7 1163.4–2326.8 581.7–1163.4 290.9–581.7 Below 290.9 without fundamental component
5 6 7 8 –
246
T. Tarasiuk / Measurement 44 (2011) 238–247
previous stage by the method described above. The sampling frequency (148917.8 Hz) was chosen in such a way that for all possible frequencies of input signal in steady-state conditions (50 Hz ±5%) the number of input data for CZT in analyzer mode would equal 4096. At the same time the number of wavelet coefficients calculated for 11 cycles of input signal would be the power of two for each decomposition level, if the actual input signal frequency was equal to the rated one (in fact 49.99 Hz; exact synchronization for 50 Hz was not possible) [29]. Finally, such a synchronization condition is checked (measured frequency within 49.99 Hz ±0.03% limit). If it is fulfilled, an FFT subroutine is called and a rectangular window is used. Then, harmonic subgroups and interharmonic centred subgroups are calculated and compared with their permissible values (specified in the relevant standard). Then the analysis is finished. If synchronization was not obtained, Hanning weighting is applied prior to FFT analysis and the fourth stage of analysis is carried out. 3.4.4. Fourth stage Because application of a Hanning window can lead to flawed results in some specific but real signals [17,18] further processing is necessary. In this stage a final method of signal processing is chosen. In Ref. [42] a windowed (Blackman) FFT was proposed for coarse spectrum analysis in order to localize frequency bands with significant spectral energy. Next, windowed CZT is applied for a medium search in the specified frequency band. CZT is applied for evaluation of the exact frequencies of components for detailed analysis. Subsequently, evaluation of exact magnitudes and phases is carried out by DFT [42]. This seems too complicated and hardly practical in the case of significantly distorted signals. So, the author proposes a modified procedure. First the relative computational complexity of DFT is estimated and compared with the relative computational complexity of CZT. The calculation of relative computational complexity takes into account the exact wavelet decomposition levels which contain components of interest as well as results of preliminary assessment by FFT with Hanning window. It can be carried out relatively fast but requires prior assumption of the appropriate threshold. Relative computational complexity of CZT depends on utilized wavelet coefficients (decomposition level), whereas relative computational complexity of DFT depends on input signal character – actual harmonic or interharmonic content or more precisely their frequencies in relation to wavelet decomposition levels (DFT can be applied for wavelet coefficients of various decomposition levels depending on order of actually calculated frequency bin [25]). Finally, a precise calculation of specified harmonic subgroups and interharmonics centre subgroups is carried out on the basis of appropriate wavelet coefficients by means of DFT or CZT, whichever would be faster. It should be added that final results are always calculated on the basis of a rectangular window, independently of actual signal frequency and used method (FFT, DFT or CZT) and despite constant sampling frequency. Moreover,
in practical cases, it should be faster than in analyzer mode (detailed analysis of time of computation will be presented in Part II). Comparison of Figs. 2 and 4 clearly indicates that most parts of developed subroutines are common in both modes. So, introduction of the estimator mode leads to only a slight increase in the size of final software and required fast RAM memory of the DSP used.
4. Final remarks The solution described in this paper is based on various methods for PQ assessment. These methods were integrated as one complete tool for profound PQ analysis. The resulting algorithm conforms to basic concepts of PQ measurement laid in the respective IEC Standard [10], although it differs from the methods directly recommended in the standard and commonly used by developers of current commercial PQ analyzers. During development of the algorithm three aims were to be achieved: increasing the sampling frequency as high as possible, increasing the algorithm immunity for relatively fast fluctuations of input signal parameters as well as saving resources of measurement device, such as computational power and RAM. The aims were obtained by complementary application of various tools of signal processing: DWT, FFT, DFT and CZT. This lead to increased performance of signal processing algorithm, in comparison to application of these tools separately. Its is also worth stressing the concept of the adaptation of processing algorithm to the actual character of the input signal (its frequency and waveform) by multistage signal processing. It was called the estimator mode and leads to a further improvement in terms of computational complexity. The introduction of this mode increases functionality of measurement device. It enables: increasing of sampling frequency or adding additional functions of measurement instrument, e.g. extending number of monitored parameters. Further, the mode can be useful in handheld devices supplied by build-in batteries. Simply put, if external power is switched-off the DSP core instruction rate can be decreased (lower computational power), which leads to energy saving and prolongs battery life. The exact gains of the implementing of the author’s proposal depend on the used hardware platform. The described solution was devised chiefly for DSP based instruments. Such a instrument, its structure and detailed results of its research will be described in Part II. Nevertheless, the idea can be also implemented in PC based instruments with standard data acquisition boards or instruments with combination of DSP and FPGA for signal processing. So, its possible applications are broader than the exemplary device described in Part II. References [1] M. Bollen, What is power quality?, Electric Power Systems Research 66 (2003) 5–14 [2] D. Katsaprakakis, D. Christakis, A. Zervos, S. Voutsinas, A Powerquality measure, IEEE Transactions on Power Delivery 23 (2) (2008) 553–561.
T. Tarasiuk / Measurement 44 (2011) 238–247 [3] R. Dugan, M. McGranaghan, H. Beaty, Electrical Power Systems Quality, McGraw-Hill Companies, 2002. [4] L. Eren, M. Unal, M.J. Devaney, Harmonic analysis via Wavelet packet decomposition using special elliptic half-band filters, IEEE Transactions on Instrumentation & Measurement 56 (6) (2007) 2289–2293. [5] L.L. Lai, W.L. Chan, C.T. Tse, A.T.P. So, Real-time frequency and harmonic evaluation using artificial neural networks, IEEE Transactions on Power Delivery 14 (1) (1999) 52–59. [6] D.M. McNamara, A.K. Ziarani, T.H. Ortmeyer, A new technique of measurement of nonstationary harmonics, IEEE Transactions on Power Delivery 22 (1) (2007) 387–395. [7] J.J. Tomic´, M.D. Kušljevic´, V.V. Vujicˇic´, A new power system digital harmonic analyzer, IEEE Transactions on Power Delivery 22 (2) (2007) 772–780. [8] L. Ferrigno, C. Landi, M. Laracca, C. Luongo, A study on the feasibility and effectiveness of digital filters approach for power quality monitoring in compliance with IEC 61000-4-7, in: 15th IMEKO TC4 Conference, Iasi, Romania, 2007. [9] P.M. Ramos, T. Radil, F.M. Janeiro, A. Cruz Serra, DSP based power quality analyser using new signal processing algorithms for detection and classification of disturbances in a single-phase power system, Metrology and Measurement Systems 14 (4) (2007) 483–494. [10] IEC Standard 61000-4-7, General Guide on Harmonics and Interharmonics Measurements for Power Supply Systems and Equipment Connected Thereto, 2002. [11] A.N. Mortensen, G.L. Johnson, A power system digital harmonic analyzer, IEEE Transactions on Instrumentation & Measurement 37 (4) (1988) 537–540. [12] M. Aiello, A. Cataliotti, V. Cosentino, S. Nuccio, Synchronization techniques for power quality instruments, IEEE Transactions on Instrumentation & Measurement 56 (5) (2007) 1511–1519. [13] A. Cataliotti, V. Cosentino, S. Nuccio, A phase-locked loop synchronisation of power quality instruments in the presence of stationary and transient disturbances, IEEE Transactions on Instrumentation & Measurement 56 (6) (2007) 2232–2239. [14] M. Sedlácˇek, M. Titeˇra, Interpolations in frequency and time domains used in FFT spectrum analysis, Measurement 23 (3) (1998) 185–193. [15] A. Testa, D. Gallo, R. Langella, On the processing of harmonics and interharmonics: using Hanning Window in Standard Framework, IEEE Transactions on Power Delivery 19 (1) (2004) 28–34. [16] Analog Devices, SHARC Processors. ADSP 21362/ADSP 21363/ADSP 21364/ADSP 21365/ADSP 21366. Data Sheet, 2007. [17] T. Tarasiuk, Hybrid wavelet-Fourier method for harmonics and harmonic subgroups measurement – case study, IEEE Transactions on Power Delivery 22 (1) (2007) 4–17. [18] T. Tarasiuk, Comparative study of various methods of DFT calculation in the wake of IEC Standard 61000-4-7, IEEE Transactions on Instrumentation & Measurement 58 (10) (2009) 3666–3677. [19] O. Poisson, P. Rioual, M. Meunier, Detection and measurement of power quality disturbances using wavelet transform, IEEE Transactions on Power Delivery 15 (3) (2000) 1039–1044. [20] L. Angrisani, P. Daponte, M. D’Apuzzo, A. Testa, A measurement method based on the wavelet transform for power quality analysis, IEEE Transactions on Power Delivery 13 (4) (1998) 990–998. [21] S. Santoso, W.M. Grady, E.J. Powers, J. Lamoree, S.C. Bhatt, Characterization of distribution power quality events with Fourier and wavelet transforms, IEEE Transactions on Power Delivery 15 (1) (2000) 247–254.
247
[22] S. Ouyang, J. Wang, A new morphology method for enhancing power quality monitoring system, Electrical Power and Energy Systems 29 (2007) 121–128. [23] D. Castaldo, D. Gallo, C. Landi, A. Testa, A digital instrument for nonstationary disturbance analysis in power lines, IEEE Transactions on Instrumentation & Measurement 53 (5) (2004) 1353–1361. [24] L.L. Rabiner, R. Shafer, C. Rader, The chirp z-transform, IEEE Transactions on Audio & Electroacoustic 17 (2) (1969) 86–92. [25] T. Tarasiuk, Hybrid wavelet-Fourier spectrum analysis, IEEE Transactions on Power Delivery 19 (3) (2004) 957–964. [26] T. Tarasiuk, A few remarks about assessment methods of electric power quality on ships – present state and further development, Measurement 42 (8) (2009) 1153–1163. [27] T. Tarasiuk, The method based on original DBFs for fast estimation of waveform distortions in ship systems – case study, IEEE Transactions on Instrumentation & Measurement 57 (5) (2008) 1041–1050. [28] T. Tarasiuk, Wavelet coefficients for window width and subsequent harmonic estimation – case study, Measurement 41 (3) (2008) 284– 293. [29] T. Tarasiuk, Method, algorithm and device for estimation of components above the harmonic frequency range up to 9 kHz, Measurement 44 (1) (2011) 219–229. [30] T. Tarasiuk, Application of CZT transform for spectrum analysis in systems with varying frequency, in: IEEE Int. Workshop on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications, 2007, Germany, pp. 329–334. [31] IEEE Standard 1159-1995. IEEE Recommended Practice for Monitoring Electric Power Quality. [32] W. Mikhael, A. Ramaswamy, An efficient representation of nonstationary signals using mixed-transforms with applications to speech, IEEE Transactions on Circuits and Systems – II Analog and Digital Signal Processing 42 (6) (1995) 393–401. [33] S. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence 2 (7) (1989) 674–693. [34] P. Gnacin´ski, J. Mindykowski, T. Tarasiuk, A new concept of the power quality temperature factor and its experimental verification, IEEE Transactions on Instrumentation & Measurement 57 (8) (2008) 1651–1660. [35] IEC Standard 61000-4-30, Testing and Measurement Techniques – Power Quality Measurement Methods, 2003. [36] C. Zang, L. Dai, H. Zheng, J. He, Using frequency zoom technology to realize high precision and adaptive frequency measurement for power system, in: IEEE Int. Conference POWERCON, Singapore, 2004, pp. 155–159. [37] T. Zielin´ski, Digital Signal Processing, WKŁ, Poland, 2005 (in Polish). [38] IEC Standard 60092-101, Electrical Installations in Ships. Definitions and General Requirement, 2002. [39] Analog Devices, ADSP-21000 Family. Applications Handbook, vol. 1, 1996. [40] EN Standard 50160, Voltage Characteristics of Electricity Supplied by Public Distribution Systems, 1999. [41] Weon-Ki Yoon, M. Devaney, Power measurement using the wavelet transform, in: IEEE Instrumentation and Measurement Technology Conference, 1998, pp. 801–806. [42] K.J. Runtz, D. Hack, A multistage DFT – FFT – CZT approach for accurate efficient analysis of sparsely distributed spectra, in: IEEE Canadian Conference on Electrical & Computing Engineering, 2002, pp. 127–132.
Contents lists available at ScienceDirect
Measurement journal homepage: www.elsevier.com/locate/measurement
Estimator-analyzer of power quality: Part I – Methods and algorithms Tomasz Tarasiuk ⇑ Gdynia Maritime University, Department of Ship Electrical Power Engineering, Morska 83, 81-225 Gdynia, Poland
a r t i c l e
i n f o
Article history: Available online 3 November 2010 Keywords: Power quality Signal processing Algorithms Digital signal processors Spectral analysis
a b s t r a c t The paper presents original methods and algorithms of measurement of power quality parameters which were implemented in a new measurement device called an estimatoranalyzer of power quality. The main idea was to overcome the shortcomings of current power quality analyzers, namely sensitivity to fast fluctuations of input signal parameters. This was attained with regard to the algorithm’s overall computational complexity. Complementary application of discrete wavelet transform (DWT) with Fourier transforms and chirp z-transform (CZT) was proposed for the aim. Next, a new concept of estimator mode was devised. It consists in multistage signal processing and applied methods depending on input signal character – its frequency and existing waveform distortions – as well as the chosen set of reference parameters. It further improves device performance in terms of time of computation. Description of the estimator-analyzer hardware and results of experimental research will be presented in the second part of the paper. Ó 2010 Published by Elsevier Ltd.
1. Introduction Nowadays, power quality and devices for its monitoring are hotly discussed topic in the field of power systems engineering. The term ‘‘power quality” (PQ) has various meaning but usually it is considered as a combination of voltage and current quality [1,2] or sometimes only the voltage quality [3]. In this paper the term power quality will be used in the latter meaning. It worth adding that the similar digital signal processing techniques can be used for determining both: the voltage and the current quality. The indispensable part of power quality assessment is the measurement of different voltage parameters. The measurement of PQ parameters is usually carried out by dedicated devices called analyzers of power quality. These devices are capable of measurement of many various parameters but arguably the most investigated part of their algorithm is devoted to waveform assessment of the input signal. It consumes a great deal of resources of the measurement device (computational power and RAM memory). Numerous studies have been carried out and ⇑ Tel.: +48 58 690 14 14; fax: +48 58 620 67 01. E-mail address: [email protected] 0263-2241/$ - see front matter Ó 2010 Published by Elsevier Ltd. doi:10.1016/j.measurement.2010.09.049
various methods of estimation of waveform parameters have been proposed. The studies concern methods of evaluation of low frequency disturbances, transient detection and evaluation or data compression. The methods for assessment of low frequency phenomena up to 9 kHz are based on various principles, e.g. methods based on wavelet packet decomposition (WPD) [4], neural networks [5], nonlinear adaptive technique [6], digital filters [7,8] or mathematical morphology [9]. But by far the most often investigated method is still discrete Fourier transform (DFT), especially its fast algorithm (FFT). It is the method recommended in current IEC Standard 61000-4-7 [10]. The FFT algorithm has been implemented in PQ analyzers for many years. In the 1980s, the method was implemented in a device equipped with a Motorola 16-bit 68,000 microprocessor clocked at 4.9152 MHz [11]. The applied constant sampling frequency was 7680 Hz. Unfortunately, FFT has some commonly known downsides. In particular, it is hard achieving both: desired window width and appropriate number of input signal samples. So, synchronization of sampling frequency with actual input signal frequency has been investigated [12,13] and implemented. Additionally different interpolation techniques have been developed [14]. The two techniques are
T. Tarasiuk / Measurement 44 (2011) 238–247
commonly applied in commercial power quality analyzers. Further, signal windowing is still considered [15] and commonly added as an option in digital oscilloscopes. The author’s research has demonstrated that all the methods based on FFT lack versatility and they can lead to dubious results for some real signals registered in electric power systems [17,18]. Simply put, observed maximum errors were above the permissible limit set in the previously mentioned IEC Standard [10]. However, microprocessor technique and especially digital signal processor (DSP) technique has developed rapidly over the last few years. Present DSPs can be clocked at hundreds of MHz (e.g. 32 bit ADSP 21364 maximum clock frequency 333 MHz [16]) and their Harvard architecture and one clock cycle multiplication enable processing of large volumes of data in real time. So, should we still use the above mentioned and recommended in relative standard methods of signal processing based on FFT? The problem of proper waveform assessment increases in difficulty if high frequency phenomena such as transients have to be analyzed. This entails an increase in sampling frequencies and the introduction of additional tools of signal processing, e.g. application of DWT [19–21], mathematical morphology [22] or subtracting the preevent steady-state waveform with sliding window and varying thresholds [23]. In the paper and its second part the original methods and prototype of the new instrument for PQ monitoring is described. It was designed for monitoring of PQ in relatively broad frequency bandwidth and the FFT as basic tool of waveform assessment was discarded. The basic tools of harmonic and interharmonic subgroups measurement are DFT and chirp z-transform (CZT – devised as a fast method of estimation of Fourier coefficients for an arbitrary number of samples [24]), applied for an exact integer number of input signal cycles, independently of its momentary frequency. Therefore, no sampling frequency synchronization or signal windowing is needed. The instrument was called estimator-analyzer of PQ since it has two main modes of operation: analyzer and estimator. Both modes are described in this paper but especially the idea of estimator mode has to be underlined. It is based on the concept of multistage signal processing with varying algorithms, depending on the input signal characteristic and user defined set of reference parameters. In this mode the conformity of PQ parameters with their permissible limits is monitored, with as limited time of computation as possible. The features and merits of methods applied in the original instrument were discussed in the author’s previous papers [17,18,25–30] separately. Therefore, in this paper only the main conclusions from these papers are summarized and related previous papers are recalled. Eventually, the methods were combined as one complex algorithm, appropriate software was developed and the estimator-analyzer was constructed. So, in the paper and its second part the details of the previously unpublished algorithms, hardware and research results of the new instrumentation are presented.
239
2. Standard framework for measuring spectral components in low frequency range IEC Standard 61000-4-7:2002 recommends methods of measurement of harmonics, interharmonics and components above the harmonic frequency range up to 9 kHz [10]. Although other principles of analysis are not excluded, DFT is specified as the reference instrument. According to the standard, new instruments are likely to use DFT, normally using the FFT algorithm. Further, synchronization of measurement time window with actual input signal periods is recommended. The window duration should be 10 periods for 50 Hz systems or 12 periods for 60 Hz systems with rectangular weighting. A Hanning window can be applied only in the case of loss of synchronization. Maximum permissible error of sampling frequency synchronization with actual duration of recommended window width should not exceed ±0.03% [10]. Further, the concept of harmonic subgroup was introduced in the discussed standard [10] in order to account for signal magnitude fluctuations and improve waveform assessment accuracy. The idea consists in calculating the square root of the sum of squares of harmonic component amplitude and amplitudes of two spectral components (interharmonics) directly adjacent to it. The result is considered as the rms value of the harmonic subgroup. The interharmonic centred subgroup is calculated in a similar way on the basis of remaining spectral bins between two harmonic subgroups [10]. However, it was concluded in Refs. [17,18] that the methods recommended in IEC Standard 61000-4-7 [10] could lead to significantly flawed results in some rare but real cases. In the discussed standard method of assessment of the components above the harmonic frequency range up to 9 kHz was proposed. They are to be calculated on the basis of a 100 ms window corresponding to approximately five or six fundamental periods. The FFT with a rectangular window can be used for the measurement, since sampling frequency need not be synchronized with the actual one. But attenuation of the fundamental component is required and it should exceed 55 dB [10]. Finally, the components calculated by FFT are to be grouped in bands of 200 Hz, according to formula [10]:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u bþ100 u X Gb ¼ t C 2f
ð1Þ
f ¼b90ðHzÞ
where Gb – group of bth frequency band, b – centre frequency of considered frequency band, Cf – rms value of spectral line. The above concepts of harmonic subgroups, interharmonic centred subgroups and groups of components in the frequency range from the 50th harmonic up to 9 kHz are utilized in the presented estimator-analyzer of PQ. 3. Estimator-analyzer of power quality 3.1. Hybrid wavelet-Fourier method of spectrum analysis The hybrid wavelet-Fourier method is the basis of all methods and algorithms behind the software of the
240
T. Tarasiuk / Measurement 44 (2011) 238–247
presented estimator-analyzer of PQ. It consists in complementary application of discrete wavelet transform (DWT) and various methods of computation of Fourier series coefficients [17,25,30]. Basically, DWT is applied for transient and notching detection and evaluation. Transients usually contain high frequency components from a few 100 Hz up to a few MHz but their duration is up to 50 ms (in real cases even less) [31]. Unfortunately, the high sampling frequency required for transients analysis would be an enormous burden for the measurement device in the case of low frequency phenomena measurement. Analysis of components in the frequency range up to 9 kHz requires lower sampling frequencies but longer time of observation (window width equals approximately 200 ms for components in the harmonic frequency range and 100 ms above this range up to 9 kHz) [10]. Fortunately, wavelet coefficients after signal decomposition contain all required information about low-frequency components [4,32,33]. Moreover, the number of wavelet coefficients is decimated after each decomposition level [33] and only coefficients representing lower frequencies should be preserved for the whole measurement window for measurement of low frequency parameters. Therefore, if one calculates coefficients of Fourier series on the basis of these wavelet coefficients (instead of original samples of input signal) a huge decrease in the required mathematical operation would be achieved (in practical terms by a factor of 24) [17,25]. Additionally, RAM memory of the used DSP would be saved as well [27]. One more feature of DWT was utilized. Wavelet approximation coefficients which represent low-frequency components with the fundamental one can be treated as a smoothed version of the original input signal. Simply put, the higher frequency components are filtered out. Then, measurement window width can be calculated by analysis of these coefficients’ zero-crossings, improved by simple interpolation. That way, the influence of multiple zero-crossing or higher order harmonics would be minimized [28]. This method of window width estimation was implemented in the described measurement device. 3.2. Estimator-analyzer of power quality – functions overview The basic measurement functions of the prototype estimator-analyzer are as follows [26]: – – – – – – – – – –
–
voltage rms value, frequency, unbalance coefficient, newly introduced distortion band factors (DBF) [27], voltage dip (value and duration), harmonic subgroups up to the 50th harmonic, interharmonic centred subgroups in the range up to the 50th harmonic, group of sub-harmonics, DC component, groups of high frequency components with bandwidths approximately 200 Hz in the range from the 50th harmonic up to 9 kHz, transient and notching detection and assessment,
– power quality temperature factor for assessment of additional temperature rise of windings of induction motors [34]. In all cases, a waveform analysis is carried out for exactly 10 periods of input signal (the prototype was constructed for 50 Hz systems) and only a rectangular window is adopted, independently of input signal character. Voltage dips are estimated on the basis of voltage rms values calculated for one cycle of input signal refreshed each half-cycle Urms(1/2) by the method described in IEC Standard 61000-4-30 [35]. Transients and notching are detected and evaluated on the basis of wavelet coefficients after filtering out low-frequency components by DWT. 3.3. Analyzer mode The discussed measurement device operates in two basic modes: analyzer and estimator. But a significant part of the developed software is common to both modes. The device operates with constant sampling frequency (148917.8 Hz) and first DWT is applied for signal samples in both modes. Daubechies filters with 10 filter coefficients and eight levels of wavelet decomposition are applied. The differences appear after DWT calculation. The estimator mode consists in multistage signal processing with an algorithm that depends on actual characteristics of the input signal. It will be described in the following subsection. In the analyzer mode the device operates on the basis of a constant algorithm. This means that all required PQ parameters are always determined independently of the actual characteristics of the input signal (frequency and waveform distortions). After DWT the following main tasks are carried out: 1. The measurement window width and frequency are calculated on the basis of zero-crossing detection of wavelet approximation coefficients of the seventh decomposition level. These coefficients represents frequency components up to 581.7 Hz (cut-off frequency of respective filters of seventh decomposition level). 2. Signal rms value is calculated on the basis of wavelet coefficients. Details are set out in the following subsection. 3. Partial synthesis of wavelet detail coefficients of third, fourth, fifth, sixth and seventh decomposition levels is carried out. FFT is applied to results of the synthesis in order to estimate frequency components above harmonic frequency range up to 9 kHz. The method enables attenuation of the fundamental component as required by the IEC Standard [10]. Details of the method were provided in Ref. [29]. 4. CZT is applied for wavelet approximation coefficients of the fourth decomposition level and subsequently harmonic subgroups and interharmonic centred subgroups are calculated. The wavelet approximation coefficients of the fourth decomposition level represent frequency components below 4653.7 Hz (cut-off
241
T. Tarasiuk / Measurement 44 (2011) 238–247
frequency of respective filters of fourth decomposition level), so it contains components of harmonic frequency range. The important part of the algorithm is CZT. It is used for calculation of components in the range up to the 50th harmonic, since it enables relatively fast calculation of DFT for an arbitrary number of samples [24]. The block diagram of CZT for the considered aim is depicted in Fig. 1 [24,29,36,37]. It consists of two FFTs and one inverse fast Fourier transform (IFFT). The W constant can be defined as follows [30]:
j 2p W ¼ exp n
ð2Þ
n value represents a theoretical (non-integer) number of wavelet coefficients of the fourth decomposition level and it is calculated by means of the formula:
n¼
Tw fs
ð3Þ
24
where Tw – duration of considered measurement window calculated in the above-described step 1, fs – sampling frequency. Since only components’ magnitudes are estimated, the last multiplication in the CZT algorithm was omitted. The input data sets y1(k) and y2(k) can be calculated by formulae (4) and (5) [29,30].
k2 2 y1ðkÞ ¼ C 4;k W 0 k2 2 W y2ðkÞ ¼ 0 ðkMÞ2 W 2
k ¼ 0; . . . ; N 1 k ¼ N; . . . ; M 1 k ¼ 0; . . . ; N 1 k ¼ N; . . . ; M N 1
ð4Þ
ð5Þ
k ¼ M N; . . . ; M 1
where c4,k – wavelet approximation coefficients of fourth decomposition level, N – number of coefficients calculated for 10 cycles (in this case n from formula (3) means theoretical (non-integer) number of wavelet approximation coefficients of considered decomposition level), M/2 – the nearest power of 2 greater than N. Further, part of the CZT algorithm is used for measurement of frequency components above the harmonic frequency range up to 9 kHz. Since these components may be calculated on the basis of 5 or 6 cycles of the input signal [10] (a window two times shorter than the measurement window for harmonics and interharmonics calculation) the software for parallel FFTs (see Fig. 1) is used for measurement of these components. Data calculated by means of wavelet synthesis during the above-mentioned step 3 for the first half of the basic 10-cycle measurement window replace array y1(k) and the respective data calculated for the second half replace array y2(k). This enables optimization of the two FFTs’ algorithm and improves performance of the algorithm for measurement of the higher frequency components in terms of computational complexity [30]. Finally, a block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels and
Fig. 1. Block diagram of DFT calculation (X(k)-bins of Fourier spectrum) by means of CZT.
main subroutines in analyzer mode is presented in Fig. 2 (excluding transients and dips detection and evaluation). The main difference between the considered device operating in analyzer mode and commercial analyzers consists in the DFT processor, as defined in IEC Standard 61000-4-7 [10]. The measurements of components in the harmonic frequency range are carried out on the basis of the combined DWT–CZT algorithm and FFT is discarded altogether, so synchronization of sampling frequency is not needed. This means that coefficients of Fourier series are always calculated on the basis of an integer number of input signal cycles. So, the basic idea of waveform assessment set out in IEC Standard 61000-4-7 [10] is still preserved. 3.4. Estimator mode The estimator mode is a completely new concept. It is based on the assumption that in many cases only conformity of PQ parameters with their permissible limits specified in related standards or defined by user are necessary. So, it would be useless to determine all PQ parameters and find out that most of them or all are well below their maximum permissible limits. This realization is the basis of the estimator mode algorithm of the discussed measurement instrument. The purpose of measurements in this mode is detection of breaching of the user defined limits, with as reduced number of processor’s instructions as possible. This aim is achieved by multistage signal processing. The algorithm applied at each stage depends on the results of the previous one. After DWT four stages of signal processing are implemented. The algorithm of multistage waveform assessment is graphically presented in Fig. 3 and a block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels and main subroutines in the estimator mode is presented in Fig. 4. 3.4.1. First stage In this stage basic PQ parameters are calculated, namely input signal rms value, its frequency, unbalance coefficient as well as THD coefficient defined as the ratio of rms value remaining after subtracting the fundamental component to the fundamental component [38].
THD ¼
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi U 2rms U 21 U1
100
ð6Þ
where Urms – input voltage rms value, U1 – fundamental component rms value.
242
T. Tarasiuk / Measurement 44 (2011) 238–247
Fig. 2. Block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels (from third to seventh levels of signal wavelet decomposition) and main subroutines in analyzer mode for measurements of low frequency PQ parameters. Dashed line marks flow of information among different subroutines and processor’s Data Memory bank DM required to carry out calculation; h – low-pass filter bank, g – high-pass filter bank, c – wavelet approximation coefficients, d – wavelet detail coefficients.
The value of U1 is calculated by means of ‘‘ordinary” DFT applied for wavelet approximation coefficients of the seventh decomposition level as the input signal. The average number of these coefficients calculated for the duration of the basic measurement window is 27 lower than the number of input signal samples taken by the A/D converter during the same time interval. Obviously, it has a positive impact on the algorithm’s computational complexity. The computation of DFT, CZT or FFT requires appropriate sine and cosine values. It should be mentioned that generally these values are taken from the sine table (4096 values) located in the processor’s data memory bank DM. But THD assessment is an exception. In this case, sine and cosine values are calculated by a min–
max polynomial approximation algorithm [39]. The former method is much faster, but it can lead to unacceptable error in case of low value THDs (below 2% in practical terms due to used formula (6)) and some singular input signal frequencies. This phenomenon was observed during the first tests of the presented device in a ship network. Because increase in volume of the sine table was not possible, other solutions were investigated. A method of calculation of sine and cosine values based on the sine table and interpolation was tested, but the results were mixed at best. Finally, the above-mentioned method of min–max polynomial approximation algorithm led to satisfactory results and it was used for calculation of only U1 in estimator mode.
T. Tarasiuk / Measurement 44 (2011) 238–247
243
Fig. 3. Algorithm of multistage waveform distortion assessment in low frequency range up to 9 kHz in estimator mode, PV – permissible value for chosen set of parameters – 0.5% in the case of EN Standard 50160 [40] or 3% in the case of IEC Standard 60092-101 [38], ncDFT – relative number of mathematical operations for DFT and actual input signal waveform distortion, ncCZT – relative number of mathematical operations for CZT and actual input signal waveform distortion.
If the calculated THD value is lower than the maximum permissible value (PV) of any required higher frequency
component the whole analysis for the considered measurement window is finished and further analysis stages are
244
T. Tarasiuk / Measurement 44 (2011) 238–247
Fig. 4. Block diagram of utilization of wavelet coefficients from respective wavelet decomposition levels (from third to seventh) and main subroutines in estimator mode for measurements of low frequency PQ parameters (third to seventh levels of signal wavelet decomposition); dashed line marks flow of information among different subroutines and processor’s Data Memory bank DM required to carry out calculation; h – low-pass filter bank, g – high-pass filter bank, c – wavelet approximation coefficients, d – wavelet detail coefficients.
omitted. PV are chosen by the user and it depends on the chosen set of permissible values of PQ parameters. There are two predefined sets of these values included in the RAM memory of the used DSP. The first set is based on EN Standard 50160 [40]. The second set is based on IEC Standard 60092-101 [38]. In the case of the former standard [40] PV was assumed to be equal to 0.5%. In the case of the latter standard [38] PV was assumed to be equal to 3.0%.
3.4.2. Second stage This stage consists in rough waveform assessment by means of distortion band factors (DBF). The concept of these new factors was described in previous papers by the author, most importantly in Ref. [27]. In shorthand, DBFs are calculated by means of the formula:
DBF miþ2 ¼
U rmsðfmi fmði1Þ Þ U1
100½%
ð7Þ
245
T. Tarasiuk / Measurement 44 (2011) 238–247
where U rmsðfmi fmði1Þ Þ – rms value of voltage components of (fmi – fm(i1)) frequency band, U1 – rms value of the fundamental component. The rms value of signal components for each frequency band can be calculated by means of wavelet detail coefficients for the respective decomposition level as follows [27]:
U rmsðfmi fmði1Þ Þ
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u i 1 u 1 NX ¼t i ðdmi;k Þ2 2 ni k¼0
ð8Þ
where dmi,k – wavelet detail coefficients of respective decomposition level, Ni – number of wavelet coefficients of considered decomposition level calculated for a time interval of basic measurement window, n – theoretical (non-integer) number of wavelet coefficients in an analysed window calculated by means of appropriately modified formula (3), i – calculation order of the respective decomposition level (i 2 1. . . m), m – number of decomposition levels (m = 8 for the described device). The last DBF for the lowest frequency band is calculated differently, by the formula [27]:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u m X u DBF 1 ¼ tTHD2 DBF 2miþ2
ð9Þ
i¼1
DBF1 comprises components up to a cut-off frequency of low-pass filter of the last decomposition level without the fundamental component. DBFs can be calculated without an additional requirement for computational power of the measurement device. The considered signal rms value was calculated concurrently with these DBFs by the formula [27,41]:
U rms
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi " # u N m1 m1 i 1 u X NX X1 1 2 t 2 ¼ cm7;k þ ðdmi;k Þ 2m1 nm1 k¼0 i¼1 k¼0 ð10Þ
where cm7,k – wavelet approximation coefficients of seventh decomposition level, m = 8. One should take into account that the applied number of wavelet coefficients N in formula (8) is decimated for each consecutive layer of wavelet decomposition. This means that calculation of DBFs together with the considered signal rms value by means of wavelet coefficients requires a roughly similar number of mathematical operations as the calculation of only the signal rms value on the basis of the original samples of the signal. So, implementation of the method does not increase the computational complexity of the algorithm for PQ estimation on the one hand. But it increases the scope of information about the signal behaviour on the other hand [27]. (In fact, DBFs are calculated independently of THD value. If THD < PV further analysis based on DBFs is skipped.) DBFs can be treated as a rough measure of waveform distortions in the whole frequency band that is determined by the cut-off frequency of the anti-aliasing filter. More information about possibilities of DBFs utilization can be found in Ref. [27]. Since each DBF represents information about waveform distortion in the respective frequency band, it can be used
in a similar way as the THD factor. Namely, the content of the respective frequency band should be analysed pffiffifficarefully only if the product of corresponding DBF and 2 exceeds the permissible limit of any higher frequency component from the analyzed frequencyp band. The necessity of multiffiffiffi plication value of DBF with 2 results from the shape of wavelet filters response [27]. The permissible limit depends on the chosen set of reference values and pffiffiffi analysed frequency band. If all products of all DBFs and 2 are lower than the respective permissible limits, the whole analysis is finished and the next steps are omitted. Otherwise a further analysis is required. In the case of chosen standards all permissible limits equal PV defined for THD assessment in the first stage. Firstly, if values of DBF6 or DBF5 corresponding with wavelet detail coefficients of the fourth and fifth decomposition levels pffiffiffi (frequencies 2326.8–9307.4 Hz) are greater than PV/ 2, a subroutine for calculation of components above the harmonic frequency range up to 9 kHz is called. (The same subroutine is always called in analyzer mode.) Secondly, an affected frequency band in the harmonic frequency range is determined by analysis of respective DBFs associated with the harmonic frequency range. These DBFs, corresponding frequency bands and number of decomposition levels are presented in Table 1 (assuming sampling frequency 148917.8 Hz). By DBFs analysis a set of wavelet coefficients for detailed analysis of harmonics or interharmonic subgroups is defined. It should be mentioned that in many cases only lower order harmonics are presented in real power system signals. It diminishes algorithm computational complexity due to decimation of the wavelet coefficient after each decomposition level. For example, DBF2 is calculated on the basis of wavelet detail coefficients of the eighth pffiffiffi decomposition level. If DBF2 has a value above PV/ 2, a computation of Fourier series coefficients would be carried out on the basis of wavelet approximation coefficients of the sixth decomposition level. This way a possible magnitude damping of Fourier coefficients (due to the shape of wavelet filters’ frequency response) would be less than 0.01 and aliasing of components above 581.7 Hz would be minimized as well. Fortunately, the number of wavelet approximation coefficients of the sixth decomposition level is decimated by a factor of 26 in comparison to the number of input signal samples. As a result a decrease in computational complexity would be achieved. 3.4.3. Third stage In this stage FFT is performed. Input data set for the procedure consists of wavelet coefficients defined in the Table 1 DBFs of harmonic frequency range and related frequency bands. DBF
Frequency band (Hz)
Number of wavelet decomposition level
DBF5 DBF4 DBF3 DBF2 DBF1
2326.8–4653.7 1163.4–2326.8 581.7–1163.4 290.9–581.7 Below 290.9 without fundamental component
5 6 7 8 –
246
T. Tarasiuk / Measurement 44 (2011) 238–247
previous stage by the method described above. The sampling frequency (148917.8 Hz) was chosen in such a way that for all possible frequencies of input signal in steady-state conditions (50 Hz ±5%) the number of input data for CZT in analyzer mode would equal 4096. At the same time the number of wavelet coefficients calculated for 11 cycles of input signal would be the power of two for each decomposition level, if the actual input signal frequency was equal to the rated one (in fact 49.99 Hz; exact synchronization for 50 Hz was not possible) [29]. Finally, such a synchronization condition is checked (measured frequency within 49.99 Hz ±0.03% limit). If it is fulfilled, an FFT subroutine is called and a rectangular window is used. Then, harmonic subgroups and interharmonic centred subgroups are calculated and compared with their permissible values (specified in the relevant standard). Then the analysis is finished. If synchronization was not obtained, Hanning weighting is applied prior to FFT analysis and the fourth stage of analysis is carried out. 3.4.4. Fourth stage Because application of a Hanning window can lead to flawed results in some specific but real signals [17,18] further processing is necessary. In this stage a final method of signal processing is chosen. In Ref. [42] a windowed (Blackman) FFT was proposed for coarse spectrum analysis in order to localize frequency bands with significant spectral energy. Next, windowed CZT is applied for a medium search in the specified frequency band. CZT is applied for evaluation of the exact frequencies of components for detailed analysis. Subsequently, evaluation of exact magnitudes and phases is carried out by DFT [42]. This seems too complicated and hardly practical in the case of significantly distorted signals. So, the author proposes a modified procedure. First the relative computational complexity of DFT is estimated and compared with the relative computational complexity of CZT. The calculation of relative computational complexity takes into account the exact wavelet decomposition levels which contain components of interest as well as results of preliminary assessment by FFT with Hanning window. It can be carried out relatively fast but requires prior assumption of the appropriate threshold. Relative computational complexity of CZT depends on utilized wavelet coefficients (decomposition level), whereas relative computational complexity of DFT depends on input signal character – actual harmonic or interharmonic content or more precisely their frequencies in relation to wavelet decomposition levels (DFT can be applied for wavelet coefficients of various decomposition levels depending on order of actually calculated frequency bin [25]). Finally, a precise calculation of specified harmonic subgroups and interharmonics centre subgroups is carried out on the basis of appropriate wavelet coefficients by means of DFT or CZT, whichever would be faster. It should be added that final results are always calculated on the basis of a rectangular window, independently of actual signal frequency and used method (FFT, DFT or CZT) and despite constant sampling frequency. Moreover,
in practical cases, it should be faster than in analyzer mode (detailed analysis of time of computation will be presented in Part II). Comparison of Figs. 2 and 4 clearly indicates that most parts of developed subroutines are common in both modes. So, introduction of the estimator mode leads to only a slight increase in the size of final software and required fast RAM memory of the DSP used.
4. Final remarks The solution described in this paper is based on various methods for PQ assessment. These methods were integrated as one complete tool for profound PQ analysis. The resulting algorithm conforms to basic concepts of PQ measurement laid in the respective IEC Standard [10], although it differs from the methods directly recommended in the standard and commonly used by developers of current commercial PQ analyzers. During development of the algorithm three aims were to be achieved: increasing the sampling frequency as high as possible, increasing the algorithm immunity for relatively fast fluctuations of input signal parameters as well as saving resources of measurement device, such as computational power and RAM. The aims were obtained by complementary application of various tools of signal processing: DWT, FFT, DFT and CZT. This lead to increased performance of signal processing algorithm, in comparison to application of these tools separately. Its is also worth stressing the concept of the adaptation of processing algorithm to the actual character of the input signal (its frequency and waveform) by multistage signal processing. It was called the estimator mode and leads to a further improvement in terms of computational complexity. The introduction of this mode increases functionality of measurement device. It enables: increasing of sampling frequency or adding additional functions of measurement instrument, e.g. extending number of monitored parameters. Further, the mode can be useful in handheld devices supplied by build-in batteries. Simply put, if external power is switched-off the DSP core instruction rate can be decreased (lower computational power), which leads to energy saving and prolongs battery life. The exact gains of the implementing of the author’s proposal depend on the used hardware platform. The described solution was devised chiefly for DSP based instruments. Such a instrument, its structure and detailed results of its research will be described in Part II. Nevertheless, the idea can be also implemented in PC based instruments with standard data acquisition boards or instruments with combination of DSP and FPGA for signal processing. So, its possible applications are broader than the exemplary device described in Part II. References [1] M. Bollen, What is power quality?, Electric Power Systems Research 66 (2003) 5–14 [2] D. Katsaprakakis, D. Christakis, A. Zervos, S. Voutsinas, A Powerquality measure, IEEE Transactions on Power Delivery 23 (2) (2008) 553–561.
T. Tarasiuk / Measurement 44 (2011) 238–247 [3] R. Dugan, M. McGranaghan, H. Beaty, Electrical Power Systems Quality, McGraw-Hill Companies, 2002. [4] L. Eren, M. Unal, M.J. Devaney, Harmonic analysis via Wavelet packet decomposition using special elliptic half-band filters, IEEE Transactions on Instrumentation & Measurement 56 (6) (2007) 2289–2293. [5] L.L. Lai, W.L. Chan, C.T. Tse, A.T.P. So, Real-time frequency and harmonic evaluation using artificial neural networks, IEEE Transactions on Power Delivery 14 (1) (1999) 52–59. [6] D.M. McNamara, A.K. Ziarani, T.H. Ortmeyer, A new technique of measurement of nonstationary harmonics, IEEE Transactions on Power Delivery 22 (1) (2007) 387–395. [7] J.J. Tomic´, M.D. Kušljevic´, V.V. Vujicˇic´, A new power system digital harmonic analyzer, IEEE Transactions on Power Delivery 22 (2) (2007) 772–780. [8] L. Ferrigno, C. Landi, M. Laracca, C. Luongo, A study on the feasibility and effectiveness of digital filters approach for power quality monitoring in compliance with IEC 61000-4-7, in: 15th IMEKO TC4 Conference, Iasi, Romania, 2007. [9] P.M. Ramos, T. Radil, F.M. Janeiro, A. Cruz Serra, DSP based power quality analyser using new signal processing algorithms for detection and classification of disturbances in a single-phase power system, Metrology and Measurement Systems 14 (4) (2007) 483–494. [10] IEC Standard 61000-4-7, General Guide on Harmonics and Interharmonics Measurements for Power Supply Systems and Equipment Connected Thereto, 2002. [11] A.N. Mortensen, G.L. Johnson, A power system digital harmonic analyzer, IEEE Transactions on Instrumentation & Measurement 37 (4) (1988) 537–540. [12] M. Aiello, A. Cataliotti, V. Cosentino, S. Nuccio, Synchronization techniques for power quality instruments, IEEE Transactions on Instrumentation & Measurement 56 (5) (2007) 1511–1519. [13] A. Cataliotti, V. Cosentino, S. Nuccio, A phase-locked loop synchronisation of power quality instruments in the presence of stationary and transient disturbances, IEEE Transactions on Instrumentation & Measurement 56 (6) (2007) 2232–2239. [14] M. Sedlácˇek, M. Titeˇra, Interpolations in frequency and time domains used in FFT spectrum analysis, Measurement 23 (3) (1998) 185–193. [15] A. Testa, D. Gallo, R. Langella, On the processing of harmonics and interharmonics: using Hanning Window in Standard Framework, IEEE Transactions on Power Delivery 19 (1) (2004) 28–34. [16] Analog Devices, SHARC Processors. ADSP 21362/ADSP 21363/ADSP 21364/ADSP 21365/ADSP 21366. Data Sheet, 2007. [17] T. Tarasiuk, Hybrid wavelet-Fourier method for harmonics and harmonic subgroups measurement – case study, IEEE Transactions on Power Delivery 22 (1) (2007) 4–17. [18] T. Tarasiuk, Comparative study of various methods of DFT calculation in the wake of IEC Standard 61000-4-7, IEEE Transactions on Instrumentation & Measurement 58 (10) (2009) 3666–3677. [19] O. Poisson, P. Rioual, M. Meunier, Detection and measurement of power quality disturbances using wavelet transform, IEEE Transactions on Power Delivery 15 (3) (2000) 1039–1044. [20] L. Angrisani, P. Daponte, M. D’Apuzzo, A. Testa, A measurement method based on the wavelet transform for power quality analysis, IEEE Transactions on Power Delivery 13 (4) (1998) 990–998. [21] S. Santoso, W.M. Grady, E.J. Powers, J. Lamoree, S.C. Bhatt, Characterization of distribution power quality events with Fourier and wavelet transforms, IEEE Transactions on Power Delivery 15 (1) (2000) 247–254.
247
[22] S. Ouyang, J. Wang, A new morphology method for enhancing power quality monitoring system, Electrical Power and Energy Systems 29 (2007) 121–128. [23] D. Castaldo, D. Gallo, C. Landi, A. Testa, A digital instrument for nonstationary disturbance analysis in power lines, IEEE Transactions on Instrumentation & Measurement 53 (5) (2004) 1353–1361. [24] L.L. Rabiner, R. Shafer, C. Rader, The chirp z-transform, IEEE Transactions on Audio & Electroacoustic 17 (2) (1969) 86–92. [25] T. Tarasiuk, Hybrid wavelet-Fourier spectrum analysis, IEEE Transactions on Power Delivery 19 (3) (2004) 957–964. [26] T. Tarasiuk, A few remarks about assessment methods of electric power quality on ships – present state and further development, Measurement 42 (8) (2009) 1153–1163. [27] T. Tarasiuk, The method based on original DBFs for fast estimation of waveform distortions in ship systems – case study, IEEE Transactions on Instrumentation & Measurement 57 (5) (2008) 1041–1050. [28] T. Tarasiuk, Wavelet coefficients for window width and subsequent harmonic estimation – case study, Measurement 41 (3) (2008) 284– 293. [29] T. Tarasiuk, Method, algorithm and device for estimation of components above the harmonic frequency range up to 9 kHz, Measurement 44 (1) (2011) 219–229. [30] T. Tarasiuk, Application of CZT transform for spectrum analysis in systems with varying frequency, in: IEEE Int. Workshop on Intelligent Data Acquisition and Advanced Computing Systems: Technology and Applications, 2007, Germany, pp. 329–334. [31] IEEE Standard 1159-1995. IEEE Recommended Practice for Monitoring Electric Power Quality. [32] W. Mikhael, A. Ramaswamy, An efficient representation of nonstationary signals using mixed-transforms with applications to speech, IEEE Transactions on Circuits and Systems – II Analog and Digital Signal Processing 42 (6) (1995) 393–401. [33] S. Mallat, A theory for multiresolution signal decomposition: the wavelet representation, IEEE Transactions on Pattern Analysis and Machine Intelligence 2 (7) (1989) 674–693. [34] P. Gnacin´ski, J. Mindykowski, T. Tarasiuk, A new concept of the power quality temperature factor and its experimental verification, IEEE Transactions on Instrumentation & Measurement 57 (8) (2008) 1651–1660. [35] IEC Standard 61000-4-30, Testing and Measurement Techniques – Power Quality Measurement Methods, 2003. [36] C. Zang, L. Dai, H. Zheng, J. He, Using frequency zoom technology to realize high precision and adaptive frequency measurement for power system, in: IEEE Int. Conference POWERCON, Singapore, 2004, pp. 155–159. [37] T. Zielin´ski, Digital Signal Processing, WKŁ, Poland, 2005 (in Polish). [38] IEC Standard 60092-101, Electrical Installations in Ships. Definitions and General Requirement, 2002. [39] Analog Devices, ADSP-21000 Family. Applications Handbook, vol. 1, 1996. [40] EN Standard 50160, Voltage Characteristics of Electricity Supplied by Public Distribution Systems, 1999. [41] Weon-Ki Yoon, M. Devaney, Power measurement using the wavelet transform, in: IEEE Instrumentation and Measurement Technology Conference, 1998, pp. 801–806. [42] K.J. Runtz, D. Hack, A multistage DFT – FFT – CZT approach for accurate efficient analysis of sparsely distributed spectra, in: IEEE Canadian Conference on Electrical & Computing Engineering, 2002, pp. 127–132.