|Estimations are free. For more information,
please send a mail
|Power Spectral Density computation (Spectral Analysis)
MicroJob Package Deal Service PD001A/B(Copyright Cygnus Research International (Apr 20, 2015))
Page 6 of PD001A/B User Guide
Table of Content
Appendix 1; About numerical error generated by computers.
(4-2) Taper ratio of Tukey window; 0~100(%) Default setting is 10%
(4-2-1) Why we need to apply a window function
This assumption connects the end point of our data to the start point of our data as shown in red ellipses in Figure 4-14(b) and if there is a 'cliff' caused by the difference of values at the end point and the start point of data as shown in Figure 4-14(c), it will raise PSD across entire frequency range. The purpose of applying a window function is to remove these cliffs so that this artificial increase of PSD would be reduced. This artificial increase of PSD is spectral leakage and we will describe more about the relationship between spectral leakage and window functions in Appendix 2.
(4-2-2) What does window function do
(4-2-3) Why Hanning window is not necessarily a good solution
Figure 4-16(a) shows PSD of this data with and without windowing. Since this graph is too complicated to make comparisons, we applied 3-bin FDS and the results are shown in Figure 4-16(b). This graph shows that the Hanning windowed PSD is consistently smaller than the others within Band B, and, to the lesser degree, within Band A. Hanning windowed PSD divided by non-windowed PSD shown in Figure 4-16(c) indicates this feature more clearly.
The left part of Figure 4-17 ((a), (c) and (e)) shows filtered time series for bands A, B and C computed by the method described in (2-1-3-5) using non-windowed amplitude spectrum and the right part of Figure 4-17 ((b),(d) and (f)) shows the result of windowing for each of these components. By using equation (3), we can express our data as a summation of 4 components each covering different frequency range as
This is not really important for the argument here, but summation of bands A and B reproduces approximate profile of data as shown in Figure 4-18. Therefore, Band A and Band B are very important in this particular case.
(4-2-4) Why we provide non-windowed and Tukey windowed result for PD001A
One of the methods to reduce dumping near the both ends of data while still using Hanning window is to use so called Welch's method. At this moment we do not have an intention to implement this method in this package deal service but we describe it briefly here.
First, we cut the data in the middle along the dashed line in Figure 4-19(a). This will produce two data segments. Then we chop off the data along the dashed lines in Figure 4-19(b) and pick up the central part of the data.
Now we have three data segments of equal length. Then, we apply Hanning window to each of these data segments individually as shown in Figure 4-20 and compute PSD of each. This will give us three sets of PSD.
After that, we compute average of PSD for each bin (ensemble average) from these three sets of PSD. Figure 4-21(a) shows the weight, which is a value of a Hanning window function, multiplied to data. This graph also shows the coverage of individual data segment. The trick is that by allowing overlapping of data segments in this way, we can reduce the portion where excessive dumping occurs. If we add three curves in Figure 4-21(a) together we will get 'Overlapped window' in Figure 4-21(b) (blue line). The red line in this graph is the weight of an ordinary (non-overlapped) Hanning window function for a comparison. When we use Welch's method, some values of data are used more than once. Half of the data are used twice in this example.
The blue line in Figure 4-22 shows PSD computed by this method and the purple line in this graph shows non-windowed PSD. We applied 5-bin FDS to the non-windowed PSD to make comparison easier. The reason why we specifically choose 5-bin for FDS bin-width here is as follows. Firstly, the data length of individual data segment when we applied Welch's method is half of original data length. Then, frequency width of each bin is twice as wide as frequency width of each bin when we did not apply Welch's method. Secondly, the width of implicit FDS caused by Hanning window is three-bin as we will explain in (4-2-6). Since we cannot chose even number for the bin-width of explicit FDS, we decided to have 5-bin instead of 6-bin for FDS bin-width.
The red line in this graph shows non-overlapped (ordinary) Hanning windowed PSD. The purple and the red lines are the same as those in Figure 4-16(b) except for the bin-width of FDS. This graph shows that the reductions of PSD within bands A and B are considerably reduced if we use Welch's method. Thus, we can conclude that the reductions of PSD within bands A and B when we applied ordinary Hanning window are artificial and Welch's method will fix a problem. However, there is a price to pay; the frequency resolution is worsened and the lowest non-zero frequency of PSD (third raw of our data file) increases since the data length of each segment is shorter than the original data.
(4-2-6) Implicit FDS
First of all, one-sided windowed continuous PSD is the square of absolute value of Fourier transform of data multiplied by a window function divided by length (time) of data and multiplied by 2. This is very mouthful but if we use a pseudo equation, it can be expressed as
Instead of using a continuous form of Fourier transform to compute continuous PSD, we usually use a discrete form of Fourier transform to compute PSD at discrete frequencies. The discrete Fourier transform of data, xm where the subscript m is the index ranging from 0 to N-1, at the frequency of k_th (k is 0 to N/2 where N is the number of data) bin is
Now, Fourier transform of data xm multiplied by a window function wm is the convolution of the Fourier transform of xm and the Fourier transform of wm according to the Convolution theorem. The important point here is that this convolution of the Fourier transform of xm and the Fourier transform of wm is the same as a moving average of the Fourier transform of xm with the non-uniform weights which are the Fourier transform of wm except that the summation of weights is not equal to 1.0. Let Fourier transform of a window function wm at the frequency of n_th bin be
Here we intentionally use 'n' instead of 'k' as an index because we want to use Wr(n) and Wi(n) as weights in the convolution we will do later. In case of Hanning window, Wr(0)=T/2, Wr(1)=Wr(-1)= - T/4 but others are all zero.
Using Convolution theorem, equation (18), equation (20) and the fact that all the Wr(n) and Wi(n) except for Wr(0), Wr(1) and Wr(-1) in equation (20) are zero, Fourier transform of data multiplied by a Hanning window function at the frequency of k_th bin can be expressed as
The coefficient 1/T in the right hand side of above equation is equal to the frequency width of bins. When k is 0, we replace k-1 with N/2 and when k is N/2, we replace k+1 with 0. Using the fact that Wr(-1)=Wr(1)= - T/4 and Wr(0)=T/2, equation (21) can be simplified as
If we consider "Wr" as weights, the right hand side of equation (22) is simple 3-point weighted moving average (again, in the strict sense it is not an average because the summation of weights is not equal to 1.0) of Xr+iXi and this is why we call this an implicit 'FDS'. Hanning windowed PSD of k_th bin is, replacing Xr(k) by Br(k) and Xi(k) by Bi(k) in equation (19),
This equation shows why the width of the peak when we applied Hanning window becomes 3 bins when data length is A (Please see Figure 2-6, Figure 4-6(a) or Figure 4-10(c)). Signal at k_th bin affects PSD of k-1, k and k+1_th bins. Although we usually do not call this widening of signal spectral leakage, the basic principles of this widening and spectral leakage are the same as we will describe in Appendix 2. Another important aspect of equation (23) is that it is no longer a simple 3-point weighted moving average unlike equation (22) and we will describe about this topic next. If we use other window functions, we generally do not get such a simple equation as above because Fourier transform of a window function (equation (20)) generally produces many non-zero Wr(n) and Wi(n).
(4-2-7) Implicit FD'Smoothing' Really?
(4-2-8) Reduction of PSD and amplitude due to the window
We multiply PSD by CP and amplitude spectrum by CA to correct reduction. The important point about this correction is that this method only uses information about a window function.
(4-2-9) Does correction of reduction of PSD work well?
In case of PSD, it is slightly more complicated because we can think of two criteria to judge the performance of correction. The first one is the same as the one we used for amplitude correction. Another one is that how well equation (1) holds. Since windowing reduces magnitude of data, magnitude of resultant PSD is also reduced. Therefore, without correction the right hand side of equation (1) becomes smaller than the left hand side of equation (1). Here we assume that the "data" in the left hand side of equation (1) is not windowed data but the original one since we usually are not interested in variance of windowed (thus reduced) data. Equation (1) is very important for PSD.
Figure 4-24(a) shows how PSD of the bin where signal exist (we called it signal bin before) varies as we add more data when data contains only a single frequency sinusoid but nothing else. The data we used here is exactly the same as the one shown in Figure 2-4(a). We computed PSD without detrending to avoid effect of artificial inclination we described in (4-1-3) and the reduction of PSD due to the window is corrected using the method described above. The reason why PSD increases as we add more data is that the frequency bandwidth of the bin decreases as we add more data (Figure 2-4(f)) while signal amplitude itself does not change. We show this graph only to demonstrate that PSD increases as we add more data even if we keep everything else the same but there is nothing wrong about it.
Since graph like Figure 4-24(a) is not convenient to interpret results, we computed variance of the signal bin by multiplying PSD of the signal bin by the frequency bandwidth. Let us call this value, signal bin variance. Figure 4-24(b) shows signal bin variance when we applied Tukey window (blue), Hanning window (red) and when we did not apply any window (purple). The black solid line is the variance of data and this value becomes 0.5 when data length is such that the data is terminated exactly at the end of the cycle of the signal (frequency match; vertical black dashed line). This graph shows that the Hanning windowed signal bin variance is consistently smaller than others. This might make us conclude that the correction of reduction is wrong based on the first criterion.
Now let us use the second criterion. As we described before, Hanning window makes width of peak of PSD wider. Figure 4-24(c) shows the results of summation of variance of all of the bins except zero frequency bin (left hand side of equation (1)). The result of non-windowed PSD (purple) in this graph becomes the same as data variance (black line in Figure 4-24(b)) by equation (1) and, thus, we did not draw data variance. This graph shows that the summed variance of Hanning windowed PSD is not particularly smaller than other cases. Figure 4-24(d) shows results of summation of variance of 3 consecutive bins, center of which is the signal bin. This graph shows that 3-bin summation of variance of Hanning windowed PSD is much closer to the variance of data and its variation is much smaller. If we look at this graph carefully we can see that the 3-bin summation of variance of Hanning windowed PSD and data variance both becomes 0.5 where frequency match occurs (vertical black dashed line). Therefore, we can say that the correction of reduction of PSD is done well based on the second criterion.
This apparent contradiction of the results of these two criteria comes from the fact that implicit FDS of window functions would modify frequency resolution of resultant PSD. In other word, we have to use larger frequency bandwidth to compute signal bin variance in Figure 4-24(b). We describe more about this issue in (4-3-3). Our conclusion here is that the correction applied to PSD works well but we need to consider the actual frequency resolution when we use a value of PSD.
We applied these corrections to all of the examples we have shown in this document. As these graphs show, correction generally works well both for PSD and amplitude spectrum. However if characteristic of data varies a lot from one moment to another as in the case of Figure 4-17, this method might not work well for specific frequency ranges.
(4-2-10) Tukey window; how taper ratio affects its characteristics
Figure 4-25 shows a value (weight) of 50% tapered Tukey window function as an example. If we applied this window to the data, the first 25% and the last 25% of data will be tapered (reduced) but 50% of the data in the middle will be unaffected (weight equals 1.0). The characteristic of Tukey window varies from the one equivalent to that of non-window by setting taper ratio to 0% to the one equivalent to that of Hanning window by setting taper ratio to 100% as we mentioned before.
Figure 4-26(a) shows part of PSD when we applied 10% tapered (blue), 50% tapered (red) and 90% tapered (purple) Tukey window as well as part of PSD when we did not apply any window (dark green; equivalent to 0% tapered Tukey window) and when we applied Hanning window (black; equivalent to 100% tapered Tukey window). The data we used for these computations contains only a single sinusoid, frequency of which is 100. This data is the same as the one shown in Figure 2-4(a) except for frequency. The data length is adjusted so that the frequency match is achieved (data length A in Figure 2-4(a)). In this case, non-windowed PSD (dark green) shows that only one bin, frequency of which is exactly the same as signal frequency, contains entire variations of data and PSDs of all other bins are zero within computational accuracy (blue line in Figure 4-5(a)). In other words, there is no spectral leakage and this is true even if we applied any window function. However, this graph shows Tukey windowed PSDs have many bumps labeled side lobes. The reason why they exist is that the Fourier transform of Tukey window produces many non-zero Wr(n)s and Wi(n)s in equation (20) unlike Hanning window, and those non-zero coefficients spread the effect of the signal to an entire frequency range. As we will describe in Appendix 2, cause of this side lobes is practically the same as the cause of spectral leakage. The magnitudes of side lobes are symmetrical with respect to the signal bin and they decrease as taper ratio increases. As we mentioned before, we will see side lobes even in the case when we applied Hanning window if we compute continuous PSD (Appendix 2).
Figure 4-26(b) shows an expanded view of amplitude spectra near the signal bin. This part is often called a main lobe. We drew dashed lines connecting bins adjacent to the signal bin to show how the 'local' width of peak varies. This graph shows that the values of the bins next to the signal bin steadily increase (dashed lines shift upward) and top to middle part of 'mountain' becomes wider as taper ratio increases. However, if we look at the width of the 'slope' part of this 'mountain', 50% tapered Tukey window gives us the widest 'slope' in this graph. Figure 4-26(a) shows that overall, or global, sharpness of the peak increases if we use a larger taper ratio but Figure 4-26(b) shows that 'local' sharpness decresases if we only look at few bins near the signal bin.
Figure 4-26(c) shows the same as Figure 4-26(a) does except that we have changed the data length so that the difference between frequency of the signal bin and the actual signal frequency becomes the largest (data length C in Figure 2-4(a)). In this case we expect large spectral leakage occurs. Non-windowed PSD (dark green) shows the largest PSD across the entire frequency range and this is caused by spectral leakage. Tukey windowed PSDs in this graph are not that much different from those in Figure 4-26(a) in term of amplitude although the shapes of side lobes are different. Hanning windowed PSD does not have bumps but the magnitude is close to the 90% tapered Tukey windowed PSD as we expected. Figure 4-26(d) shows the same as Figure 4-26(b) does except for the data length. Amplitude of signal bin approaches closer to the actual signal amplitude (1.0) as taper ratio increases but the amplitude of a bin to the left of the signal bin also increases and the top of the 'mountain' becomes flatter. Variations of amplitudes at 'slope' are somewhat complicated as before.
In short, Figure 4-26 shows that smaller taper ratio gives us higher frequency resolution (sharper main lobe) but larger taper ratio gives us better PSD (closer to actual PSD).
Figure 4-27(a) shows PSD of real data when we applied 10% tapered (blue), 50% tapered (red) and 90% tapered (purple) Tukey window as well as PSD when we did not apply any window (dark green) and when we applied Hanning window (black). The time series of this data is shown in Figure 4-15(a). Since this graph is too complicated to make comparisons, we applied a 3-bin FDS and the results are shown in Figure 4-27(b). The differences among PSDs are most prominent within Band B as expected because large variations within band B are concentrated near both ends of data (Figure 4-17). However, in general, we feel differences are relatively small for casual use of PSD. These graphs suggest that, although the side lobes of Tukey window in Figure 4-26 are annoying, their actual effects may be negligible for majority of practical applications.
Figure 4-28(b) shows how the maximum and the minimum amplitudes change as taper ratio increases. Here the maximum amplitude is the amplitude when the frequency of the signal bin matches actual signal frequency (VH in Figure 4-28(a)) and the minimum amplitude is the amplitude when difference between frequency of the signal bin and actual signal frequency becomes the largest (VL in Figure 4-28(a)). As we mentioned before, Tukey window becomes Hanning window when taper ratio is 100%. When we deal with real world data, it is usually difficult to achieve frequency matches. Therefore, it would be better if the difference between VH and VL were small. This graph shows that larger taper ratio gives smaller difference between VH and VL. Thus, again, larger taper ratio would likely give us more accurate PSD.
Finally we show how PSD and amplitude vary when we change taper ratio from 10% to 90% with 10% increment in Figure 4-29. The first row ((a), (b) and (c)) shows the shapes of Tukey window. The second row ((d), (e) and (f)) shows PSDs similar to the one shown in Figure 4-26(a). The third row ((g), (h) and (i)) shows amplitude spectra similar to the one shown in Figure 4-26(b). The fourth row ((j), (k) and (l)) shows amplitude spectra similar to the one shown in Figure 4-26(d). The fifth row ((m), (n) and (o)) shows PSDs similar to the one shown in Figure 4-27(a). We grouped results of three different taper ratios together. The taper ratios and corresponding color codes are shown in the first row for each column. We hope this document will provide enough information so that our customers can determine proper taper ratio.