|Estimations are free. For more information,
please send a mail
Power Spectral Density computation (Spectral Analysis)
Page 5 of PD001A/B User Guide
Table of Content
Appendix 1; About numerical error generated by computers.
(4) Guide for setting proper selectable parameters
(4-1) Detrending and removing average; Yes or No? Default setting is Yes.
(4-1-1) What is detrending? Why do we have to do that?
(4-1-2) Example of detrending
and the intercept is the vertical position of the purple line at t=0, which is one point left (outside) of the left edge of this graph. After computing slope and intercept, we compute linear approximation of data, namely (Slope x tm + Intercept) at each data point. The purple line in this graph is this linear trend. Finally, we subtract this linear trend from data and the result is shown in Figure 4-1(b).
(4-1-3) When detrending becomes harmful
The blue curve in Figure 4-2(a) shows first one period plus one more point of test data we used in Figure 2-4(a) and this corresponds to V(tm) in equation (13). This is just a minor detail but the reason we added one extra point which is actually the first point of the second period is that the end point of the first period is not the point where blue curve crosses zero but one point before that occurs. We want to have zero-to-zero data here and since our data is discrete (not continuous), we need to have one extra point. Now, our first guess of the linear line approximation, a linear trend, of blue curve is black line; zero slope and zero intercept.
Unfortunately our rather innocent guess turned out to be wrong and an actual linear trend of blue curve using equation (13) is the purple line in this graph. The red curve in this graph shows detrended data. The cause of the problem we describe in this section is the fact that, by detrending a pure sinusoid, we actually add linear 'inclination' component (in the graph, it is subtracted but that is equivalent to adding a negative of computed linear trend) to a pure sinusoid. Please note that PSD computation is based on the premise that our data can be described by the summation of pure sinusoids of discrete frequencies and a constant (zero frequency component) as shown by Equation (3).
Making data a bit longer does not solve this problem. Even if we have 8 period lengths of data, corresponding to data length A in Figure 2-4(b), we still have non-zero 'inclination' as shown in Figure 4-2(b). This can be mathematically explained by letting sampling interval very small and then approximating summation in equation (13) by integral. By doing so, slope for the sinusoidal curve of amplitude G, frequency F and phase P, that starts from time=0 and ends at time=T, becomes
We skip all the algebraic detail of getting this equation because it is a tedious book-keeping job and we do not think showing the process of getting this equation will match the purpose of this document.
Figure 4-3(a) shows how the slope varies as we add more data when amplitude G is 1.0 and frequency F is 1.0. The unit of horizontal axis of this graph is the period of sinusoid (=1/F). The data lengths A, B and C indicated by the vertical short blue lines are the same as those in Figure 2-4(b). If we set start point of our data half a period ahead, result will become a mirror image of this graph with respect to the horizontal black line drawn at 0.
This indicates that slope is almost always negative and it only becomes approximately zero when data length is extremely long or approximately n plus half a period where n is any integer greater than 1. These are clearly seen in Figure 4-3(a). This equation also shows that the magnitude of slope increases as the amplitude of sinusoid increases and decreases as frequency of sinusoid and/or data length increases. With a fixed data length T, slope decrease as frequency of sinusoid increases because there will be more and more oscillations within the data.
It is ironical enough that we will add 'inclination' to a pure sinusoid by detrending it but the most ironical part of this result is that both slope and intercept of 'inclination' become large when we achieve frequency match (data length A). The result of PSD computation is supposed to be the most accurate when we achieve the frequency match as shown before.
Adding 'inclination' to a pure sinusoid would cause two effects on PSD. Let us call 'a pure sinusoid', signal. The first one is addition of PSD of 'inclination' to the PSD of the signal. If data length is A (end of 8th period; T=8/F), then only a single bin, frequency of which is equal to frequency of the signal (we call it frequency match before), has non-zero PSD if we do not detrend (Figure 2-5(a)) and this is the real PSD of this signal. PSD due to the added 'inclination' at k_th bin for data length A is from equation (14),
Here, the amplitude of the signal is equal to 1.0 (G=1.0). We give bin number one to the lowest non-zero frequency bin (third raw of our product file). PSD of zero frequency bin (k=0; second raw of our product file) is zero if we detrend data. This equation shows that PSD due to added 'inclination' increases if frequency of the signal decreases. PSD of the bin where signal exist in this particular case (data length=A) will be simply adding the result of equation (17) to the real PSD of the signal itself which is the variance of signal divided by the frequency bandwidth of bin.
Equation (17) also shows that if we draw this PSD on the logarithmically scaled graph we will see a straight-sloped line since frequency, f, of k_th bin is k/T. Taking the logarithm of both sides of equation (17),
Figure 4-4(a) shows PSD obtained from above equations and Figure 4-4(b) shows PSD numerically computed by our program of this package deal service when data length is A. The inclination of the straight part of PSD in Figure 4-4(b) is about -1.9994341, which is very close to the value in above equation (-2.0). This linear inclination of -2.0 appearing in a logarithmically scaled graph is a signature of the existence of a linear inclination in our data. Now suppose there is actually a linear trend in our data and we do not detrend it before computing PSD. Then, we will see a similar signature in a resultant graph. When we cannot achieve frequency match, things are going to be slightly more complicated because the addition of the effect of 'inclination' must be done in the complex domain, which is one step before finishing the PSD computation. We also will have spectral leakages of signal across the entire frequency range in that case and that effect might overwhelm the effect of an artificial inclination as we will show next.
Before showing what the PSD looks like when weak signals coexist with a strong signal, we will show the differences between detrended and non-detrended PSD when we have only a single frequency strong signal for three different data lengths and how window functions affect those differences first. Figures 4-5(a) shows PSD when data length is A and no window is used. The blue line in this graph is the result when data is not detrended and average is not removed. This PSD is the same as the one shown in Figure 2-8(a). The red line is the result when data is detrended and this PSD is the same as the one shown in Figure 4-4(b) except for the vertical scaling. The blue line of this graph is what we want to see as the PSD of this data. The effect of detrending is annoying, but it would not obscure the signal itself. Figures 4-5(b) and (c) show the same as Figure 4-5(a) but when data length is B and C (please see Figure 4-3(a)), respectively. These two graphs show that the effect of detrending is much smaller than the effect of spectral leakage and we cannot see red lines at all because blue lines are practically drawn over red lines. Therefore, it does not matter whether we detrend data or not for these two cases. This is what we meant by condition (4). Figures 4-5(d), (e) and (f) show the same as Figures 4-5(a), (b) and (c) except that we applied a Tukey window before computing PSD. We can see small differences between non-detrended and detrended PSD when data length is A or B but it is apparent that the implicit FDS, effect of which appears as many bumps, plays a great roll determining PSD at frequencies other than the signal frequency. If we applied a Hanning window, Figures 4-5(g), (h) and (i) show that the effects of both detrending and spectral leakage are reduced. In any case, if we cannot match the frequency, it will be difficult to see any adverse effect of detrending due to the spectral leakage when signal frequency is low. When signal frequency is high we might be able to see the effect of detrending even if we have a spectral leakage as we will show later.
The cause of the difference we can see at the lowest frequency part of PSD when data length is B and C, marked by red circles in Figures 4-5(h) and (i), is slightly complicated. Figure 4-6(a), (b) and (c) show the same PSD as in Figure 4-5(g), (h) and (i) except that we removed average from non-detrended data before computing PSD. PSD at the lowest non-zero frequency (left side edge of this graph) becomes larger than the PSD at one point right of it when data length is B and C (red circles; Figure 4-6 (b) and (c)). This is caused by removal of average before computing PSD combined with the effect of an implicit FDS of a Hanning window and this increase of PSD is artificial. We will show that an implicit FDS of a Hanning window occasionally un-smooth PSD in (4-2-7). Please remember that detrending also removes the average of the data.
We do not see similar behavior when we applied a Tukey window with a default setting of 10% taper. This is because the weights of implicit FDS caused by a Tukey window at a point next to the central point is much smaller than those caused by a Hanning window as shown by the fact that the width of peak is narrower (peak itself is sharper; Please see Figure 2-6). This fact causes the PSD to which a Tukey window is applied at the lowest non-zero frequency not to be affected that much by value at zero frequency, which is practically zero because we removed the average.
Removing an average before computing PSD is often considered to be a least requirement if we do not apply detrending to get a better estimate of PSD at a low frequency range. The reason behind it is that if we have a relatively large average and thus have a large PSD at zero frequency, PSD at lowest few bins are badly affected by zero frequency PSD because of implicit or explicit FDS. Thus, it may sound contradictory that we have artificially large PSD at the lowest non-zero frequency precisely because we removed average. However, the particular behavior we have shown here is due to relatively rarely occurring un-smoothing effect of an implicit FDS of a Hanning window and, thus, it is not universal. Also, our test data is rather special and when we deal with real world data, we also have signals and noises at many frequencies and it will be rather rare to see an adverse effect of removing an average. Thus, we still recommend removing an average if we do not detrend.
Figure 4-7(a) shows how the amplitude of signal bin varies as we add more data in the case when we detrend data before computing PSD. Comparing this graph with Figure 2-4(b) (No detrending), we can see that amplitude of signal bin approaches to the actual signal amplitude much slower. This can be seen more clearly if we compare Figure 4-7(b) with Figure 2-4(c). Both of these graphs show amplitude of signal bin only at positions where frequency match occurs. Figure 4-7(c) shows that variance of data deviates considerably from 0.5 when data length is shorter than 2.0 unlike variance shown in Figure 2-4(d). Although we can see these differences, the adverse effect on amplitude estimation caused by detrending is relatively minor as long as data length is longer than a few periods of signal.
Figure 4-8 shows PSD when there are signals at multiple frequencies to demonstrate that the problem we described so far is not limited to the case when we have a single frequency signal. The data we used here is the noiseless data we used in Figure 2-22(red lines) and we did not use any window function. In this graph we can see that there are differences between detrended and non-detrended PSD at lower frequency range even when we have spectral leakages (data lengths B and C). This is because effect of spectral leakage is relatively small at far left of these graphs simply because frequency of that part of graph is far enough from the signal frequencies.
Now the question is what would happen if we have weak signals coexisting with a strong sinusoidal signal (Let us call it S1) as we have been showing so far. To demonstrate this, we generated test data by adding part of data shown in Figure 1-1(a) to S1. In the following, let us call these added weak signals SR. Figure 4-9(a) shows time series of SR. We adjusted amplitude of SR so that the variance ratio of SR to S1 is 0.0001. Significant digit of our computation is large enough to keep SR faithfully. The unit of horizontal axis (Data Length) of this graph is not the same as Figure 1-1(a) but we chose this one for convenience. Figure 4-9(b) shows PSD of SR. We detrended SR and used a Tukey window before computing PSD. We made the unit of horizontal axis (Frequency) of this graph consistent with previous graphs. There are two peaks labeled S2 and S3 in this graph. We use these two peaks as targets to detect because variations represented by these peaks exist relatively steadily throughout the entire record of SR.
First, let us demonstrate when the frequency of S1 is relatively low. The frequency of S1 is 5 cycles per 1 time unit of PSD (Data Length 5 is equivalent to 1 time unit). Figure 4-9(c) shows synthesized data, namely S1+SR. We cannot see the existence of SR in this graph at all because SR is much weaker than S1. In term of amplitude instead of variance, S1 is about 3 million (3.0E+6) times larger than S3. Please note that the range of vertical axis of Figure 4-9(c) is much larger than the range of vertical axis of Figure 4-9(a). We will show PSD when data length is A and C indicated by vertical red lines in Figure 4-9(c). Data length A is the case when frequency of S1 matches one of the frequencies of PSD and the data length C is the case when difference between frequency of S1 and the nearest frequency of PSD becomes the largest as before. We omitted data length B to save some spaces.
Blue lines of upper three graphs of Figure 4-10 show PSD without detrending and without removing average when data length is A. This is the special case (we called frequency match before) that we do not have a spectral leakage of S1 at all. However, when we need to deal with real world data, it is quite likely that strong signal is composed of several different frequency components and it will be very difficult, if not impossible, to achieve frequency match like this. Therefore, we believe this case is good only for demonstration. Having mentioned that, the blue and red asterisks (*) in these graphs show the locations of S2 and S3. The black lines in Figure 4-10 show PSD of SR alone. For the sake of comparison, we did not detrend SR for these graphs and we moved these black lines upward to make comparison easier. The blue line in Figure 4-10(a) shows that PSD computation without detrending and without using any window did do a superb job. It shows that the S1 affects PSD only at the frequency where S1 exists. When we applied a Tukey window or a Hanning window, Figure 4-10(b) and (c) show that PSD within a certain frequency range around S1 is dominated by S1 due to the implicit FDS and it is difficult to see PSD of SR within that frequency range.
The lower three graphs of Figure 4-10 show PSD when we detrend data. Figure 4-10(d) shows that if we do not use any window function, we are no longer able to see anything except S1 and the effect of added 'inclination' caused by S1. The reason why it is practically impossible to see PSD of SR when we detrend is not that PSD of SR is buried under the PSD of the added 'inclination' generated by detrending (please see Figure 4-4(b)). PSD in this case is close to the summation of PSD of S1 and SR. There are two possibilities. The first one is truncation error. Showing the entire significant digits, PSD of the 'inclination' at the frequency of S3 is 0.1180623344177745E-9. There are 16 digits below decimal point. If PSD of S3 (Figure 4-9(b)) is less than 0.0000000000000001E-9 (15 zeros below decimal point), then adding PSD of S3 to the PSD of the 'inclination' would not change total PSD at all because of the limitation of the size of memory assigned to store that value. In our example this does not happen because the PSD of the peak of S3 is 0.00000681348403258001E-9 (5 zeros below decimal point). The second reason, and this is what happened here, is that the PSD of S3 is too small to make any detectable visual difference. Adding PSD of S3 to PSD of the 'inclination' would not change first 5 digits.
Figure 4-10(b) and (e) show that the difference between detrended and non-detrended PSD is relatively small when we used a Tukey window. The implicit FDS applied to S1 has a very strong influence on PSD across the wide frequency range and S2 is no longer detectable. We will describe more about this implicit FDS in the next section. We can see the hint of the effect of added 'inclination' at very low frequency range. If we use a Hanning window instead, detrended PSD (Figure 4-10(f)) shows a straight-sloped line at low frequency range caused by added 'inclination'. The effect of an implicit FDS is limited to one bin up and one bin down from the bin where S1 exists and we can get more information about SR than the case when we used a Tukey window.
Figure 4-11 shows PSD when data length is C (Figure 4-9(c)) as an example when we have a spectral leakage of S1 across the entire frequency range and this is more likely the case we would encounter. If we do not apply any window function, the effect of spectral leakage overwhelms the effect of added 'inclination' as expected from Figure 4-5(c) and it is not possible to see anything except S1. Detrend or not does not make any practical differences. If we apply a window function, effect of a spectral leakage is reduced but the effect of an implicit FDS is added and the difference between detrended and non-detrended PSD is difficult to see except at the lowest frequency. This difference at the lowest frequency is the same kind of difference we described using Figure 4-5 and Figure 4-6 (red circle). When frequency of S1 is higher we can clearly see the difference between detrended and non-detrended PSD especially when we used a window function as we will show next.
Now, let us change the frequency of S1, but keep anything else the same. Figure 4-12 shows PSD when data length is A. S1 is now located at 1000. Again, PSD without applying any window function and without detrending is the best. S1 in this graph appears to be sharper (width is narrower) than S1 in Figure 4-10(a) but this is only a result of visual effect caused by logarithmically scaled horizontal axis. If we compare Figure 4-12(d) with Figure 4-10(d), we can see that the PSD due to added 'inclination' becomes much smaller. This is expected from equation (17), which shows that PSD due to added 'inclination' is inversely proportional to the frequency of S1. We can even detect S2 (barely) and S3 in the detrended PSD. We will describe how the frequency of S1 affects detectability of SR in Appendix 4.
If we compare Figure 4-12(b) with Figure 4-10(b), we might think that effect of an implicit FDS when we applied a Tukey window is considerably narrowed but, again, this is a result of visual effect of logarithmically scaled axis. The width and the magnitude of the effect of an implicit FDS do not depend on the frequency of S1.
Figure 4-13 shows PSD when data length is C. If we do not apply any window, we will practically miss SR almost completely although non-detrended PSD shows hint of it at low frequency range. If we applied a window function, non-detrended PSD gives us more information about SR than detrended PSD. Differences between non-detrended and detrended PSD are most prominent at low frequency range partially because SR has a larger PSD at that frequency range. When we applied a Tukey window, SR deforms some of the typical bumps caused by S1 and it is slightly deceptive.
As a summary, if it is important to detect SR, these results suggest that using a Hanning window without detrending might be the best choice although there are cases when a Hanning window might cause a problem for a quite different reason as we will describe in (4-2-3). However, if our data shows an obvious trend, we have to detrend data. The problem here is that the amplitude of S1 is too large compared with the amplitude of SR and ordinarily small effect of added 'inclination' due to detrending of S1 is enormous for SR. If amplitudes of most of the signals in SR are comparable to amplitude of S1, we will not see the effect of added 'inclination' just like blue lines of Figure 2-22. If our strong signals are limited within a relatively narrow frequency range, it might be better to apply a digital filter to our data and suppress strong signals before computing PSD. That method also reduces the effect of an implicit FDS and spectral leakage caused by strong signals. If our customers wish us to apply a filter before computing PSD, please, let us know. These results seem to show an awful effect of detrending but in reality we think the majority of us would never encounter data like this.