|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
As an example of an alternative window function we use a Hanning window function here. The Fourier transform (discrete) of a Hanning window function is
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).
(4) Alternative method to reduce spectral leakages (Not supported in this package deal service)
We will start showing that overall general characteristics of PSD of data that contain various signals and noises such as SR computed by MEM are quite similar to the one computed by PD001A/B. Figure A2-9(a) and (b) show PSDs of SR computed by PD001A/B (blue and red lines) and by MEM (magenta lines). We moved blue and red lines upward to make comparison easier. The data, SR, we used here are the same as the one shown in Figure 4-9(a). The number of data is 8514 and this is the case when spectral leakage caused by S1 would be close to the maximum if we add S1 and use PD001A/B to compute PSD (similar to "Window-width B" in Figure A2-3 or "data length C" in Figure 2-4(b)). We did not detrend data for the sake of consistency. There are two peaks labeled S2 and S3 in these graphs. We use these two peaks as targets to detect as before because variations represented by these peaks are tidal components and they exist relatively steadily throughout the entire record of SR.
The blue line in Figure A2-9(a) shows the case when we did not apply any window function and it is the same as the black line in Figure 4-13(a) except that we used linear horizontal axis this time and we moved the line vertically. The red line in Figure A2-9(a) shows the case when we applied a Hanning window and it is the same as the black line shown in Figure 4-13(c). The difference between Figure A2-9(a) and A2-9(b) is that the results shown in Figure A2-9(a) are the cases when we did not apply FDS and the results shown in Figure A2-9(b) are the cases when we applied FDS. FDS has an effect of smoothing PSD as we described in (4-3).
The magenta lines in Figure A2-9(a) and A2-9(b) show PSDs computed by MEM. There is a parameter called autoregression order (AR order) when we compute PSD by MEM and it affects practical frequency resolution among others. We adjusted AR orders so that resultants PSDs have similar practical frequency resolutions to others in Figure A2-9. We can compute PSD at arbitrary frequencies if we use MEM but we computed PSD at the same frequencies as others (blue and red lines) for the sake of comparison.
The point we would like to make here is that Figure A2-9(a) and A2-9(b) show that overall general characteristics of PSDs of SR computed by these different methods are quite similar. After all we are trying to compute same thing, PSD. However, if we take a close look at these graphs we will see the differences. For example, PSD computed by MEM suggests S3 might include signals of several slightly different frequencies, which is quite possible because S3 corresponds to semi-diurnal tide, and semi-diurnal tide is known to have several components at slightly different frequencies. Likewise, some of these differences of PSDs of SR might become very important for serious spectral analyses but we will ignore these differences in this Appendix.
The blue lines in Figure A2-10 show PSDs when S1 is added to SR and the black lines in Figure A2-10 show PSDs when there is only SR. Figure A2-10(a) shows PSDs computed by PD001A/B when we did not apply any window (equivalent to applying a box car window) and Figure A2-10(b) shows PSDs computed by PD001A/B when we applied a Hanning window. Figure A2-10(c) shows PSDs computed by MEM. We did not detrend data for these computations and we moved black lines downward to make comparison easier. The blue line in Figure A2-10(a) is the same as the blue line in Figure 4-13(a) and the blue line in Figure A2-10(b) is the same as the blue line in Figure 4-13(c) except that we used linear horizontal axis this time. The black lines in Figure A2-10(a), (b) and (c) are the same as the blue, red and the magenta lines in Figure A2-9(a) except that we moved some lines vertically.
Figure A2-10(a) shows that spectral leakage of S1 practically makes it impossible to detect SR if we did not apply any window. Figure A2-10(b) shows that spectral leakage of S1 is reduced considerably by a Hanning window but still existence of it is obvious. MEM, by principle, produces no spectral leakage by windowing and, Figure A2-10(c) clearly shows it.
Now, let us take a look at S1 instead of SR. The PSDs of S1 shown in Figure A2-10(a) and A2-10(b) are reasonably close to the actual PSD of S1. On the other hand, PSD of S1 shown in Figure A2-10(c) is considerably smaller (about 1/500000000) than the actual PSD of S1 although PSDs of S2, S3 and others of SR are comparable to those in Figure A2-10(b). This is mostly because we did not compute PSD at the frequency of S1 (800) although S1 is a constant amplitude, constant frequency, pure sinusoid without any phase variations. To make this point clearer we recomputed PSD with a much smaller frequency interval. The blue (S1+SR) and the black (SR alone) "x" marks in Figure A2-11 show the results of these recomputations. Please note that the frequency range of this graph is only between 785 and 815.
The vertical red lines in this graph indicate frequencies where we computed PSD in Figure A2-9(c) and the vertical black line at the center of this graph indicates frequency of S1. One of the frequencies where we computed PSD this time exactly matches the frequency of S1 and the blue "x" mark at the center of the blue circle labeled S1 shows PSD at the frequency of S1. This value is reasonably close to PSD of S1. Theoretically, we should not be able to see S1 except at the exact frequency of S1 but we can see S1 in Figure A2-10(c) because MEM gives non-zero frequency width to S1 (S1 is not a "line" but a "mountain") as shown in Figure A2-11. Nevertheless, this graph suggests that we need to choose frequencies where we compute PSD carefully if knowing magnitude of PSD of signals is important when signals are similar to S1, a constant amplitude, constant frequency, pure sinusoid without any phase variations. If amplitude of S1 varies, theoretical spectral will have some frequency width as we described in (Case 3) of (2-1-3-3).
The horizontal scale of "hills" and "valleys" shown by "x" marks in Figure A2-11 is about the same as horizontal distance between adjacent vertical red lines. Thus, the frequency interval of PSD computations we chose for Figure A2-11 is apparently much smaller than the actual frequency resolution of PSD. Frequency resolution of PSD computed by MEM is determined by AR order and there is a limit determined by a sampling interval. Computing PSD with a finer frequency interval as we did may help obtaining peak frequencies accurately but it would not give us PSD of a finer frequency resolution. Another thing we would like to mention is that, while MEM does not leak spectrum, we still can see the effect of S1 (red ellipse) away from the "mountain" of S1.
Since MEM does not assume that data are composed of pure sinusoids unlike PD001A/B (please, see equation (3)) results of MEM occasionally give surprisingly less accurate information when data are consist of pure sinusoids as in the case shown in (2-1-3-3). We generally recommend computing PSDs by several different methods and comparing results for serious spectral analyses.
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)