Transcript Spectra of random processes

Spectra of random processes
Signal, noise, smoothing and
filters
Gaussian processes are defined
by…
• Their mean at each time, m(t)
• Their covariance matrix, Cm(t,s)
Stationary Gaussian processes are
defined by…
• Their mean at each time, m
• Their covariance function, C(t)
– Note that t is now a time interval, not an
absolute time
– Cm(t,s)=C(t-s)
Stationary processes are defined
by…
• Mean and covariance are not sufficient in
general
• Let’s keep that in mind for a while…
The covariance function
• We assume that the process is Gaussian
and stationary, and that its mean is 0
MATLAB
Spectral representation of a
random process
• We have our samples of the random
process
• We can transform each of them into the
frequency domain
• We can collect now all the values of the
frequency components at a specific
frequency
• These are samples of a random variable
whose distribution can be computed
Spectral representation of a
random process
• We do it because the spectral
representation may be often simpler than
the time representation
• Frequency components may be far less
correlated than time samples
Go to Matlab
Spectral representations
• What we saw is that the process can be
rather nicely described in the following
way:
– Select real and imaginary parts,
independently and normally, with variance
that is frequency-dependent
– Synthesize the appropriate time sequence
Spectral representations
• It turns out that all Gaussian processes
follow this prescription
• The function that gives the variance as a
function of frequency is called the spectral
density function of the process
Spectral representations
• And if the process is stationary but nongaussian?
• The frequency components are still
uncorrelated, but they are not normal any
more
• Therefore, there may dependence
between them (but no correlation!)
Spectral representations
• Non-stationary processes
– Correlations between frequency components, nonzero means
• Stationary processes
– No correlations between frequency components, zero
means
– Dependences between frequency components
• Gaussian processes
– Independence and normality in the frequency domain
Spectral density functions and
autocorrelations
• The spectral density function of a
stationary process is the Fourier transform
of its (non-normalized) autocorrelation
function
– Autocorrelations are very special: they have a
real positive Fourier transform
Go To Matlab
Spectral density functions and
autocorrelations
• Why use spectral density rather than
autocorrelation?
– Because Fourier components of stationary
processes are uncorrelated and even
independent (in the gaussian case)
– Because of dynamic range issues
– Because the spectral density function is a
variance decomposition, good for any
stationary process…
Spectral density functions and
autocorrelations
• The spectral density function can be seen as a
variance decomposition of the total energy of the
process
• Proof:
– The energy of the process is the (non-normalized)
autocorrelation at lag 0
– This is the 0-th component of the inverse Fourier
transform of the spectral density
– Which is the DC of the spectral density, i.e. the sum of
the variance of all spectral components
Additive models
• Usually, the signal is not
stationary/Gaussian because there is a
special shape in it (an action potential,
evoked potential in the LFP, a crosscorrelation peak, …)
• We will usually assume that the shape is
exactly the same every time, and that
whatever is not shape is gaussian noise
Additive models
• The independence of frequency
components in Gaussian processes is the
frequency domain expression of their lack
of structure
• Generally, we will consider Gaussian
contributions to our signal as noise
• Sometimes this is not true – it’s the
variance that is the signal
– EMG
Getting rid of noise
• If we know that our signal doesn’t have
much energy in a certain frequency range,
we can filter out that range of frequencies
without hurting much the signal
• By this process, we reduce the energy of
the noise more than we reduce the energy
of the signal, improving signal to noise
ratio
Filtering in the frequency domain
• What we would like to do is to set all
frequency components outside the
interesting range to 0
• We will see that this is not necessarily the
best way of doing it, but let’s proceed…
Go To Matlab
Filtering in the frequency domain
• Two problems:
– Large fluctuations (when signal is transient)
– Non-causality
• Partial solution to non-causality – zero-padding
Filtering in the frequency domain
• Fluctuations are due to the introduction of
sharp discontinuities into the Fourier
domain
• Time domain then becomes wide (high
‘frequencies’)
• We can try smoother cutoffs
Go To Matlab
Filtering in the frequency domain
• Decide which frequencies should be
suppressed (‘stop band’)
• The other frequencies are called the pass
band
• Build a window that is 1 inside the pass
band, 0 inside the stop band, and
interpolates between 0 and 1 gradually
• Do it!
• Beware of edge effects…
Filtering in the frequency domain
• The fine point here is the selection of the
window to use for multiplying the spectra
Filtering in the frequency domain
• Additive model:
– s(t)=x(t)+n(t)
– x(t), n(t) are stationary processes and independent of each other with
known autocorrelation functions
• In that case, the constraint of making s(t) as similar as possible to x(t)
can be formalized to get an optimal filter (Wiener filtering)
• The multiplier is roughly
1
1  N( f ) / X ( f )
•
•
•
•
When X(f)>>N(f), this is about 1 (passband)
When X(f)<<N(f), this is about 0 (stopband)
Smooth interpolation between the two extremes
The same formula works for constant x(t)
Go To Matlab
Filtering in the time domain
• We already did this to remove highfrequency noise:
– Select a short window
– Take a short piece of signal, average, replace
the middle point of the piece of signal by the
new point
– This is non-causal!
Filtering in the time domain
• Causality:
– We want to implement such filtering on-line
– This means that at each point in time, we
know the current sample and all past ones,
but not future ones
Filtering in the time domain
• Causal time-domain filtering:
– Select a short window
– Take a short piece of signal, average, replace the
last(!) point of the piece of signal by the new point
• The choice of the window will determine what we
do
– Lowpass filtering (weighted averages)
– Highpass filtering (highly fluctuating weights)
– Bandpass filtering
• (MATLAB)
Filtering in the time domain
• Problems:
– Edge effects
• The implementation assume that the signal is 0
before it begins
– Time shifts
• Impulses peak later (causality)
– At the moment, we have only a rough idea
what this is doing in the frequency domain
Time-domain filtering and sine
waves
• We want to understand what time-domain
filtering does to the frequency components
of a signal
• We will start by testing the effect of a timedomain filter on sine waves
• (MATLAB)
Time-domain filtering and sine
waves
• When putting in a sine wave, a timedomain filter will output a sine wave (up to
edge effects)
• The output will have the same frequency
as the input
• The output will differ in amplitude and
phase
Time-domain filtering and sine
waves
• The change in amplitude and the shift in
phase depend only on the frequency of the
input
• This can be concisely summarized by a
complex number (‘gain’)
Time-domain filtering and sine
waves
• The complex gain is in fact the Fourier
transform of the time-domain window
• Proof: (whiteboard)
Time-domain filtering and general
signals
• If we know what a time-domain filter does
to sine waves, we know what it does to
any other signal
• The reason is linearity
Time-domain filtering and general
signals
• Linearity:
• For two signals,
– First summing and then filtering
– First filtering and then summing
– Will give the same result
Time-domain filtering and general
signals
• We know that Fourier synthesis (inverse
Fourier transform) represents a stimulus
as a sum of sine waves
• The effect of a time-domain filter is then
the sum of the effects on each frequency
component
Time-domain filtering and general
signals
• Do a Fourier transform, to get amplitude
and phases of each frequency component
• Change the amplitude and phase of each
frequency component as implied from the
Fourier transform of the filter
• Do an inverse Fourier transform
Time-domain filtering and general
signals
• Note that the frequency domain step is
equivalent to a multiplication by a (possibly
complex) window
Comparing frequency domain and
time domain filtering
• Time domain:
• Do an FT
• Multiply by a window
(FT of time-domain
window)
• Do an iFT
• Frequency domain:
• Do an FT
• Multiply by an
(arbitrary) window
• Do an iFT
Comparing frequency domain and
time domain filtering
• Frequency domain and time domain
filtering are in fact the same thing (almost,
up to edge effects – see soon)
• The time-domain window used for filtering
and the frequency-domain window
multiplying the frequency representation
are Fourier-transform pairs
• (MATLAB)
Comparing frequency domain and
time domain filtering
• Why do we have a difference?
• The reason is the assumptions about the
edges
Comparing frequency domain and
time domain filtering
• In the time-domain filter, we assumed that before
stimulus start, there were zeros
• In the frequency-domain operation, it is implicitly
assumed that the stimulus is periodic
– Before stimulus start, the last samples come from the
end of the stimulus
• Solution: increase stimulus duration, padding
with 0’s
• (MATLAB)
Comparing frequency domain and
time domain filtering
• The same equivalence occurs when
starting with a frequency-domain window
and going back to a time domain window
• Need to keep the right symmetry of the
frequency-domain window in order to get a
real-valued time-domain window!
• (MATLAB)
Comparing frequency domain and
time domain filtering
• Why do our time-domain filters shift the
signal whereas our frequency-domain
filters keep peaks in place?
• The property of keeping peaks in place is
called ‘zero phase’
– Because all frequency components ‘stay in
place’ (the phase of the complex gain is 0)
• The result is a non-causal time-domain
window, as we saw previously
Comparing frequency domain and
time domain filtering
• We shifted the resulting time-domain
windows in order to see their actual shape
• This is equivalent to shifting the location of
filtered peaks
Non-causal time-domain filtering
• We can do non-causal filtering by first
doing a causal filtering, then delaying back
the signal to compensate for the time shift
• Or by filtering the signal twice: once going
forward (introducing a delay) and then
filtering the signal reversed in time
(introducing an equivalent delay in the
other direction)
• (MATLAB)
Nomenclature
• The time-domain window is called the
impulse response
• The frequency-domain window is called
the transfer function
Practical issues
• EDGE EFFECTS!
– Know what you are assuming
• Time domain or frequency domain?
– Causal vs. non-causal
– Computational efficacy
– Shape of time-domain filters
– Shape of frequency-domain windows
Practical issues
• EDGE EFFECTS!
– For frequency domain filtering, pad to double
the size in order to avoid circular leaks
– Don’t generate extra jumps – make sure the
average of the padded region is the same as
that of the true data
– Keep only the relevant part (first half of the
buffer)
Practical issues
• EDGE EFFECTS!
– For time domain filtering, Matlab assumes
zero initial conditions
– This may not be appropriate when the
stimulus has a non-zero average (DC)
– The function filter has explicit control of initial
conditions, see me for details
Practical issues
Causal vs. non-causal
• The literature is biased towards causal
filtering because of real-time applications
• We usually want non-causal filtering
• Either shift back the signal to compensate
for delays, or use filtfilt that does backand-forth filtering
Practical issues
Computational efficiency
• Time-domain filtering is more efficient
when the stimulus is much longer than the
window length
• Issue for very long stimuli (e.g. filtering
millions of samples)
Practical issues
Shape of time-domain windows
• You usually want your time-domain window short
– For computational efficiency (minor concern!)
– For not mixing up far segments of the stimulus (Major
concern)
• You usually want your time-domain windows
symmetric
– For keeping the general shape of peaks
– These are called ‘linear phase filters’
– When using filtfilt, the window is always effectively
symmetric
Practical issues
Shape of frequency-domain windows
• You usually want your frequency-domain
windows very sharp, 1-0 type
– To reject as fully as possible the noise
Practical issues
Incompatibility of time-domain constraints
and frequency-domain constraints
• Short time-domain filters have gradual,
non-zero frequency-domain windows
• Sharp, 1-0 frequency domain windows
have very long time-domain windows
• Must compromise!
Practical issues
• Filter design tools
– FDATool
– Fir1 + kaiserord