**NON-STATIONARY SPECTRA :**
This tutorial covers the spectral analysis capabilities of AutoSignal for non-stationary data. This includes the Short-Time Fourier Transform (STFT) and the Continuous Wavelet Transform (CWT) spectral procedures. These procedures are useful for studying the changing properties within a data stream as well as identifying stationary regions for conventional Fourier and parametric analysis.
**Generating A Test Signal Containing Four Sequential Sinusoids**
Select the Generate Signal option in the Edit menu or Main toolbar.
For this tutorial, we will create a non-stationary data set with four data regions each containing a sinusoid of different frequency.
Click Read and select the file tutor5.sig from the Signals subdirectory.
The following signal expression is imported:
SRATE=5000
NYQ=SRATE/2
AMP1=1
AMP2=1
AMP3=1
AMP4=1
FREQ1=NYQ*0.1
FREQ2=NYQ*0.2
FREQ3=NYQ*0.3
FREQ4=NYQ*0.4
PHASE1=PI/2
PHASE2=PI/4
PHASE3=PI
PHASE4=3*PI/2
F1=IF(X<=0.05115,AMP1*SIN(2*PI*X*FREQ1+PHASE1),0)
F2=IF(X>0.05115 && X<=0.1023,AMP2*SIN(2*PI*X*FREQ2+PHASE2),0)
F3=IF(X>0.1023 && X<=0.15345,AMP3*SIN(2*PI*X*FREQ3+PHASE3),0)
F4=IF(X>0.15345,AMP4*SIN(2*PI*X*FREQ4+PHASE4),0)
Y=F1+F2+F3+F4
The X (time) values vary from 0 to 0.2047 with a 0.0002 sample increment. The Nyquist frequency is thus 2500 (half the 5000 sampling rate). The first sinusoid in the data stream is at frequency 250, the second is at frequency 500, the third is at frequency 750, and the last sinusoid to appear is at frequency 1000. The four components are limited to the lower half of the Nyquist band. The amplitudes are set to a constant to illustrate an important aspect of CWT analysis. 10% random Gaussian (white) noise is added.
Click OK to process the current signal.
An AutoSignal graph is presented containing the 1024 point generated data. The stepwise increase in frequency is apparent.
Click OK to accept the generated data. Click Yes when asked to update the main data table with the revised data.
The Loss Of Time Information In A Fourier Transform
We will begin by exploring a basic Fourier spectrum of this data set.
Select the Fourier Spectrum option in the Spectral menu or toolbar. Use the Best Exact N algorithm with Nmin set to the 1024 data length. Set the plot to dB and the signal count (sig) to 4.
The Fourier spectrum readily confirms the frequency of each of the four components. All time-localization information is lost. Each component is assumed to be infinite in duration.
**Isolating Components**
We will not make the claim that wavelet or STFT analysis is the only way to analyze non-stationary signals. Indeed, it is often simpler to isolate components by frequency and process the reconstructed data.
Select the Apply Fourier Decomposition and Reconstruction to Local Copy of Data option. Be sure that the window is set to None, that the Best Exact N algorithm is selected with Nmin set to the 1024 data length, and the plot set to dB. Using the left mouse, box the third spectral peak including the first couple of adjacent harmonics.
This approach clearly reveals the range of time where this signal is present. Although this is a very simple illustration, the same approach can be applied to far more complex data analysis.
Click OK to accept the local filtering. Change sig to 1. Note the true brickwall filtration.
Select the Section Local Copy of Data option.
Enable the XY Sectioning mode and box the region shown.
Click OK to accept the sectioning.
Click the Numeric Summary button and inspect the Sine Component Fit section.
**Sine Component Fit (Suboptimal)**
**Frequency** |
**Amplitude** |
** Phase ** |
**Power** |
**%** |
**Rel %** |
749.939176 |
1.01778903 |
3.19475563 |
0.02052905 |
100.000000 |
100.000000 |
Given the 10% random noise, it is safe to say that the 750 frequency, the 1.0 amplitude, and the 3.14159 phase were accurately recovered. Due to the random noise component and variable sectioning, the results you see will vary somewhat from those shown here.
**Close the Numeric Summary and exit the Fourier Spectrum procedure.**
Because the Fourier filtration and sectioning were local options to the spectral procedure, the original data are restored. To permanently alter the data, options such as these must be invoked from the main program menu or toolbars.
Although it is possible to isolate signal components by frequency range or signal strength, and to then reconstruct them for analysis using spectral procedures that require stationary data, it is often convenient to analyze a time-frequency representation of the data.
**Short-Time Fourier Transform**
The Short-Time Fourier Transform (STFT) is similar to the periodogram or segmented-overlapped FFT, except that a central time value is assigned to each segment and unless there is a need to directly retrieve the amplitude or power from the spectral peaks, a data tapering window is usually used to create the best possible time localization.
Time-frequency spectra are intrinsically fuzzy. The better the resolution in time, the poorer is the frequency resolution. The sharper the resolution in frequency, the poorer is the time resolution. The STFT is the best way to get a good feel for this effect.
Select the Short-Time Fourier Spectrum option from the Spectral menu or toolbar. Set the window to cs2 Hann, the segment length (Seg n) to 256, the overlap % to 90, the FFT length (FFT n) to 256, the plot to Int=PSD TISA, and the 3D Graph Profile to Gradient.
Optimizing the STFT parameters can be challenging. The above settings generate 31 separate segments. The frequency localization is very good, due to the 256 data length segments that are transformed. On the other hand, these 256 data points span 25% of the time range. The time localization is poor. A contour plot can help in visualizing this.
Set the 3D Graph Profile to Contour.
The sharp definition in frequency is apparent as is the extremely fuzzy time localization. If the time localization were well defined, there would be no overlap in time between the components. Note also that half of the first and fourth components are not mapped since the first and last central time values are 128 data points in from the edges of the data.
Change the overlap % to 50, Seg n to 64, and FFT n to 64.
Note that the 64 data length segments with a 50% overlap also result in the same 31 segment count. Some frequency localization has been lost due to the reduced length transforms. The time localization, however, is markedly improved. The edge effect zones have also been diminished by a factor of four. Most of the first and fourth components are now mapped.
The STFT is a single resolution technique. The same time-frequency tradeoff is present everywhere in the spectrum regardless of frequency. Although the STFT does require careful fine tuning, an optimized spectrum can readily yield spectral amplitudes and powers directly from the height of spectral peaks.
Set the window to None, Seg n to 70, overlap% to 70, FFT n to 256, plot to Amplitude, and 3D Graph Profile to Shaded.
Click the 3D View button. Adjust the Z View Angle to 0 degrees using the scrollbar to the right of the plot. Click OK.
Note that the amplitude peak for all four spectral components is very close to the true 1.0 value.
AutoSignal's 3D graphs offer four profiles that initially produce commonly used plots. These profiles are fully customizable and are automatically saved across sessions. You can thus save any four 3D formats of value to you in assessing time-frequency spectra.
In the STFT, a separate FFT is generated for each segment. The spectral surface is saved and rendered using a B-spline interpolant. If the total rendering grid would be greater than 16384 elements, a decimated spectrum is used instead. All of these elements consume significant memory. AutoSignal limits the number of STFT segments to a maximum that will be in the vicinity of 512. Zero padding increases the frequency density and can drastically increase memory requirements. The various STFT spectra in this tutorial should have been rendered very swiftly. If you experienced excessive disk activity and extremely slow processing and procedure closure times, the available physical memory is insufficient for efficient computation and rendering of STFT spectra.
**Continuous Wavelet Transform Spectra**
The STFT uses phase-bearing sine basis functions. The infinite duration of sine functions requires the transform of windowed segments of the data stream. If the basis functions were themselves of finite duration, a spectrum could be created from the correlation between shifted and stretched/compressed versions of this limited duration mother "wavelet" and the data stream. This is the essence of the Continuous Wavelet Transform (CWT).
Despite the awesome array of wavelet shapes, using wavelets for time-frequency spectral analysis is generally easier than optimizing the STFT. And in most instances, the differences between the various mother wavelets and settings do not greatly impact the resulting time-frequency spectrum. The mother wavelet and adjustable parameter selection tytut5_imgally involve selecting a time-frequency tradeoff.
Select the Continuous Wavelet Transform Spectrum (3D Surface) option in the Spectral menu or toolbar. Be sure the wavelet is the Morlet, that the Complex box is checked, that the wavelet's adjustable parameter (Adj) is set to 8, that Nmin is set to 2048, that Full Range is checked and Ln Steps is not checked, that the number of frequencies (n) is set to 60, that the Int=PSD TISA plot option is selected, and the 3D Graph Profile is set to Gradient.
This is a very good CWT time-frequency spectrum. In fact, respectable results are attained regardless of the mother wavelet or the value of its adjustable parameter. You could likely use the Morlet wavelet with a wavenumber of 8 for all time-frequency analysis work, and never miss the small refinements that would derive from selecting the wavelet whose shape is most compatible with the oscillations in the data and whose oscillation count yields the optimum time-frequency resolution. To see this fact, we will explore the extremes offered in AutoSignal.
Set the Adj value to 20. Change the 3D Graph Profile to Contour.
The Morlet with an adjustable parameter (wavenumber) of 20 produces the best frequency localization and the worst time localization. Note, however, that this worst case time fuzziness is much less than the previous ill-designed STFT example.
Set the Wavelet to GaussDeriv and the Adj value to 2.
The Gaussian Derivative wavelet with an adjustable parameter of 2 produces the sharpest time localization but as is apparent, very fuzzy frequency localization.
In all cases, though, the power of the four spectral features is conserved. The difference rests completely in the extent to which the 3D spectral peak is smeared in time and the extent to which it is smeared in frequency.
**Wavelet Shapes**
Click the View Mother Wavelet button.
Wavelets are usually complex time domain functions. The white curve in the lower graph is the real component and the blue is the complex component. This shape produces a very compact single positive frequency peak in the frequency domain.
Select the Morlet wavelet and set the Adj value to 8.
This is AutoSignal's default for CWT analysis. This option is used for both graphical wavelet selection and wavelet visualization.
Click OK to update the spectrum with this selection.
**Multiresolution Analysis**
The CWT is a multiresolution technique. If you inspect the previous plots, you will note that the lowest frequency peak has the highest CWT spectral magnitude and is least smeared in frequency of the four signal components. The highest frequency peak consistently evidences the lowest spectral magnitude and is most smeared in frequency of the four components.
This property is a byproduct of the CWT algorithm rather than a design consideration. Multiresolution analysis has its benefits and drawbacks. If high frequency components are intermittent with a very short time duration, the optimum position would be to tradeoff frequency resolution in order to accurately capture its appearance and disappearance in time. If low frequency components are long-lived, the optimum state would be to tradeoff time resolution for more accurate frequency estimation. When these conditions are true, the multiresolution analysis in the CWT is an asset.
If persistence is independent of frequency, or if the primary aim is to retrieve amplitudes and power directly from the spectral magnitudes, the CWT's multiresolution property is a drawback.
**Wavelet Peak Significance**
Because of the multiresolution property, critical limits based on spectral peak magnitudes are more complicated. The CWT magnitude is a property of the smearing in both the time and frequency dimensions, and this varies with frequency. Further, the smearing in time and frequency is a product of the wavelet and its adjustable parameter.
AutoSignal offers peak-based critical limits based upon the wavelets, adjustable parameters, and algorithms in the product. Monte-Carlo data were converted into bivariate models specific to each wavelet, taking into account the wavelet's adjustment as well as the data size being analyzed. The critical limits are thus accurate and valid, although their appearance is a bit different from those used in the various frequency spectra.
Change the 3D Graph Profile to Gradient. Check the AR(1) Bkgrnd box.
The critical limits for the time-frequency spectra are plotted as color gradients. These override any Z-gradient coloring. The wavelet critical limit gradients are the following colors by default: 8-level grayscale from 10 to 50%, 8-level cyanscale from 50% to 90%, 8-level greenscale from 90% to 95%, 8-level yellowscale from 95% to 99%, and 8-level redscale from 99% to 99.9%. The 99% critical limit thus exists at the transition boundary between yellow and red. A 99% critical limit means that in only 1 of 100 white noise data sets of equivalent length and variance did a CWT peak attain this magnitude due to chance.
Keep in mind that critical limits apply only to the largest magnitude peak in the spectrum and their purpose is to address whether or not a spectrum can be considered statistically distinct from white noise. Although AutoSignal renders the critical limit gradient throughout the spectrum, its interpretation should rest with only the dominant peak. Here the largest magnitude peak is very close to the 99.9% upper limit. In only 1 of 1000 random noise sets would a peak with this wavelet power be encountered from chance.
The twist in the interpretation is that the peaks in white noise CWT spectra are not uniformly distributed across all frequencies. Due to the multiresolution property, the largest peak in a given spectrum is often more likely to occur at a lower rather than higher frequency. In a basic Fourier spectrum, four equal magnitude spectral peaks will evidence equivalent power and equivalent positions relative to the critical limits. In wavelet spectra, this is not true.
In a Fourier spectrum, a secondary peak at a 95% critical limit can be assumed significant to at least this level, understanding that its significance may be appreciably higher if the primary component were mathematically removed. The same is true for wavelet spectra except that the secondary spectral magnitude is likely to be reduced even further by this multiresolution effect.
Click OK to close the wavelet procedure.
**Continuous Wavelet Transform Spectra as 2D Contour**
If you prefer to view time-frequency spectra mainly using 2D contour plots, AutoSignal offers a separate option that uses less memory and renders the contour more swiftly and with a greater host of display options.
Select the Continuous Wavelet Transform Spectrum (2D Contour) option in the Spectral menu or toolbar. Be sure the wavelet is the Morlet, that the Complex box is checked, that the wavelet's adjustable parameter (Adj) is set to 8, that Nmin is set to 2048, that Full Range is checked and Ln Steps is not checked, that the number of frequencies (n) is set to 60, that the dB plot option is selected, and the dBlim value is set to 24.
Select the Modify Contour Properties option. Be certain the evaluation grid is set at 60 and that the 48 Spectrum contour is selected. Click OK.
When the spectrum is plotted with a dB contour, each color in the gradient corresponds to a specified dB delta. In this case, there are 48 colors, each covering a range of 1/2 dB.
If the YR (Y-right) labels and titles are on in the 2D View settings, and provided the Y title is "Frequency", the right axis will automatically display the Period (wavelength).
**CWT Frequency Redundancy**
Unlike the Discrete Wavelet Transform (DWT) traditionally used for wavelet analysis, the CWT can be evaluated at any set of frequencies desired. These frequencies can overlap with an arbitrary measure of redundancy and they can be spaced linearly, logarithmically, or in any sequence. Because a separate FFT must be stored for each frequency evaluated, AutoSignal limits the total frequency count to 100. It is this frequency redundancy that allows the CWT to offer a far higher resolution than the DWT.
Uncheck the Full Range box and enter 150 for the starting frequency (f start) and 400 for the ending frequency (f end).
The full analytical power of the CWT spectral analysis is now concentrated in the 150 to 400 frequency band.
**FFT Zero Padding And Edge Effects**
Unlike a conventional Fourier analysis where the spectrum is generated by the transform, the FFT is used in the CWT is strictly for performing fast convolution of the scaled and translated wavelet with the data stream. In the CWT the basis functions are scaled and shifted wavelets, not the Fourier sines and cosines. Zero padding has only a single purpose in the CWT, that being to avoid wraparound effects in the fast convolution procedure.
Wraparound is fully avoided at all frequencies when the FFT length is twice the length of the data stream. It is often possible to zero pad to the next power of two and avoid any perceptible wraparound and at the same time achieve the fastest possible convolutions.
Check the Full Range box and change Nmin to 1024 using the drop down box.
When an exact N FFT is used for the fast convolution, the wraparound at low frequencies is readily apparent. This wraparound effect diminishes with frequency. This explains the absence of a wraparound effect with the highest frequency component. If there is negligible low frequency content, or if the data are cyclic in the band processed, there is no need to use zero padding.
When zero padding is used to avoid this effect, the power associated with the wraparound is discarded. This results in a diminished power near the edges of the data record. In CWT terminology, this zone where the power may be lessened is known as the cone of influence.
Check the Add Data box.
The data points within the spectrum that fall within this cone of influence are grayed (by plotting them using the current inactive point color). The cone of influence is particularly important when the frequency spacing is logarithmic.
**Surface Interpolation**
The wavelet spectrum is rendered using a bivariate B-spline interpolant. This interpolant is used for rendering both the 3D surface and 2D contour wavelet options. It is also used for the various power analysis options. Since the wavelet spectrum is represented by a continuous model, a couple of evaluation options are also available.
Click the Quick Evaluation option. Click the Surface Max box. Observe that the maximum in the interpolated surface is in the center of the time and frequency associated with the first component. Close the Quick Evaluation window.
Click the Evaluation option. Click Generate Table. Be sure the Data Source is Generate, Write Data is To Evaluation Table, and X and Y Input, Z=Fn at X,Y is checked. Click OK. Click OK to accept the X and Y ranges and increments. Note the generated data as well as the manual evaluation options such as roots, derivatives, and integrations available in this procedure. Click OK to close the evaluation option.
Click OK to close the 2D wavelet procedure.
**Power Analysis - Frequency Range**
Because CWT peaks do not offer a visual indication of power, it is necessary to integrate this interpolant in order to generate meaningful power information.
Select the Continuous Wavelet Spectrum Frequency Range option in the Spectral menu or toolbar. Be sure the wavelet is the Morlet, that the Complex box is checked, that the wavelet's adjustable parameter (Adj) is set to 8, that Nmin is set to 2048, that Full Range is checked and Ln Steps is not checked, that the number of frequencies (n) is set to 60, that the Int=PSD TISA plot option is selected, that ndt is set to 100, f init is set to 2.44141, and that f final is set to 2500.
The ndt of 100 specifies 100 different time snapshots of power. For each of these 100 time windows, the surface interpolant for the CWT spectrum is integrated between f init and f final. In this case, we are integrating the entire frequency range. Note that the power is relatively constant over time. Note the diminished power at the start of the data due to the cone of influence. Note also that the right hand power scale's initial value is well above zero. Only the upper 40% of the power range is plotted.
This option can be used to investigate the power present within a non-stationary signal across time within any frequency band desired.
Click OK to close this procedure.
**Power Analysis - Time Range **
A time range option is also available.
Select the Continuous Wavelet Spectrum Time Range option in the Spectral menu or toolbar. Be sure the wavelet is the Morlet, that the Complex box is checked, that the wavelet's adjustable parameter (Adj) is set to 8, that Nmin is set to 2048, that Full Range is checked and Ln Steps is not checked, that the number of frequencies (n) is set to 60, and that the Int=PSD TISA plot option is selected. Set ndf to 256, and be sure t init is set to 0 and t final is set to 0.2046.
This produces the "global wavelet spectrum". By integrating across all times, the time localization information is discarded just as is true with a basic FFT. The global wavelet spectrum offers neither high resolution nor high accuracy frequency estimation, and the multiresolution property makes it impossible to assess power from the heights of the peaks.
Although it may seem foolhardy to go to such great lengths to create a time-frequency representation of the data only to discard that time localization by a subsequent integration, this option does have merit. The global wavelet spectrum is a power spectrum based upon a set of wavelet basis functions rather than sinusoids.
The most important utility may rest with those instances where power is present at approximately the same frequency at different times. Having a tut5_imgture of time-frequency space, it is straightforward to generate several spectra for the different time bands in order to observe the subtle differences in frequencies and power of the various occurrences. To generate power spectra for different segments of time, you need only change the t init and t final values.
Click OK to close the procedure. |