Introduction to Wavelet Analysis using SAS/IML: Time-Frequency Decompositions and Wavelet Scalograms
The purpose of this blog is to learn how wavelets can be used to analyze signal and sensor data. We learn what wavelets are and how they can be used to decompose periodic signals into their constituent frequencies. In this blog we focus on the basics of wavelet analysis and compare the time-frequency decomposition generated by wavelet analysis (we will explore a plot of this decomposition called the ‘scalogram’) with the spectrogram generated by traditional Fourier analysis (using the short-time Fourier transform). This is just one application of wavelets, with other popular applications including denoising and data compression. We will explore other wavelet applications in subsequent blog posts.
Wavelets are short oscillations with a finite duration in time. There are a variety of families of functions that serve as wavelets, with the common theme being a function that starts at zero and ends at zero with some oscillations in between. This is in contrast with the basis functions used in Fourier analysis, the sine and cosine waves that have infinite extent. Thus, from the perspective of Fourier analysis, wavelets are already “windowed” in that they are zero valued outside of a small, finite duration. Plotted below are two examples of simple wavelets (more complex wavelets are often used in practice). The first plotted wavelet is a Morlet wavelet and the second is a Haar wavelet.
Morlet Wavelet
Select any image to see a larger version.
Mobile users: To view the images, select the "Full" version at the bottom of the page.
Haar Wavelet
These “mother wavelets” are used to define a collection of similar wavelets at different scales and locations in time. From a given mother wavelet, we derive wavelets at different shifts and scales. Correlations between these wavelets and our signal of interest identify the signal components at different frequencies (scales) and at different points in time (shifts). This idea of ‘scaling’ and ‘shifting’ the wavelet is key to calculating the continuous wavelet transform (although for digital signals we will have to calculate the discrete wavelet transform), which is the time-frequency decomposition of our signal. Given a mother wavelet, ψ, we can define the child wavelet at scale a and shift b as follows:
We can think of the scale a as representing frequency; as a increases, the frequency of the oscillation in the wavelet increases. We can think of the shift b as representing time; as b changes, we move the wavelet forwards or backwards in time relative to the original signal. Of course, this all applies to continuous wavelets. In the discrete case (anytime we are dealing with digital signals we must use the discrete wavelet transform), we define the child wavelets based on integers m and n as follows:
In this case the integer m is associated with the scale and the integer n is associated with the shift. The most common implementation of the discrete wavelet transform involves shifting and scaling by powers of two, enabling much more efficient algorithms for computing the transform. We thus can simplify the equation above:
In this case the integer j is the scale parameter, and the integer k is the shift parameter. The set of child wavelets for integer values of j and k establish a complete orthonormal basis for square-integrable functions of a real variable. In other words, they form a basis for the kinds of functions we encounter when working with signal data. This is analogous to how sines and cosines form a complete basis when working with the Fourier transform. The basis used in wavelet analysis can do a better job balancing time and frequency resolution when decomposing our signal as compared to the Fourier transform basis. We’ll talk about why this is in more detail below.
When we perform the short-time Fourier transform, window functions are applied to break the input signal into a bunch of different, finite duration portions. The Fourier transform is then applied to each portion to calculate the frequency spectrum. Since each portion represents a different time point in the original signal, we can assemble the frequency spectra into a two-dimensional spectrogram where we plot time on the x-axis (corresponding to the different windowed portions of the signal) and we plot frequency on the y-axis (the information that was contained in each frequency spectrum). The time resolution we plot on the x-axis are determined by the location and number of time windows we apply to the input signal while the frequency resolution is determined by the length of signal in the windows (the size of the windows). Although we can seemingly improve frequency resolution in the Fourier transform by zero-padding, this just interpolates values in the frequency spectrum; the fundamental frequency resolution of the transform is determined by the length of the signal. Thus, there is a fundamental tradeoff in time-frequency resolution in the short-time Fourier transform. Longer windows mean more frequency resolution, but with a fixed-length signal longer windows means we have fewer windows overall, and thus the time resolution is limited. Shorter windows allow us to improve the time resolution by including more windows to cover the signal, but of course degrade the frequency resolution. This fundamental limitation is the same basic idea as the Heisenberg uncertainty principle in quantum mechanics (we cannot simultaneously identify the exact location and momentum of a particle) and is often called the Gabor limit in signal processing. Wavelet analysis avoids some of the consequences of this limitation by enabling us to vary the time-frequency resolution at different locations in the spectrogram (which we will now call a scalogram).
To calculate the wavelet scalogram for a digital signal, we find wavelet coefficients for the signal associated with each of the shifted and scaled wavelets we derived previously from the mother wavelet. These coefficients are analogous to the power in the time-frequency bins from the STFT and can be calculated by finding the dot-product between the input signal and the corresponding child wavelet (at a specific value of the shift and scale parameters). In this way we don’t have to window the input signal at all, since each of the child wavelets are already finite in duration (as opposed to the sines and cosines used in the STFT). The magnitude of these wavelet coefficients are plotted in a two-dimensional scalogram with the shift parameter representing the different times and the scale parameter representing different frequencies. The key advantage over the STFT results from the construction of the child wavelets from the mother wavelet. We have longer wavelets at lower frequencies (thus we have bad time resolution but good frequency resolution) and we have shorter wavelets at higher frequencies (thus we have better time resolution but worse frequency resolution). This allows us to more easily tell when unwanted high-frequency noise is included in a signal, while still allowing us to accurately identify long-lasting low-frequency oscillations in the signal. A toy visualization of the time-frequency resolution tradeoff in the STFT versus the wavelet transform is displayed in the following figures:
In the STFT the time-frequency resolution is fixed across all frequency scales.
In the wavelet transform we can have longer wavelets at lower frequencies and shorter wavelets at higher frequencies, allowing the time-frequency resolution to vary across frequency scales.
Let’s explore some examples of plotting time-frequency decompositions in SAS/IML. We will start by simulating a digital signal with a linearly increasing frequency and then plot the Fourier transform spectrogram and the wavelet transform scalogram.
proc iml;
/*start by simulating the chirped signal*/
f_s = 1000; /*sampling frequency*/
f_s2 = f_s/2; /*Nyquist frequency*/
t_s = 1/f_s; /*time between samples*/
N_samp = 2000; /*number of samples*/
pi = constant("pi");
t = do(0, (N_samp-1)*t_s, t_s); /*array of times*/
/*create the chirp*/
f0 = 10; /*starting frequency*/
c = 50; /*chirp, i.e. rate of change of the frequency*/
f_t = f0 + c*t; /*here is the linear chirp*/
/*f_t is the instantaneous frequency*/
phi_t = 2*pi*(f0*t + c/2*t#t); /*we integrate the linear chirp w.r.t. time and multiply by 2*pi to get the phase*/
/*phi_t is the instantaneous phase*/
x_in = sin(phi_t);
create work.input_signal var {"t","x_in","f_t"};
append;
close work.input_signal;
/*plot the instanaeous frequency and our signal in the time domain*/
submit;
title "Instantaneous Signal Frequency";
proc sgplot data=work.input_signal;
series x=t y=f_t;
xaxis grid label="Time (s)";
yaxis grid label="Frequency";
run;
title "Input Signal with Linear Chirp";
proc sgplot data=work.input_signal;
series x=t y=x_in;
xaxis grid label="Time (s)";
yaxis grid label="Amplitude";
run;
endsubmit;
/*now we do a short time fourier transform on the data to explore how the frequency changes*/
winLen = 500; /*this will have a big impact on the time/frequency resolution*/
window = tfwindow(winLen, "Gaussian");
noverlap = floor(0.6*winLen);
nfft = 512; /*increasing this increases the number of frequency bins*/
ramp = palette("Spectral",7)[7:1]; /* reverse "spectral" color ramp */
title "Short-Time Fourier Transform Spectrogram";
call spectrogram(x_in,window,noverlap,nfft) scale="linear"
colorramp=ramp sampRate=f_s;
/*next we do a wavelet transform and plot a scalogram
/*the wavelet module expects a column vector*/
x_in = x_in`;
/*initialize the wavelet graphics module*/
%wavinit;
/*set wavelet options using macro variables from %wavinit macro*/
/*in this case let's use a predefined Daubechies 3-tap wavelet*/
optn = &daubechies3;
/*now perform the discrete wavelet transform to find the decomposition*/
call wavft(decomp, x_in, optn);
/*plot a scalogram of the wavelet decomposition*/
title "Wavelet Transform Scalogram";
call scalogram(decomp);
quit;
This yields the following spectrogram and scalogram plots:
In this example the spectrogram does a better job of revealing the simple linear trend in the frequency, but we can still see the linearly increasing pattern the wavelet coefficients. Notice how the highest frequency detail coefficients in the wavelet scalogram have a high time resolution but are mostly empty of energy. This simulated signal is nice for illustrating how spectral decompositions works but is not representative of realistic scenarios. Most real signals involve high-frequency components riding on top of low frequency signals, so let’s explore another simulated signal with a temporally limited burst of high-frequency “noise” added to a low-frequency signal. The code for simulating the signal and plotting the spectrogram and scalogram is as follows:
proc iml;
/*We start by simulating a digital signal with a short burst of high-frequency noise*/
f_s = 500; /*sampling frequency*/
f_s2 = f_s/2; /*Nyquist frequency*/
t_s = 1/f_s; /*sampling period (time between samples)*/
N_samp = 1000; /*number of samples*/
t1 = do(0, (N_samp/2-1)*t_s, t_s);
t2 = do((N_samp/2)*t_s, (N_samp/2+5)*t_s, t_s);
t3 = do((N_samp/2+6)*t_s, (N_samp-1)*t_s,t_s);
f1 = 10;
f2 = 200;
pi = constant("pi");
x1 = sin(2*pi*f1*t1); /*signal for the first half*/
x2 = 0.1*sin(2*pi*f2*t2) + sin(2*pi*f1*t2); /*short piece of HF noise*/
x3 = sin(2*pi*f1*t3); /*signal for second half*/
t = t1||t2||t3; /*combined time from N=0 to N=200*/
x_in = x1||x2||x3; /*combined signal with short HF noise*/
create work.input_signal var {"t","x_in"};
append;
close work.input_signal;
/*plot our signal in the time domain*/
submit;
title "Input Digital Signal";
proc sgplot data=work.input_signal;
series x=t y=x_in;
xaxis grid label="Time (s)";
yaxis grid label="Amplitude";
run;
endsubmit;
/*Now we obtain our signal's frequency domain representation using STFT*/
winLen = 50; /*this will have a big impact on the time/frequency resolution*/
window = tfwindow(winLen, "Gaussian");
noverlap = 0;
nfft = 64; /*increasing this increases the number of frequency bins*/
ramp = palette("Spectral",7)[7:1]; /* reverse "spectral" color ramp */
title "STFT Spectrogram";
call spectrogram(x_in,window,noverlap,nfft) scale="log"
colorramp=ramp sampRate=f_s panel=1;
x_in = x_in`;
/*initialize the wavelet graphics module*/
%wavinit;
/*set wavelet options using macro variables from %wavinit macro*/
/*in this case let's use a predefined Daubechies 3-tap wavelet*/
optn = &daubechies3;
/*now perform the discrete wavelet transform to find the decomposition*/
call wavft(decomp, x_in, optn);
/*plot a scalogram of the wavelet decomposition*/
title "Wavelet Transform Scalogram";
call scalogram(decomp);
quit;
This yields the following spectrogram and scalogram:
We can see in the spectrogram that there is some amount of high-frequency noise added to the signal at about 1 second. The noise is low power, so it’s exact frequency is ambiguous in the spectrogram. The width of the time window also makes it hard to pinpoint when this noise was added to the signal.
In the wavelet scalogram we still don’t have good frequency resolution at the higher frequencies, but we can see that there was some high-frequency noise added about halfway through the signal. At the 9th level we can see a very narrow yellow line, corresponding to the exact portion of the signal where the high-frequency noise was added. The wavelet scalogram does a much better job of localizing the high-frequency noise in time than the STFT spectrogram.
Wavelets are a powerful tool for analyzing signal and sensor data. We have seen how we can use wavelets to decompose periodic signals into time-frequency representations that often have more useful resolution properties than the STFT spectrogram, especially for real-world signals. The examples in this blog just scratch the surface of what we can do with wavelets, and we will explore applications of wavelets such as denoising and data compression (for one-dimensional signals and two-dimensional images) in future blog posts.
References:
- Previous article on time-frequency analysis using Short-Time Fourier Transform
- SAS Documentation describing wavelet analysis
- SAS Documentation page listing the wavelet options and symbolic names
Don't miss out on SAS Innovate - Register now for the FREE Livestream!
Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.
