NDTools
Tutorial for Matlab
Phillip M. Gilley & Nicholas K. Walker
Neurodynamics Laboratory
University of Colorado, Boulder
Overview
The CWT is a method for transforming 1-dimensional signals (e.g., time-amplitude
waveforms) into the time-frequency space. While several methods for computing the
CWT exist, we will be using a method first proposed by Torrence & Compos (1997),
which makes use of the convolution theorem.
In the following tutorial, we will examine the properties of the CWT while
demonstrating the functionalities of the NDTools Toolbox. In current form, the
toolbox requires two Matlab toolboxes also be installed: signal processing & image
processing. This overview contains XX Sections, each corresponding to sections in
the accompanying Matlab m-file, a script that performs all functions that generated
figures in this tutorial.
i
1
Generating a Signal
mkpsin( )
The examples in this tutorial will use a variety of signals and noises that are generated by
functions included in the NDTools toolbox. These functions can be used to create basic
sinusoidal signals, enveloped sinusoidal signals, and randomly generated noises. In this
first section, the function mkpsin() will be used to generate a frequency sweep, which will
be further analyzed in subsequent sections.
The function mkpsin( ) is used to generate sinusoidal signals. The function allows
specification of frequency, number of cycles to include, and options to pad the output to a
specified length, and to apply an amplitude envelope around the signal. In the first
example, mkpsin( ) will be used to generate 5 cycles of a 10 Hz wave. To compute
frequency, the sampling rate of the signal must be determined a priori. The following
Matlab code can be used to generate this signal.
% define the input variables:
F=10;
% frequency in Hz
nC=5;
% number of cycles
Fs=1000;
% sampling rate in Hz
% generate the signal:
signal = mkpsin(F,Fs,nC);
% plot the signal:
figure;
plot(signal);
axis tight;
xlabel('points');
ylabel('amplitude');
title('Signal: 10Hz, 5cycles');
2
The above code will generate the following figure:
The Matlab workspace now contains the following variables:
The signal returned by mkpsin() has a length of 501 points. The function can also return a
signal with a specified length, padding the signal with zeros on each side. In the following
example, the above signal will be generated and forced to have a length of 1000 points.
With a sampling frequency of 1000 Hz, the length of the signal will be 1 second, with one
point for each millisecond in the signal.
3
% define the input variables:
F=10;
% frequency in Hz
nC=5;
% number of cycles
Fs=1000;
% sampling rate in Hz
pnts=1000; % number of cycles
% generate the signal:
signal = mkpsin(F,Fs,nC,pnts);
% plot the signal:
figure;
plot(signal);
axis tight;
xlabel('points');
ylabel('amplitude');
title('Signal: 10Hz, 5cycles');
The result of specifying a signal length can be seen in the above figure. That feature can
be useful when generating multiple signals for comparison, allowing for control over
multiple signal features. In addition to a basic sine wave, mkpsin( ) can generate an
amplitude envelope over the sine wave. In the following example, the above signal will be
generated to include an amplitude envelope modeled by a Hamming window.
4
% define the input variables:
F=10;
% frequency in Hz
nC=5;
% number of cycles
Fs=1000;
% sampling rate in Hz
pnts=1000;
% length of the signal
env='hamming'; % window type for envelope
% generate the signal:
signal = mkpsin(F,Fs,nC,pnts,env);
% plot the signal:
figure;
plot(signal);
axis tight;
xlabel('points');
ylabel('amplitude');
title('Signal: 10Hz, 5cycles, Hamming window');
The result is an amplitude modulated sine wave, but with only one cycle of the modulation
envelope. Enveloped signals are useful for modeling bursts of activity, and transient
activity. The same signal can be returned without the padding, by setting the pad length
to -1 (the variable pnts in the above code).
5
To generate a frequency sweep, a signal for each frequency can be generated by
mkpsin( ) and concatenated to create one signal. In this example, a frequency sweep will
begin at 90 Hz and sweep downward to 10 Hz and back up to 90 Hz in 1 Hz steps, and
with 1 cycle at each frequency.
% define the input variables:
f=[90:-1:10 11:1:90];
% frequencies in Hz
nc=1;
% number of cycles
Fs=1000;
% sampling rate
% initiate signal:
signal=[];
% for each frequency, add a cycle to the signal:
for i=1:length(f)
y=mkpsin(f(i),Fs,1);
signal=[signal y];
end
clear i y
% plot the signal
figure;
hold on;
plot(signal);
axis tight;
ylim([-1.1 1.1]);
xlabel('points');
ylabel('amplitude');
title('Frequency Sweep: 90-10-90 Hz');
6
2
Frequency Spectrum
plotspec( )
A primary tool in signals analysis is the Fast Fourier Transform (FFT), a method to quickly
determine the frequency content of a signal. The NDTools Toolbox makes use of Matlab’s
native FFT function and its inverse: fft( ) and ifft( ). NDTools contains a function called
plotspec( ), which computes the FFT of a signal, and returns the one-sided frequency
spectrum (also known as the Fourier spectrum). The function plotspec( ) can perform
several tasks. First, if no output is specified when using plotspec( ), then a figure is
generated with a plot of the frequency spectrum. Additional options to limit the frequency
range of the plot can be specified in the input. The plotspec( ) function also has the
option to return the results of the transform, including the array of frequencies and the
spectral power at each frequency. An additional option to disable plotting is available, so
that plotspec( ) can be used to quickly compute the one-sided spectrum of a signal, the
results of which can be used in subsequent analyses.
In the following example, a spectrum plot of the frequency sweep from Chapter 1 will be
generated. Inputs to the plotspec( ) function include the signal to be analyzed and the
sampling rate. In addition, an optional input to limit the frequency range from 1 to 100 Hz
will be included.
7
%
Create a new figure:
figure;
hold on;
%
Plot the spectrum:
plotspec(signal,Fs,'on',[1 100]);
title('Spectrum of 90-10-90 Sweep');
For comparison, the following is the spectrum of the 10 Hz, amplitude modulated sine
wave from Chapter 1:
8
To retrieve the results of the analysis used by plotspec( ), outputs can be specified, and
the plotting function turned off. The following code will return two variables to the
workspace.
% get the Fourier spectrum:
[magnitudes, frequencies] = plotspec(signal,Fs,'off');
The workspace now contains two new variables:
The plots that are returned by plotspec( ) are created by plotting the magnitudes as a
function of frequency. In standard Matlab code, the same plot can be created as follows:
% create new figure
figure;
hold on;
% plot spectrum
plot(frequencies,magnitudes);
axis tight;
xlabel('Frequency (Hz)');
ylabel('Relative Magnitude');
title('Spectrum of 10 Hz AM Signal');
% limit frequency range from 1 to 100 Hz:
f1=abs(frequencies-1);
f2=abs(frequencies-100);
f1=find(f1==min(f1),1);
f2=find(f2==min(f2),1);
xlim([f1 f2]);
9
3
Continuous Wavelet Transform
mkwobj( ) qwave( )
A limitation of the FFT, demonstrated in the previous Chapter, is that it is not possible to
simultaneously examine time and frequency. That is, we can only observe the signal’s
content over time, or over frequency, but not both. However, in many applications it is
desirable to examine changes in frequency over time. One approach to resolving this
time-frequency tradeoff is to compute the spectrogram of a signal, essentially a moving
averaging of many FFTs. However, such an approach has limitations with resolution in the
time and frequency space. Wavelet analysis is a popular method for time-frequency
decomposition of a signal that has a superior time-frequency resolution. While several
methods are available to examine time-frequency, the NDTools Toolbox employs a method
for the continuous wavelet transform (CWT) described by Torrence & Compos (1997).
The CWT is an extension of the more classical wavelet decomposition known as the
discrete wavelet transform (DWT). The DWT decomposes, or “breaks down” a signal into
a set of orthonormal basis functions, each function being a scaled version of a “mother
wavelet”. With the CWT, at least the method used here, we relax the assumption of
orthonormality, and use a wavelet computed in the Fourier domain, i.e., by using the FFT.
NDTools uses two functions to compute the CWT: mkwobj( ) and qwave( ). The benefit
of two separate functions will become more clear when processing data with multiple
observations. The function mkwobj( ) creates a Matlab structure, or wavelet object,
containing information about the wavelet transform to be applied. The features of the
wavelet object are dependent on the length of the signal to be analyzed and the sampling
rate of the signal to be analyzed. The function generates a set of daughter wavelet
spectra, each a scaled version of the mother wavelet. The function allows construction of
N different types of wavelets: Morlet, Complex Morlet, Morlet10, Derivative of Gaussian
10
(DOG), Paul, and Mexican Hat. The examples used here will use the Morlet wavelet, which
has excellent frequency resolution and well known properties. Further, the Morlet wavelet
has been consistently used in the literature for EEG analysis, thus making it suitable for
interpreting EEG results.
The construction of the daughter wavelet spectra is determined by a pre-defined set of
wavelet scales, which mkwobj( ) computes from an array of center frequencies. Each
daughter wavelet has a “center frequency”, the frequency in the wavelet with the most
energy and that corresponds with the period of a Fourier frequency. The default frequency
array in mkwobj( ) is a linearly spaced array from 1 to 100 Hz. The inverse of the
frequency values in that array are the Fourier periods that will be used to construct each
daughter wavelet. The corresponding center frequencies that were actually generated are
returned in the wavelet object. Although the length of each wavelet is dependent on the
size of the signal to be analyzed, the length of the wavelets actually constructed will be the
next highest power of 2 above the signal length. Other features of the wavelet object be
discussed throughout the analysis.
Once a wavelet object has been constructed, the signal and the wavelet object can be
passed to the function qwave( ). This function performs the CWT on the signal using the
wavelets constructed by mkwobj( ) and returns a structure containing the wavelet
coefficients, the mean of the signal, and an indexing array used by other functions.
The following example will construct a wavelet object and perform the CWT on the
90-10-90 Hz frequency sweep from the prior example. The wavelet object will be
constructed using a Morlet wavelet.
%
First, construct a wavelet object
w=mkwobj('morl',length(signal),Fs);
%
Next, compute the wavelet coefficients:
q=qwave(signal,w);
11
The workspace will now contain the following variables:
Note that w (the wavelet object) and q (the wavelet coefficients returned by the CWT) are
both structures, which each contain useful information. The wavelet object contains the
following fields:
The wavelet object contains two basic classes of information: wavelet properties and CWT
inputs. The wavelet properties are variables that provide information about the timefrequency plane that is being analyzed. The CWT inputs are used to compute the CWT
and/or the iCWT. These variables (fields in the wavelet object) include:
12
Field in wavelet object
Description
type
String specifying the type of wavelet that was constructed
Fs
Sampling frequency of the signal to be analyzed
cF
Array of center frequencies for each wavelet scale
scales
Array of the periods for each of the wavelet scales
coi
Array defining the cone of influence
daughter
Matrix of the daughter wavelet spectra
The q-structure returned by qwave( ) contains the following fields:
Field in q-structure
Description
cfs
Matrix of the wavelet coefficients
sigm
Mean of the signal being analyzed (used for reconstruction)
oidx
Array of indices for the original signal after padding
The cfs matrix in the q-structure is the result of the CWT, and is the data to plot for
visualizing the CWT results. Because the data is a 2D matrix (note that cfs can also be 3dimensional; this will be discussed in a later section), each row will correspond to one of
the analyzed wavelet scales, and each column to a time point in the signal. Therefore, the
data should be visualized in a 3-dimensional plot.
Visualizing CWT Results
qsurf( )
NDTools contains a special plotting function designed for these data, called qsurf( ). This
function creates a 3-dimensional surface plot of the coefficients. If qsurf( ) is only given
the coefficients, then the size of each dimension (rows and columns) are used for the axes.
However, the function can also accept inputs specifying the axes. In this case, the arrays
needed for each dimension are the center frequencies for each wavelet scale (found in
13
w.cF) and an array of times corresponding with each point. In this example, the signal is
sampled at 1000 Hz, so each point in the signal represents 1 millisecond of time. If we set
the first point in the signal to time 0, then there will be 4591 milliseconds in the signal (the
length of the signal minus one). An array corresponding to these values can be generated
quite simply in Matlab:
% create an array of time points
atimes=0:1:length(signal)-1;
Next, the array of wavelet scale center frequencies (w.cF) and the new times array (atimes)
can be passed to qsurf( ) to generate a surface plot of the data. Before generating a plot,
there is an important consideration regarding how to plot the data. The wavelet
coefficients in the q-structure are complex numbers, meaning that each point in the matrix
actually contains two numbers, a real part and an imaginary part. So what should be
represented in the plot? The representation most commonly used is the complex
modulus, essentially the absolute value of a complex number, which transforms the
number into a real, positive value. If a complex-valued matrix is passed to qsurf( ), then
the modulus is automatically computed for plotting. However, if another representation is
desired, then that transform must be performed before passing the data to qsurf( ).
%
Create a new figure:
figure;
hold on;
%
plot the data; let qsurf take the absolute value of q.cfs
qsurf(w.cF,atimes,q.cfs);
% label the axes:
xlabel('points');
ylabel('pFrequency (Hz)');
title('CWT of 90-10-90 Hz Sweep');
% add a colorbar:
colorbar;
14
The figure to the left shows
time (in milliseconds) on
the abscissa and center
wavelet frequency on the
ordinate axis. The
magnitude of each
frequency activation is
shown by the color
gradient. Note that the
lower frequency
information (near the
bottom of the figure)
reveals a stronger overall
magnitude than the higher
frequencies, and this
change in magnitude appears to be negatively correlated with frequency.
CWT Normalization
qnrg( )
The CWT result might seem at odds with an expectation of equal energy at all frequencies,
given that each frequency is represented by only one cycle. To better understand this
effect, consider the content of the wavelets being used for analysis. Each daughter
wavelet is represented by its one-sided Fourier spectrum (this is slightly different when
working with complex wavelets). The CWT works by convolving each daughter wavelet
with the original signal, an operation that is performed by multiplying the FFTs of each
signal, and then taking the inverse-FFT of the result. The daughter wavelet spectra are
each designed to specifically represent the chosen wavelet, in this case the Morlet
wavelet. Consider, then, that the CWT results actually contain the original signal, plus each
of the daughter wavelets. The Matlab code below can be used to plot the daughter
wavelets for examination.
15
In the following example, two daughter wavelet spectra will be plotted: the wavelets
closest to 20 Hz and 70 Hz in center frequency.
%
Define the center frequencies to plot
dd=[20 70];
%
Find the indices for those frequencies
f1=abs(w.cF-dd(1));
f2=abs(w.cF-dd(2));
f1=find(f1==min(f1));
f2=find(f2==min(f2));
%
Retrieve the daughters:
d1=w.daughter(f1,:);
d2=w.daughter(f2,:);
%
Define the array of frequencies for representation
freq=Fs/2*linspace(0,1,length(d1)/2+1);
%
Create new figure
figure;
hold on;
%
Plot the two daughters
plot(freq,d1(1:length(freq)),'b');
plot(freq,d2(1:length(freq)),'r');
axis tight;
xlabel('Frequency (Hz)');
ylabel('Magnitude');
title('Daughter Wavelet Spectra: 20 & 70 Hz');
The plot of the wavelet spectra
reveal that the lower frequency
wavelet has a larger magnitude,
but also has a smaller bandwidth
(the range of frequencies under
the curve). This pattern exists
exists for the entire range of
scales, and reflects the tradeoff
between time and frequency
localization.
16
To view these wavelets in the time-domain, the inverse-FFT can be applied to the wavelet
spectra. To visualize the wavelet, a shift must also be applied to the result. The following
Matlab code will transform each of the wavelets plotted above, and plot them in the timedomain. The wavelet spectra are not complex, so the result of the ifft( ) will be. Therefore,
only the real part of the result will be retained.
% compute inverse-FFT
y1=ifft(d1);
y2=ifft(d2);
% apply shift
y1=ifftshift(y1);
y2=ifftshift(y2);
% retain the real part:
y1=real(y1);
y2=real(y2);
% Create new figure
figure;
hold on;
% plot the wavelets
plot(y1,'b');
plot(y2,'r');
axis tight;
xlabel('points');
ylabel('amplitude');
title('Morlet Wavelets: 20 & 70 Hz');
The plot to the left shows both
of the wavelets in the time
domain. At this scale, it is a bit
difficult to visualize either of the
wavelets. By changing the limits
on the x-axis the wavelets can
be viewed more closely.
17
The following code will rescale the x-axis.
% rescale the axis
xlim([3900 4300]);
Notice that the peak amplitude of the 70 Hz wavelet (red wave) is larger than the peak of
the 20 Hz wavelet (blue wave). Further, the 70 Hz wave appears more narrow, overall.
Examining the amplitude envelope of the wavelets reveals this tradeoff. To amplitude
envelope of a signal is the magnitude of the analytic signal. In Matlab, the analytic signal
is computed by the Hilbert transform, using the function hilbert( ). The magnitude of the
analytic signal is the modulus (the analytic signal is a complex numbered signal with both
a real and imaginary part).
% compute the amplitude envelope
env1=abs(hilbert(y1));
env2=abs(hilbert(y2));
18
The same code used to plot the wavelets can be used to plot the wavelet envelopes.
% Create new figure
figure;
hold on;
% plot the wavelets
plot(env1,'b');
plot(env2,'r');
axis tight;
xlabel('points');
ylabel('amplitude');
title('Morlet Envelopes: 20 & 70 Hz');
% rescale the axis
xlim([3900 4300]);
The above figure reveals that when in the time domain, the higher frequency signal (the 70
Hz signal plotted in red) has a larger magnitude and smaller bandwidth (the range of points
under the curve, while the lower frequency signal (20 Hz, plotted in blue) has a smaller
magnitude and larger bandwidth. Together this pattern is the inverse of the pattern found
in the frequency domain. Again, this demonstrates the tradeoff between time and
frequency localization.
19
The previous exercise was meant to demonstrate the tradeoff in time-frequency
localization. Having now been demonstrated, these properties of the wavelets can be
used to explain why higher frequencies in the sweep appeared with less energy despite
having the same number of cycles in the time domain. Recall that each daughter wavelet
is a scaled version of a mother wavelet, which implies that each daughter wavelet has the
same number of cycles. Given that period is inversely proportional to the frequency of a
signal, then it stands that one cycle of a higher frequency wavelet (e.g. 70 Hz) occupies a
smaller area of time than one cycle of a lower frequency wavelet (e.g., 20 Hz). Conversely,
one cycle of a higher frequency wavelet will occupy a smaller area of frequency than one
cycle of lower frequency wavelet.
Extended from this property of the CWT, the areas occupied in either the time or frequency
domain can serve as an energy constant, by which time and frequency information can be
represented with equivalent energies between signal components. In other words, the
energy constant can be used to normalize the signal representations. In fact, there are
several ways in which these properties can be used for normalization. In a subsequent
example, the estimated Fourier period will be used for normalization.
In the the present example, the estimated total energies for each wavelet will be used for
normalization. Essentially, the total energy for each wavelet is the area under the curve of
the wavelet spectrum, which can be computed as the sum of the values from each
daughter wavelet. When normalizing the frequencies/scales, each scale is divided by the
scale’s energy. In NDTools, the function qnrg( ) will normalize a set of wavelet coefficients
to the wavelet energies.
% normalize scales
nq=qnrg(q,w);
The output from the above function is a structure that is identical to the q-structure, but
with the wavelet coefficients normalized. Now if the normalized results are visualized
using the qsurf( ) function, the energies at each scale should be similar.
20
% create new figure
figure;
hold on;
% subplot the original coefficients
subplot(1,2,1);
qsurf(w.cF,atimes,q.cfs);
colorbar;
xlabel('points');
ylabel('Frequency (pHz)');
title('Original Coefficients');
% subplot the normalized coefficients
subplot(1,2,2);
qsurf(w.cF,atimes,nq.cfs);
colorbar;
xlabel('points');
ylabel('Frequency (pHz)');
title('Normalized Coefficients');
21
Visualizing Multiple Signal Properties
qprop( ) qprops( )
As predicted, normalizing the coefficients equates each scale to the maximum relative
energy (i.e. the maximum point in the original coefficients matrix). Because normalizing
the coefficients can affect interpretation of the data, it is good practice to consider both
the raw and normalized coefficients when determining features of a signal. To facilitate
visualization of a signals properties, NDTools includes a function called qprop( ) that takes
a signal, its q-structure, and its wavelet object (i.e., signal, q, and w) and returns a figure
with four subplots, the time-amplitude function, the frequency spectrum, the original
coefficients, and the normalized coefficients.
% qprop automatically opens a new figure
qprop(signal,q,w);
22
The next example will demonstrate a verification of the CWT and its visualized
components. That is, a verification that the representation of both time and frequency are
accurate. To do this, the properties of the wavelet energies can once again be exploited
as a normalizing function. Recall the tradeoff between the magnitudes and bandwidths
between the time and frequency domain of the wavelets, and that the energy constant can
be used to scale the data between these domains. In that case, the coefficients can be
normalized toward the time domain by the inverse of the energy constant. Then, the sum
of the real part of the normalized coefficients will represent a hyperscaled representation of
the original signal in the time domain. Likewise, the sum of the scale normalized
coeffcients will represent a (smoothed) frequency spectrum of the original signal. It is
important to note that the quality of the representation is dependent on the resolution of
the wavelet object. A function called qprops( ) operates in the same way as qprop( ) but
returns a six-paneled figure that includes the sums of the domain normalized coefficients,
as well as the sum of the raw coefficients in the frequency domain.
% qprops automatically opens a new figure
qprops(signal,q,w);
The figure returned by qprops( ) shows that time-frequency representations are relatively
accurate in both domains. The time domain appears to have some missing power around
the edges of the of waveform. The frequency domain also appears to be well represented.
In the bottom right panel, the blue line is the sum of the raw coefficients, which appears to
mirror that of the original spectrum. The red line is the sum of the normalized coefficients,
demonstrating that the energy for active frequencies has been equated.
23
The Cone of Influence
qmask( )
As demonstrated by the sums of normalized time energies, some of the energy at the
edges of the waveform appears to have been lost. Indeed, this is an effect of the wavelet
transform. The extent of information lost is dependent on the center frequency of the
wavelet scales and the length of the signal being analyzed. Specifically, this area of
24
uncertainty can be defined by conical boundaries called the “cone of influence” (COI).
It
is possible to suppress, or “turn off” the values outside the COI when visualizing the data.
This is achieved with the function qmask( ), which creates a two-dimensional mask from
the COI and applies it to the coefficients in the q-structure. The initial cone was computed
by the function mkwobj( ), and is therefore a variable within the wavelet object. By
default, qmask( ) masks the coefficients by multiplying values outside of the cone by zero.
The function also accepts an optional input to make this any value. The following Matlab
code will implement a COI mask using NaNs instead of zeros.
% mask the coefficients
mq=qmask(q,w,NaN);
% create new figure
figure;
hold on;
% subplot the original coefficients
qsurf(w.cF,atimes,nq.cfs);
xlabel('points');
ylabel('Center Frequency (Hz)');
title('Normalized & Masked Coefficients');
Although the difference is not easily noticeable in this case, the COI will be useful when
analyzing signals in a shorter time window in subsequent sections.
25
4
CWT Filtering
iqwave( )
The frequency sweep used to demonstrate the CWT is not very representative of the
signals analyzed for EEG (or many other applications). Most signals of interest in the EEG
are small, transient signals; signals that are only present for a very brief time. Most often
these signals are evoked by some known event, typically a controlled sensory or cognitive
stimulus.
To demonstrate the effects of the CWT on such signals, the following example will
demonstrate how to generate a simulated signal and analyze using the CWT. For this new
example, the sampling rate will be set to 1000 Hz, and the signal window will be 1 second
long, and will contain an 18 Hz wave with 2.5 cycles and a Hamming window function for
the amplitude envelope.
% Sampling rate
Fs=1000;
% Number of points (for 1 second)
pnts=1000;
% Generate the signal
y=mkpsin(18,Fs,2.5,pnts,'hamming');
Next, a wavelet object will be generated to match this signal’s length, and the CWT will be
computed.
% make wavelet object
w=mkwobj('morl',pnts,Fs);
% compute CWT
q=qwave(y,w);
% visualize results
qprops(y,q,w);
26
The above figure shows a near-perfect reconstruction of the signal in the time-frequency
domain. However, much like our previous examples, the signals being used have only one
frequency activation at any point in time. In most systems, there is a considerable overlap
between frequency components. Rather than demonstrate the CWT with a single
component, the following demonstrates the CWT properties of a signal with two,
overlapping components. To create the signal, a second component will be generated
using mkpsin( ) and added to the 18 Hz component, above. The new component will
have a frequency of 55 Hz with 3 cycles and a hamming window for the amplitude
envelope.
27
% create new signal
y2=mkpsin(55,Fs,3,pnts,'hamming');
% create new signal
y2=mkpsin(55,Fs,3,pnts,'hamming');
% normalize the new signal
y=qnorm(y+y2,1,'max');
% re-compute the CWT with new signal
q=qwave(y,w);
% visualize results
qprops(y,q,w);
28
The figure of the complex signal demonstrates that the two components can, indeed, be
separated in the time-frequency domain. When comparing the two activations, the
tradeoff between time and frequency localization at different scales is apparent. For
example, the higher frequency component (55 Hz) appears with a more oblong
morphology, spanning many frequencies around the center. The lower frequency
component (18 Hz) appears more conical in shape, spanning fewer frequencies. These
results also suggest that separating overlapping components very close in frequency may
prove more difficult.
Signals & Noise
mkpnoise( )
Unlike the signals analyzed above, most real data will contain a considerable amount of
noise, or unwanted frequency components. It stands to reason, then, that when
examining the effects of the CWT on a signal, the presence of noise should also be
modeled into that signal. NDtools contains a function called mkpnoise( ) that can be used
to generate a wide variety of spatial and temporal noise patterns.
mkpnoise( ) can generate 1D or 2D signals with varying noise properties. Four different
noise color patterns are available: white, pink, brown, and psuedo-white. Additionally, an
option for custom FFT filtering of uniformly distributed white noise is available. 2D spatial
patterns can be generated with periodic patterns (the pattern is similar across all rows;
good for modeling spatial EEG noise) or with aperiodic patterns (the pattern is random
across all rows). Further, the rows of noise can be spatially weighted so that the each row
contributes a different amount to the overall energy of the matrix. The output can be
specified to have an RMS amplitude of 1 for each row, or for the total matrix (spatially
weighted noise can only have a total RMS amplitude of 1).
In the following example, one second of pseudo-white noise (random phases under a
normal distribution) will be generated and added to a signal.
29
% generate one second of noise
noise=mkpnoise([1 1000],'color','white','cdist','norm');
% compute wavelet and show the properties
nq=qwave(noise,w);
qprops(noise,nq,w);
30
Re-running the previous code will generate a different result, as new random noise is
generated each time. This noise can now be mixed with a signal of the same length by
simply adding the two vectors. However, it is better to have some control over the signalto-noise ratio (SNR) before analyzing a noisy signal, resulting in a more controlled
examination. By their defaults, the function mkpsin( ) returns a signal with an amplitude of
1, and the function mkpsoise( ) returns a noise with an RMS amplitude of 1. Because the
signal is padded with zeros, the computation of the signal RMS will be skewed because of
any extra padding in the signal. It is possible to have mkpsin( ) return a signal with an
RMS amplitude of 1, which is computed and scaled before padding the output. The
following example will recreate the two component signal from before, but with each
component scaled to an RMS amplitude of 1.
Because both the signal and noise will
have equivalent RMS amplitudes, a signal with a controlled SNR can be created by scaling
the noise (or the signal) to the desired SNR before adding. The following example creates
a noisy signal with SNR=2 (signal is twice the power of the noise).
% generate one second of noise
noise=mkpnoise([1 pnts],'color','white','cdist','norm');
% Generate the two signals with RMS=1
y1=mkpsin(18,Fs,2.5,pnts,'hamming','rms');
y2=mkpsin(55,Fs,3,pnts,'hamming','rms');
% Make the signal by averaging
% (maintain mean RMS)
y=mean([y1;y2]);
% Scale the noise
% SNR = 2 (signal is 2x the noise)
snr=2;
noise=noise./snr;
% Create a noisy signal
yn=y+noise;
% compute CWT of noisy signal
q=qwave(yn,w);
% plot the properties
qprops(yn,q,w);
31
The time-amplitude waveform (top left panel) of the signal can still be visualized within the
noise, but it is clearly more difficult to make out. Likewise, the frequency spectrum (top
right panel) is becoming less distinguishable, especially for the higher frequency
component. However, the CWT seems to have performed well amongst the noise, and the
normalized time magnitudes (bottom left panel) seem to reveal the signal quite nicely. It
should be noted that the choice of CWT parameters affects the normalized time
magnitude representation; the CWT is essentially acting as a filter by not including
frequency information above ~100Hz. Similar results can be achieved by FFT filtering. The
function mfftfilt( ) will be used to filter the noisy signal for comparison.
32
% get the CWT time magnitudes
cmag=sum(real(q.cfs));
% normalize to size of ynf
cmag=cmag./max(abs(cmag(:)));
cmag=cmag-mean(cmag);
cmag=cmag.*max(ynf);
% find the min and max of CWT frequencies
f1=min(w.cF);
f2=max(w.cF);
% apply a (sharp) FFT filter in the same range
ynf=mfftfilt(yn,[0 f1 f2 f2+10],Fs,'hamming');
% plot the results and compare
figure;
hold on;
plot(ynf);
plot(cmag,'r');
axis tight;
xlabel('points');
ylabel('amplitude');
title('CWT vs FFT Filter');
The results of summing the CWT time magnitudes (red) are almost identical to those
achieved by a the FFT filter (blue). These results suggest that the CWT can be used as a
basic filter for the signals if the signal amplitudes can be reverted to the scale of the
original signal (in this case, the scales were forced to be equal by normalizing).
33
Inverse CWT
iqwave( )
As it turns out, Compos & Torrence (1997) solved for a normalizing function by correcting
each scale for the total energy before computing the sums. By adding the mean of the
original signal to the result, the time waveform is re-scaled to original units. The qstructure contains this value, which was retained during the initial CWT computation. The
function iqwave( ) can be used to perform this inverse operation, called the inverse CWT.
% compute iCWT
iq=iqwave(q,w);
% create new figure
figure;
hold on;
% plot the original
plot(yn,':b','LineWidth',0.5);
plot(iq,'r','LineWidth',1.5);
axis tight;
xlabel('points');
ylabel('amplitude');
title('Inverse CWT');
34
The above filtering was achieved with no transforms other than computing the CWT and
then computing the iCWT. Also, the noisy signal construction embedded a relatively lower
frequency signal (18 and 55 Hz) into a randomly generated noise that has energy up to
~500 Hz. In most signal processing applications, some filtering of the signal will occur
online, meaning the noise will also be filtered within some known range. Therefore, when
modeling the noise, it may be a good idea to filter the noise before computing the SNR.
% generate one second of noise
noise=mkpnoise([1 pnts],'color','white','cdist','norm');j
% filter the noise
noise=mfftfilt(noise,[0 f1 f2 f2+10],Fs,'hamming');
% rescale the RMS
nr=sqrt(mean(noise.^2));
nr=1/nr;
noise=noise.*nr;
% now make new noise with SNR=2
snr=2;
noise=noise./snr;
yn=y+noise;
% get new CWT
q=qwave(yn,w);
% visualize signal properties
qprops(yn,q,w);
The result of creating a noisy signal in-range can be seen in the properties figure
generated by the above code (see figure, next page). There is a noticeable difference
when the noise has been filtered; it is more difficult to distinguish signal components from
noise components (even at 2 to 1 SNR).
If the CWT is found to be just as susceptible to in-range noise as other methods, then
what do we really gain by using the CWT? This question will be examined in the next
chapter. The issue of low-SNR signals will be addressed again in a subsequent chapter.
35
36
5
Time-Frequency Masking
qmask( )
In the previous chapter, it was shown that the CWT acts in a very similar manner to an FFT
filter. However, some limitations of this were also demonstrated. What, then, are the
advantages to using the CWT? In that previous example, there were no further
transformations of the data before summing the coefficients. In other words, the filtering
only occurred over the frequency domain. If signal features can be identified in the timefrequency plane, and isolated over both dimensions, then the CWT would have an
advantage as a filtering approach.
To explore this possibility, a new signal with three, overlapping components will be
generated and embedded in noise. The signal components will be three frequency bursts,
each windowed with different window functions (to demonstrate the use of the amplitude
functions). Rather than the components being exactly on top of each other, they will be
slightly staggered throughout the signal.
% sampling rate
Fs=1000;
% number of points
pnts=1000;
% wavelet object
w=mkwobj('morl',pnts,Fs);
% generate three waveforms:
y1=mkpsin(38,Fs,7.5,-1,'hamming','rms');
y2=mkpsin(18,Fs,3.5,-1,'hamming','rms');
y3=mkpsin(8,Fs,3.0,-1,'hamming','rms');
% generate three empty arrays:
a1=zeros(1,pnts);
a2=zeros(1,pnts);
a3=zeros(1,pnts);
37
% first signal starts at 231 (arbitrary)
a1(1,231:231+length(y1)-1)=y1;
% second signal starts at 297 (1/3rd of the way out)
a2(1,297:297+length(y2)-1)=y2;
% third signal starts at 395 (1/2 of the way out)
a3(1,395:395+length(y3)-1)=y3;
% add them up
y=sum([a1;a2;a3]);
% clear used variables
clear y1 y2 y3 a1 a2 a3
% get cwt
q=qwave(y,w);
% visualize signal properties
qprops(y,q,w);
38
Next, a complex noise will be generated, which will be added to the signal. For this noise,
a mixture of normally distributed noise (pseudo-white) and brown noise will be filtered from
0.1 to 100 Hz, and normalized to an RMS amplitude of 1. The noise will then be added to
the signal with an SNR of 2. As we will see, this noise will begin to resemble the kind of
noise observed in EEG.
% generate noises
wnoise=mkpnoise([1 pnts],'color','white','cdist','norm');
bnoise=mkpnoise([1 pnts],'color','brown','cdist','norm');
% make a weighted mixture
wnoise=wnoise.*0.75;
bnoise=bnoise.*0.25;
noise=wnoise+bnoise;
% filter
noise=mfftfilt(noise,[0 0.1 100 120],1000,'hamming');
% scale to RMS=1
nr=sqrt(mean(noise.^2));
nr=1/nr;
noise=noise.*nr;
% visualize the noise
figure;
hold on;
plot(noise);
axis tight;
xlabel('points');
ylabel('amplitude');
title('Complex Noise');
39
% set SNR
snr=2;
% make a noisy signal
yn=y+(noise./snr);
% get CWT
q=qwave(yn,w);
% visualize properties
qprops(yn,q,w);
40
As seen by the properties figure, the CWT did a much better job of recovering the signal,
but it is also clear that the signal was not as affected by the noise, as the components are
mostly distinguishable in the time-amplitude waveform, and certainly in the frequency
spectrum. Taken together, this suggests that different types of noise can affect the signals
quite differently. Compare this result to the same signal and noise, but with an SNR of 1.
41
Even with more noise added, the CWT performs well at recovering the components.
Recall that the wavelet coefficients are actually complex numbers, with both a real and
imaginary part. The surface plots of the coefficients are plotted as the magnitudes
(absolute values) of those complex numbers. Further, the inverse CWT relies on the
amplitudes; the real part of the complex number. Therefore, it would be ideal to “highlight”
the features observed by magnitude in the raw coefficients in order to filter only those
features by suppressing the uninteresting information.
IF the features of the signal are the most salient features in the time-frequency plane, then
the wavelet power can be used to create a mask (a 2-dimensional filter). The mask can
then be applied to the raw coefficients before computing the iCWT. The result would be
that information outside of these features is suppressed by some exponential value before
computing the inverse. This can all be accomplished with the function qmask( ),
introduced earlier. In addition to masking the cone of influence, the qmask( ) function can
accept a string specifying that the mask should be the power of the time-frequency plane.
%% mask with power
mq=qmask(q,'pow');
The above code returns a structure identical to the original q-structure, but with the
coefficients (q.cfs) masked by the time-frequency power. The following code shows, stepby-step, how this mask was computed and applied.
% get the power
pow=abs(q.cfs).^2;
% normalize to unit 1
pow=pow./max(pow(:));
% pow is now the mask!
% copy q and apply the mask
mq=q;
mq.cfs=mq.cfs.*pow;
To visualize the effect of the mask, replace the original q-structure with the masked qstructure in the qprops( ) function.
42
% visualize masked properties
qprops(yn,mq,w);
The effects of masking appear to allow a much cleaner reconstruction of the original components, as seen by
the sums of time energies (lower left panel). Further, the sums of frequency magnitudes do a relatively decent
job of identifying the three frequency components, but is still dominated by the very low frequencies of the
brown noise. As noted before, this filtering approach is severely limited by the ability of CWT to construct
43
those features. The following figure shows the same procedure on the noisy signal from the previous chapter
(two components, plus pseudo-white noise with an SNR of 2.
Despite the noisier signal, the power masking method does a good job of recovering the original components.
Ideally, these properties can be further exploited to detect time-frequency components in very noisy signals
(e.g., with negative SNRs), but with multiple observations; something more like an evoked potential data set
extracted from EEG.
44
6
Signal Averaging
In many signal processing applications, the signal of interest is evoked by a stimulus and
can be recorded many times. The signal is estimated by averaging the multiple trials, such
that each contribution to the average is time-locked to the evoking stimulus. In this way,
the signal presumably appears at the same location & time each time it is evoked.
Therefore, the average of the signal will be a (mostly) stable representation of the evoked
response. The noise, which is not time-locked to the stimulus, will begin to converge
toward zero as more trials have been contributed. Because the average of the noise
approaches zero, and the signal average remains stable, an average of many trials should
reveal a signal that is likely to reflect the signal of interest. This procedure is the basis of
computing evoked and event-related potentials from multiple trials of EEG (as well as
other recordings, such as otoacoustic emissions, subterranean structure maps, etc.).
In the following example, multiple trials of an identical signal will be created, and each trial
will contain a noise mixture with an SNR of 1 (equal signal and noise). The trials and the
resulting averages will be created cumulatively, to demonstrate the convergence. Note
that when the following code is run with Matlab, a figure window will be generated, and
the cumulative average will be interactively displayed. The signal will be the same signal
used in the previous chapter, however the code to generated the signal is repeated in this
example.
% sampling rate
Fs=1000;
% number of points
pnts=1000;
% wavelet object
w=mkwobj('morl',pnts,Fs);
45
% generate three waveforms:
y1=mkpsin(38,Fs,7.5,-1,'hamming','rms');
y2=mkpsin(18,Fs,3.5,-1,'hamming','rms');
y3=mkpsin(8,Fs,3.0,-1,'hamming','rms');
% generate three empty arrays:
a1=zeros(1,pnts);
a2=zeros(1,pnts);
a3=zeros(1,pnts);
% find locations for each signal
% first signal starts at 231 (arbitrary)
a1(1,231:231+length(y1)-1)=y1;
% second signal starts at 297 (1/3rd of the way out)
a2(1,297:297+length(y2)-1)=y2;
% third signal starts at 395 (1/2 of the way out)
a3(1,395:395+length(y3)-1)=y3;
% add them up
y=sum([a1;a2;a3]);
% clear used variables
clear y1 y2 y3 a1 a2 a3
% create an empty array for the averages
psum=zeros(1,pnts);
avg=zeros(1,pnts);
% select number of trials
trials=100;
% set SNR
snr=1;
% create a new figure with handle
h=figure;
hold on;
% initiate averaging loop
for i=1:trials
% create the noise
wnoise=mkpnoise([1 pnts],'color','white','cdist','norm');
bnoise=mkpnoise([1 pnts],'color','brown','cdist','norm');
% make a weighted mixture
wnoise=wnoise.*0.75;
bnoise=bnoise.*0.25;
noise=wnoise+bnoise;
46
% filter
noise=mfftfilt(noise,[0 0.1 100 120],1000,'hamming');
% scale to RMS=1
nr=sqrt(mean(noise.^2));
nr=1/nr;
noise=noise.*nr;
% make a noisy signal
yn=y+(noise./snr);
% update the average
if i==1
% first trial, not an average
psum=yn;
avg=psum./i;
% set up the plot
hp=plot(yn);
axis tight;
limval=max(abs(yn))+1;
ylim([-limval limval]);
xlabel('points');
ylabel('amplitude');
% make a title with current N
title(sprintf('Cumulative Average: N=%i',i));
else
% add the new noisy signal
psum=psum+yn;
% get the current average
avg=psum./i;
% change data on the axes
set(hp,'YData',avg);
% update the title
title(gca,sprintf('Cumulative Average: N=%i',i));
end
% slow down the plot
pause(0.15);
end
% plot the raw signal on top
plot(y,'r');
% add a figure legend
legend('Average','Raw Signal');
% clear temporary variables
clear i hp hh wnoise bnoise noise nr limval psum
47
The above figure is the final frame of the animated signal average. Note that averaging
multiple trials does an excellent job of revealing the original signal and of improving the
SNR. However, this averaging procedure assumes that the signal components are exactly
the same each time they are present. Such a simple linear system is certainly fun to
model, but real data are rarely so precise. Another assumption is that all noise in the
signal is random; yet, most real data will also contain unwanted, non-random
components. In other words, in most situations there will be some variability in the signal
from trial to trial. The extent of one’s analysis will depend primarily on whether the
average of the signal can satisfactorily explain the experimental conditions, or whether
modeling the inter-trial variability may further explain the results. There is no correct
answer to that question; rather, the answer depends entirely on the research question and
what information is necessary to answer it. However, understanding how a signal is
effected by averaging may provide some guidance. In the following examples, several
aspects of signal averaging will be explored, including the effects of temporal and
amplitude jitters, and some ways of estimating the SNR.
48
Amplitude Jitter
Amplitude jitter refers to the variation in a signal’s amplitude over repeated observations.
Clearly, if the amplitude varies substantially, then the resulting average will be an average
representation of the signal. For the following example, assume that the signal, y, has a
mean peak-amplitude of 2 (e.g., the raw signal in used in the previous example), but the
peak amplitude of y in any random trial will fall somewhere around this mean; let’s say with
a standard deviation of 0.4. The following code will generate an array of peak-amplitude
values that follow the above rules.
% generate 100 random amplitudes
% with mean=2 and s.d.=0.4;
r=2+0.4.*randn(100,1);
% plot the distribution
figure;
hist(r,20);
xlabel('Peak Amplitude');
ylabel('Count');
title('Pseudo-randomly Distributed Amplitudes');
49
If these amplitudes represented the peak amplitudes of the signal for each of 100 trials,
then we would expect the signal average to have a peak amplitude close to 2. The
following example tests this assumption:
% create a normalized y
yt=y./max(abs(y));
% replicate yt for 100 trials
yt=repmat(yt,100,1);
% replicate the amplitudes
rt=repmat(r,1,pnts);
% scale each trial
yt=yt.*rt;
% compute the average
ya=mean(yt);
% plot
figure;
hold on;
plot(ya);
plot(y,':r','LineWidth',1.001);
xlabel('points');
ylabel('amplitude');
title('Amplitude Jittered Average');
legend('Jittered Average','Raw Signal');
50
There is no apparent difference between the amplitude jittered average and the raw signal.
If the amplitude distribution is wider, then some difference can be observed. Consider the
following average with a mean amplitude of 2 and a standard deviation of 4.
% generate 100 random amplitudes
% with mean=2 and s.d.=0.4;
r=2+3.*randn(100,1);
% can't have negative amplitudes
r=abs(r);
% plot the distribution
figure;
hist(r,20);
xlabel('Peak Amplitude');
ylabel('Count');
title('Pseudo-randomly Distributed Amplitudes');
51
Now, apply the new amplitude jitters and compare the means.
% create a normalized y
yt=y./max(abs(y));
% replicate yt for 100 trials
yt=repmat(yt,100,1);
% replicate the amplitudes
rt=repmat(r,1,pnts);
% scale each trial
yt=yt.*rt;
% compute the average
ya=mean(yt);
% plot
figure;
hold on;
plot(ya);
plot(y,'r');
xlabel('points');
ylabel('amplitude');
title('Amplitude Jittered Average');
legend('Jittered Average','Raw Signal');
52
At the end of the day, even a large amount of amplitude jitter does not negatively affect the
average (assuming the number of trials can accurately represent the signal mean).
However, this does not negate or discount the effects of amplitude variation in an
experimental condition. For example, consider a system in which the amplitude of the
signal systematically decreases over repeated observations. Computing only the mean
will not represent the systematic change in amplitude. The case of systematic deviation
will be revisited in a subsequent chapter.
Temporal Jitter
More so than amplitude jitter, temporal jitter will have a destructive effect on the signal’s
average. Temporal jitter refers to the variation in timing of a signal relative to the timing of
the evoking stimulus. Temporal jitter can be the result of natural variation in the system,
variation in the timing of the stimulus, or variation in sampling properties. To understand
how temporal jitter can be destructive, simply consider two identical sine waves, but with
opposite phases; that is, they are 180º out of phase. The sum of the two waves would
equal zero.
% create two signals, out of phase
y1=mkpsin(20,1000,1,-1,'none',1,0);
y2=mkpsin(20,1000,1,-1,'none',1,180);
% plot them
figure;
hold on;
plot(y1,'b');
plot(y2,'r');
axis tight;
xlabel('points');
ylabel('amplitude');
title('20 Hz, 180-degrees out of phase');
% plot the sum
plot(y1+y2,'m','LineWidth',2);
% add a legend
legend('y1','y2','sum');
53
In context of temporal jitter, if a signal is jittered by a time equal to one-half of the period of
the signal, then the mean would be near 0; the signals would cancel. This also implies
that different frequencies may be more effected by jitter than others, depending on the
jitter present in the system. To examine the effects of temporal jitter, a new base signal will
be generated. In this new signal, the components will be the same as the previous
example, but instead of overlapping will occur sequentially.
% Sampling rate
Fs=1000;
% Points per trial
pnts=1000;
% Number of trials
trials=100;
% make three signals
y1=mkpsin(38,Fs,7.5,-1,'hamming','rms');
y2=mkpsin(18,Fs,3.5,-1,'hamming','rms');
y3=mkpsin(8,Fs,3.0,-1,'hamming','rms');
54
% string them together
y=[zeros(1,78) y1 zeros(1,51) y2 zeros(1,51) y3 zeros(1,51)];
% plot the signal
figure;
hold on;
plot(y);
axis tight;
xlabel('points');
ylabel('amplitude');
title('3 Components: 8, 18, & 38 Hz');
The figure above shows the three, sequential signal components of varying frequency. In
the following example, a cumulative average of will be built from multiple trials of the
above signal, but with temporal jitter modeled in the trials. This jitter will be achieved by
shifting each signal by a random amount.
55
© Copyright 2024 ExpyDoc
