|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
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 equal to the actual signal amplitude (1.0) as shown in Figure 4-32(b). Thus, simple 3-bin average is 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 signal bin matches actual signal frequency and there is no spectral leakage. Figure 4-32(c) shows time series of data when we applied Hanning window. If we applied Hanning window, width of signal 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(0) and Po(5).
(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.
What we would like to say here is that it is reasonable to think that the values of purple and red lines outside of the signal peak in Figure 2-5(a) are zero although actual numbers are not zero.
Appendix 2; Spectral leakage of power spectrum density function and window functions.
In (4-2-1) we described that the 'cliff' at boundaries 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 leakage is generated and how large they would be. Finally we describe how to reduce spectral leakage. The root cause of spectral leakage 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) to "Our data" x'(n). All the values of z(n) are zero. Since Fourier transform of z(n) is also zero for any frequency, adding z(n) to "Our data" x'(n) would not modify spectral characteristic of x'(n) at all. By adding z(n) to x'(n), we will get "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 in order 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.
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 deviation of spectral of "Our data 2" shown in Figure A2-1(c) (and Figure A2-2(c)) from spectral of infinitely long data shown in Figure A2-1(a). The latter of which is consdered to be a "true" spectral. This approach makes sense because, although "Our data", limited length data shown in Figure A2-2(a), is what we usually have, spectral characteristic of it is the same as the spectral characteristic of "Our data 2". 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 would 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 would affect spectral of resultant windowed data in this appendix.
We cannot evaluate power spectrum density of "Our data 2" at discrete frequencies anymore because data length of "Our data 2" is infinite and frequency bandwidth of bins becomes zero (=1/infinite). Instead of using PSD we introduce similar quantity called PS for infinitely long data as follows. First we extend the range of summation in equation (18) as
For 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 at negative frequencies. Another apparent difference is that equation (A2-3) is not divided by T but this is simply because T is infinite for 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 Figure A2-1(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 leakage is generated and how large it would be
Equation (A2-5) shows that the Fourier transform of x(t) is 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 caused by 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 have used. Each data is shown as a small dot in this graph and we connected dots by straight blue lines.
The first case (Window-width A) is when window-width is 40 (data points). This is the case when frequency match is achieved and this case is equivalent to the data length A shown in Figure 2-4. 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 is actually the start point of the next cycle and does not belong to the second cycle. Another case (Window-width B) is when window-width is 50 and this case is similar to the case of 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 A in equation (A2-4) (amplitude of x(t)) is 1.0 and 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 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 no accident that the tip of this arrow touches PS (red line). The values of modified PSD (blue arrow) and PS (red line) at the same frequency must be the same 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 frequencies 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 PSDs of all the bins are non-zero due to the spectral leakage, 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, values of modified PSD and PS at the same frequency are the same (tips of arrows always touch red line). Figures A2-3(b) and (c) suggest 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 labeled fT 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(t) and the lower panel of Figure A2-4 shows the Fourier transform of a Boxcar window function.
These are graphical representations of equations (A2-5) and (A2-7), and the fundamental frequency ranges (horizontal width of 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, the first thing we do is flipping lower graph left and right. This is already done in this graph so that frequency (horizontal axis) of the lower graph decreases from left to right unlike upper graph (red circles). 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 parts in these graphs.
Then, we draw lines vertically from the upper graph to the lower graph at frequencies f0 and -f0 (vertical dashed lines H and L). Since the frequency difference between H and C (rightmost vertical dashed line) is fT-(-f0)=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 are 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 will 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 signal bin frequency becomes zero and we do not see spectral leakage at those frequencies. However, if we compute PS at any other frequency, we will pick up 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.
Equation (A2-8) shows magnitude of spectral leakage at any frequency is determined by the Fourier transform of a window function. Lower part of Figure A2-4 shows that the value of the Fourier transform of a Boxcar window function oscillates and the amplitude of this oscillation decreases as frequency deviates away from zero in either direction. This results that the magnitude of spectral leakage tends to become smaller as the difference between fT and f0 becomes larger. However, spectral leakage will never be zero in this case.
Figure A2-7 shows convolution process for general cases when the Fourier transform of infinitely long data (X in Figure A2-7(a)) 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, frequency width (horizontal width in this graph) of which is the same as frequency width 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) as before. We shift horizontal position of the graph of the Fourier transform of a Boxcar window function (W in Figure A2-7(b)) after flipping it left and right as before. Here, the frequency width of W is the same as the frequency width of segment A.
For the convolution we chose segment E of X. The frequency width of segment E is the same as frequency width of segment A. 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 (sum) 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 similar 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 at frequencies f that satisfy fT - f = +/- k/T (second line of equation (A2-7)).
(3) How to reduce spectral leakage
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 over 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 logarithmic scaling for this graph unlike others. From this graph we can see that Hanning window would considerably reduce spectral leakage and reduction rate would progressively increase as the difference between fT and f0 (Figure A2-4) becomes larger. This graph also shows that, although applying Hanning window removes 'cliff' described in (4-2-1), we will still have spectral leakage.
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 are the cause of widening of signals 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 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 leakage (Not supported in this package deal service)
We will start showing that the overall general characteristic of PSD of data, that contains various signals and noises such as SR, computed by MEM is quite similar to the overall general characteristic of PSD computed by PD001A/B. Figures A2-9(a) and (b) show PSDs of SR computed by PD001A/B (blue and red lines) and by MEM (purple lines). We moved blue and red lines upward to make comparison easier. The data, SR, we used here is 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 shown 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 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 purple 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 would have similar practical frequency resolutions to those of others in Figure A2-9. We can compute PSD at arbitrary frequencies if we use MEM but we computed PSD at the frequencies where we would compute PSD if we use PD001A/B (blue and red lines) for the sake of comparison.
The point we would like to make here is that figures 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 of SR alone. Figure A2-10(a) shows PSDs computed by PD001A/B when we did not apply any window (equivalent to applying Boxcar window) and Figure A2-10(b) shows PSDs computed by PD001A/B when we applied 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 purple 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 Hanning window but existence of it is still 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 and constant frequency pure sinusoid without any phase variations. To make this point clearer we recomputed PSD with a much smaller frequency interval. The blue (SR+S1) 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. Ideally, 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 steep "mountain") as shown in Figure A2-11. Nevertheless, this graph suggests that we need to choose frequencies where we compute PSD carefully if it is important to know magnitudes of PSD of signals when signals are similar to S1, a constant amplitude and constant frequency pure sinusoid without any phase variations. If the magnitude of S1 is comparable to those of S2 and S3, we might even miss S1. However, if amplitude of S1 varies, theoretical peak of PSD due to S1 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 sampling interval and number of data. Computing PSD with smaller frequency interval as we did may help obtaining peak frequencies accurately but it does not necessarily 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 is composed of pure sinusoids unlike PD001A/B (please, see equation (3)) MEM occasionally gives us unexpected results when data is 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 data is short relative to the period of S1, f1 will be included within zero-frequency bin because frequency bandwidth of bins (frequency resolution) is so large. 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 very, very slowly varying. Please note that this data length is short relative to the period of S1 but it is fairly long relative to the periods of SL and SH as shown in Figure 2-11.
As we increase data length, frequency resolution of the Fourier transform becomes finer and f1 leaves from zero-frequency bin but moves toward higher frequency bin. Then, convolution of the Fourier transform of S1 and the Fourier transform of S0 starts generating twin-peaks center of which is at f0.
The alternative way of understanding why we need more than certain length of data to have twin-peaks is considering this as a problem to separate two neighboring singnals, frequency difference between them is 2 x f1.
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)