|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 7 of PD001A/B User Guide
Table of Content
Appendix 1; About numerical error generated by computers.
(4-3) Bin-width of FDS; 1,3,5,7,... (Odd integer); Default setting is None, 3 and 7
(4-3-1) FDS works well for PSD
(4-3-2) FDS does NOT work well for amplitude spectrum
On the other hand, when frequency of the signal bin matches actual signal frequency (data length A in Figure 2-4), the amplitudes of Bin 1 and Bin 3 become zero. Figure 4-32(a) shows time series of data in this case. Amplitude of Bin 2 is almost equal to the actual signal amplitude (1.0) and amplitudes of all other bins are zero as shown in Figure 4-32(b). Thus, simple 3-bin average is about 1/3 in this case.
Now, let us show why Hanning windowed amplitude spectrum becomes vastly larger than non-windowed amplitude spectrum when we applied FDS. For the sake of simplicity, let us chose data length to be data length A, so that frequency of the signal bin matches actual signal frequency and there are no spectral leakages. Figure 4-32(c) shows time series of data when we applied a Hanning window. If we applied a Hanning window, width of peak becomes wider as we described in (4-2-6) and now amplitudes of bin 1 and bin 3 are no longer zero as shown in Figure 4-32(d). The correction we described in (4-2-8) brings amplitude of bin 2 to the value almost equal to actual signal amplitude, 1.0, but it also amplifies amplitudes of Bin 1 and Bin 3 as shown in Figure 4-32(e). As a result, simple 3-bin average becomes about 2/3, which is about two times larger than the value when window is not applied. This is why Hanning windowed amplitude spectrum becomes vastly larger than non-windowed amplitude spectrum when we applied FDS.
Thus, in general, FDS does not produce meaningful result for amplitude spectrum and this is why we do not write amplitude in our product file when we applied FDS.
(4-3-3) Frequency resolution when we applied FDS
The primarily reason why we do not discard any of Po is that it would be convenient when we want to import our product file into applications for drawing graphs and/or doing further analyses. Another reason is that we can use intermediate Po such as Po(1) and Po(2) as interpolated Po between Po(1) and Po(6).
(4-3-4) Alternative method of FDS
(4-4) Percentage of confidence interval; any positive value larger than 0.0 and smaller than 100.0; Default is 95.0
Appendix 1; About numerical error generated by computers.
Appendix 2; Spectral leakage of power spectrum density function and window functions.
In (4-2-1) we described that the 'cliff' at the boundary generated by the assumption of repetition of data (Figure 4-14(b)) is a cause of spectral leakage. The assumption we described there is consistent with the method we and many other people use to compute PSD and the explanation of the cause of spectral leakage is reasonable. However, it only gives us an idea how spectral leakage is generated but we cannot evaluate how large generated spectral leakage would be. Also we will have a problem explaining why we still have spectral leakage even after we applied a window function to remove 'cliffs'.
In this appendix we first introduce alternative way of looking at our limited length data. After that we show how spectral leakages are generated and how large they would be. Finally we describe how to reduce spectral leakages. The root cause of spectral leakages is the fact that the length of our data is limited.
(1) Alternative way of looking at our limited length data
We do not assume repetition of data described in (4-2-1). Then we apply a window function w(n) shown in Figure A2-1(b) to this infinitely long data. This window function is called a Boxcar window function due to its shape. The result is "Our data 2" shown in Figure A2-1(c). At this point all we need to know about "Original data" is the values between t=0 and T and we know them because our data covers that duration.
We can construct "Our data 2" in Figure A2-1(c) in a different way. We start from our limited length data x'(n) shown in Figure A2-2(a).
This is the exact copy of the one we showed in Figure 4-14(a). Then we add infinitely long "Additional data" z(n) shown in Figure A2-2(b). The value of z(n) is zero for any n. This results that Fourier transform of z(n) is also zero for any frequency and, thus, adding z(n) to "Our data" x'(n) would not modify spectral characteristic of x'(n) at all. The result of adding z(n) to x'(n) is "Our data 2" y(n) shown in Figure A2-2(c). This "Our data 2" is exactly the same as the one shown in Figure A2-1(c). As we described before, we can add any number of data, value of which is exactly zero, to data we want to analyze to adjust data length so that we can achieve frequency match. We called it "zero padding". What we did here is the extreme form of zero padding. Since data length of "Our data 2" is infinite, frequency bandwidth of bins becomes zero (=1/infinite), and, thus, we cannot evaluate power spectrum "DENSITY" at discrete frequencies anymore.
The reason why we introduced these two methods to generate identical infinitely long data ("Our data 2") is that we will describe spectral leakage as a difference of spectral characteristic of infinitely long data shown in Figure A2-1(a) and "Our data 2" shown in Figure A2-1(c) (and Figure A2-2(c)). This approach makes sense only because spectral characteristic of "Our data 2" is the same as the spectral characteristic of "Our data" shown in Figure A2-2(a) because "Our data" is what we usually have, i.e. a limited length data. In other word, we consider limited length "Our data" as a result of the application of a Boxcar window function to an infinitely long data in this appendix. This leads us that if we study how changing the width of a Boxcar window function affect spectral of resultant windowed data (Figure A2-1), then we would know how changing the data length of "Our data" would affect spectral (Figure A2-2). Thus we will show how changing the width of a Boxcar window function affects spectral of resultant windowed data in this appendix.
Instead of using PSD we introduce similar quantity named PS for an infinitely long data as follows. First we extend the range of summation in equation (18) as
For a Boxcar windowed data, "Our data 2 (y)" in Figure A2-1(c), equation (A2-1) becomes
The transition from the one marked by blue underline to the one marked by red underline is possible because wm equals zero if m is less than 0 or greater than N-1 and equals 1 otherwise. Please note that rightmost part (red underline) of equation (A2-2) is the same as equation (18) except for the frequency. Since we can compute Yr and Yi at any frequency, we can compute them at frequencies where we computed PSD in section (2). In that case, the rightmost part of equation (A2-2) will be exactly the same as equation (18).
This quantity PS is very much similar to PSD (equation (19)). One of the apparent differences is that there is no multiplication factor 2.0 in equation (A2-3). This is because we want to define PS both at positive and negative frequencies. Another apparent difference is that equation (A2-3) is not divided by T but this is simply because T is infinite for an infinitely long data. Aside from these apparent differences, the most important difference between PS and PSD is that PS is continuous in frequency and this becomes very convenient to understand why spectral leakage occurs.
Actual values of PS of "Our data 2 (y)" in figures A2-1(c) and A2-2(c) divided by T (equals width of a Boxcar window function) and multiplied by 2 at frequencies where we ordinarily compute PSD are exactly the same as PSD of "Our data (x')" in Figure A2-2(a). This is because equation (A2-2) will be exactly the same as equation (18) in that case as we described before.
(2) How spectral leakages are generated and how large they would be
Equation (A2-5) shows that the Fourier transform of x(t) is not continuous but non-zero only at two specific frequencies. Using the result of this and the convolution theorem we described in (4-2-6), we can evaluate middle part (blue underline) of equation (A2-2) as follows.
The last part of equation (A2-6) (red underline) show that Fourier transform of Boxcar windowed x(t) is quite different from Fourier transform of x(t) shown by equation (A2-5). This difference is the spectral leakage. Decomposing W in equation (A2-6) into a real part (Wr) and an imaginary part (Wi) and substituting them into equation (A2-3), PS of the Boxcar windowed data is,
To show how the spectral leakage occurs by showing the convolution process graphically we have chosen two different cases. Figure A2-3(a) shows time series plot of data (Equation (A2-4)) we use.
Each data is shown as a small dot in this graph and straight blue lines connect dots. The first case is when window-width is 40 (data points). This case is equivalent to the data length A shown in Figure 2-4 in the sense that frequency match is achieved. From this graph we might get an impression that the second cycle is not completed in this case but the zero value point next to the end point of this case is actually the start point of the next cycle and does not belong to the second cycle. Another case is when window-width is 50 and this case is close to the data length C in Figure 2-4. We intentionally made these window-widths rather narrow so that following graphs are easy to see. In the following we set amplitude of x(t) (A in equation (A2-4)) 1.0. The sampling interval is 0.05, which results that fundamental frequency range is -10 to 10.
Before showing how the convolution process works we show PS for infinitely long but Boxcar windowed data similar to "Our data 2" in Figure A2-1(c) and PSD for limited length data similar to "Our data" in Figure A2-2(a) first. The red line in Figure A2-3(b) shows PS of Boxcar windowed infinitely long data x(t) for the case A. Since data length is infinite, we can evaluate continuous PS. We omitted PS in the negative frequency range because PS is symmetrical with respect to the zero frequency. The vertical blue arrow in this graph is PSD of limited length data multiplied by T and divided by 2 and we call it modified PSD. It is not an accidental that the tip of this arrow touches PS (red line). The modified PSD (blue arrow) and PS (red line) at the same frequency agree as we described before. Since this is the case when frequency match is achieved, PSDs of bins are zero except for a signal bin as we showed in Figure 2-5(a) (Purple line; "No window" case). Thus, there is only one blue arrow in this graph and the signal frequency shown by a short green solid line matches frequency of that arrow. The vertical dashed blue lines indicate frequencies of other bins. PS at those frequency is zero from the second line of equation (A2-7) because both f1 and f2 of equation (A2-6) become integer multiples of 1/T at those frequencies.
Figure A2-3(c) shows PS and modified PSD for the case B. Since PSD of all the bins are non-zero due to the spectral leakages, we omitted vertical dashed blue lines in this graph. The solid blue line (labeled as "Envelope") connecting tips of the arrows shows that the result shown in this graph is similar to the PSD shown in Figure 2-6(c) (Purple line). Again, modified PSD and PS at the same frequency agree. Figures A2-3(b) and (c) indicate that spectral leakage always occurs but we do not see it on graphs of PSD when we achieved frequency match because the effect of spectral leakage is zero at frequencies where we ordinarily compute PSD for that case. It is slightly difficult to see but the PS of the main lobe is not symmetrical with respect to the center of the main lobe. Due to this asymmetry the bin of the maximum amplitude may not be a signal bin when we cannot achieve frequency match as we described in (2-1-3-3).
Now let us show how the spectral leakage occurs by showing the convolution process graphically. We picked up bins marked by purple solid triangles in figures A2-3(b) and (c) for this purpose. The upper panel of Figure A2-4 shows the Fourier transform of an infinitely long data x and the lower panel of Figure A2-4 shows the Fourier transform of a Boxcar window function.
These are graphical presentations of equations (A2-5) and (A2-7), and the fundamental frequency ranges (horizontal width in the graph) of these two are the same because sampling interval of infinitely long data x(t) and Boxcar window function w(t) must be the same. There is no real part component in the upper graph of Figure A2-4. To calculate PS at frequency fT (purple solid triangle) using convolution, first thing we do is flipping lower graph left and right. This is already done so that frequency (horizontal axis) of the lower graph decreases from left to right unlike upper graph (red circle). The next thing we do is moving this lower graph horizontally so that the horizontal position of the center of lower graph marked by a black solid triangle matches the horizontal position of fT in the upper graph. The result of this manipulation is shown in Figure A2-5.
We omitted unnecessary part of graphs. Then, we draw lines vertically from the upper graph to the lower graph at frequencies f0 and -f0 (dashed lines H and L). Since the frequency difference between H and C is fT+f0=f1 and the frequency difference between L and C is fT-f0=f2, frequencies of these two dashed lines in the lower graph is f1 and f2 (Please see equation (A2-6)). If f1 goes beyond of the edge shown in Figure A2-4, we use periodicity of the Fourier transform as we show later,
As a fourth step, we read values of blue (real part) and red (imaginary part) lines at places where dashed lines H and L intersect those lines. Then, we calculate X multiplied by W as
In case of window-width A shown in Figure A2-5, Wr(f1), Wi(f1), Wr(f2) and Wi(f2) are all zero as we described before. Thus PS at fT is zero from equation (A2-8). Likewise, PS at frequencies where we ordinarily compute PSD except for a signal bin becomes zero and we do not see spectral leakages at those frequencies. However, if we compute PS at any other frequency, we will pick up a signal no matter how far away fT is from f0 (signal frequency). Figure A2-6 shows the case when window-width is B.
In this case, none of Wr(f1), Wi(f1), Wr(f2) and Wi(f2) is zero, and as a result, PSD becomes non-zero even if fT is far away from f0. This is how the spectral leakages occur.
Figure A2-7 shows convolution process for general cases when the Fourier transform of infinitely long data ((a):X) is more complicated than the one we have shown so far.
In this graph we show segment A covering fundamental frequency range and segment B which has a equal frequency width (horizontal width in this graph) of segment A. The contents of these two segments are identical due to the periodicity of X we described before. Let us assume we want to compute PS at frequency fT (purple triangle). We shift horizontal position of the graph of the Fourier transform of a Boxcar window function ((b):W) after flipping it left and right as before. For the convolution we chose segment E of X. The frequency width of segment E is the same as frequency width of segment A as well as frequency width of W. The segment E lacks the segment C of segment A at left but the segment D of segment B added at right compensates this missing segment. Again, due to the periodicity of Fourier transform contents of segments C and D are identical. In this way, the extent of convolution covers entire fundamental frequency range of X no matter where fT exists. We multiply segment E of X by W in the same way we did in figures A2-5 and A2-6; multiply two values at the same horizontal position together. The result of this multiplication is shown in Figure A2-7(c). Then we integrate this result to obtain Yr(fT) and Yi(fT) to compute PS at frequency fT using equation (A2-3). These graphs clearly show that the convolution is equivalent to the weighted moving average (in the strict sense it is not an average because the integral of weight is not equal to 1.0) of X with the weight W and the width of averaging is equal to the fundamental frequency range. As a result, PS at any frequency fT would pick up all the signals on X except for those on frequencies f that satisfy fT - f = +/- k/T (second line of equation (A2-7)).
(3) How to reduce spectral leakages
Figures A2-8(a) and (b) show Fourier transform of a Boxcar window function and a Hanning window function, respectively.
While we can see blue and red lines clearly at entire frequency range in Figure A2-8(a), it is difficult to see those lines beyond of +/- 2.0 in Figure A2-8(b). Figure A2-8(c) shows PS of these two windows. Please note that we used a logarithmic scaling for this graph unlike others. From this graph we can see that a Hanning window would considerable reduce spectral leakages and reduction rate would progressively increase as the difference between fT and f0 becomes larger. This graph also shows that, although applying a Hanning window removes 'cliff' described in (4-2-1), we will still have spectral leakages.
Another thing we can see in this graph is that while the PS of a Boxcar window function at frequencies marked by blue asterisks is zero, PS of a Hanning window function is not zero at those frequencies. Also, PS of a Hanning window function at zero frequency (black dashed vertical line) is smaller than the PS of a Boxcar window function at that frequency. These characteristics result the widening of a signal on the graph of PSD when we applied a Hanning window function as we described in (4-2-6).
These results suggest that it would be better if we apply a Hanning window before computing PSD. However, a Hanning window function and many other window functions considerably dump significant portion of data. In certain cases, applying these window functions may cause a problem as we described in (4-2-3).
Appendix 3; Why we see a single peak instead of twin-peaks when data length is short.
When time series is short relative to the period of S1, the maximum of the coefficient of Fourier transform of S1 is at zero frequency because frequency resolution is so large that f1 is included within zero frequency bin. Then, convolution would not modify frequency of S0 and the peak of PSD of our data will be at f0. At this stage, data is practically S0 except that amplitude is not constant, but is slowly varying. Please note that this data length is short relative to the period of S1 but it is fairly long relative to the period of SL and SH as shown in Figure 2-11.
As we increase data length, frequency resolution of the Fourier transform decreases and f1 leaves from zero frequency bin but moves toward higher frequency. Then, convolution of this result and the Fourier transform of S0 generates twin-peaks center of which is at f0. It is not our intention to write any educational material, so we will stop here. We will describe more about convolution in (4-2-6).
Appendix 4; How the frequency of S1 affects the detectability of SR
END of User Guide of PD001A (Copyright: Cygnus Research International, Apr. 21, 2014)