[Eeglablist] Units in output of timefreq - wavelet normalization
Mike X Cohen
mikexcohen at gmail.com
Fri Aug 19 14:37:25 PDT 2016
Hi Makoto and list. I'm sorry if I wasn't clear enough in that email. I'll
try to explain it better, this time using better code, pictures, and a
longer email. And for the lazy readers or those without easy access to
Matlab, you can follow along the linked pictures, although I recommend
using the code so you can change parameters and see their effects.
Figure 1: http://mikexcohen.com/figure1.png
Let's first look at this qualitatively. The first cell of code will produce
a Morlet wavelet at a frequency and number-of-cycles that you specify (in
the first picture, set to 40 Hz and 20 cycles). In the lower plot, you see
the power spectrum of random noise (red; curiously enough called 'signal'),
the wavelet (black), and their point-wise multiplication (blue). The blue
trace reflects the convolution between the wavelet and the signal. The
normal procedure from here would be to take the inverse Fourier transform
of the blue line; the phases allow reconstruction of the time course and
that's how we get back a time series of frequency-band-specific activity.
If this were a filter-Hilbert procedure, then the black line would be the
power spectrum of the filter kernel, the IFFT would be real-valued, and
then you would apply the Hilbert transform to add the phase quadrature,
which gives the analytic (complex) signal. (Technically, FIR filtering
doesn't work exactly like this, but it's conceptually similar.) For STFFT,
you would just compute the average of the blue trace (not going back to the
time domain), but you wouldn't use the entire time series for the FFT, just
a piece of it (e.g,. 500 ms). Note that in this figure I've normalized all
power spectra to facilitate visualization.
Now the interactive part: Change the frequency of the wavelet (don't change
the number of cycles) and re-run. You'll see that as the wavelet creeps up
in frequency, it gets wider. That's what I meant before.
I hope this gives a qualitative idea of how the wavelet gets wider in the
frequency domain when the time-domain wavelet increases in frequency. Now
let's look more quantitatively.
Figure 2: http://mikexcohen.com/figure2.png
The second cell in the code goes through two loops that vary the frequency
and the number-of-cycles in the wavelet, and will compute the FWHM of the
wavelet spectrum (in the frequency domain), and the total power in the
result of convolution between each set of wavelet parameters and the
signal. Let's go through each plot in turn.
Top row: Given a fixed number of time-domain-defined cycles, the wavelet
power spectrum gets wider with increasing frequency. Consequently, the
total power in the convolved signal also increases. In fact, this is a
trivial result: in this case, because our noise signal has a flat spectrum,
the upper-right plot is really just the same thing as the upper-left plot
plus some noise. To link this back to Figure 1, each pixel in the
upper-left plot is the full-width-at-half-maximum (FWHM) of the black line
in the lower panel of Figure 1. And each pixel in the upper-right plot is
the average of the blue line in the lower panel of Figure 1. You can see in
the upper plots of Figure 2 that in order to maintain constant power (that
is, to follow a ray of solid coloring from the origin going out), you would
need to increase the number-of-cycles as the frequency increases. That's
why this is standard procedure in time-frequency analysis.
So is this an artifact? No, it's not. It's just the way the Fourier
transform works. Trying to correct for this while keeping the number of
cycles constant is actually underestimating power at higher frequencies,
because it involves trying to make the red bars more blue, even though the
redness is accurate. (But note that this has no effect on dB or % change,
because the underestimation affects the baseline period equally.)
Now for the lower-left plot. This one shows the same analysis as the
upper-right plot (signal power for each parameter pair), except the wavelet
has *not* been max-value normalized. Not only is the entire pattern
completely different, but the values are way different, overestimating the
power by many many orders of magnitude (you can see this by adding a
colorbar; I didn't do this because only the upper-left and lower-right
panels have values that are directly comparable). Yikes. But again, this
overestimation affects the baseline power as well, so dB or % change is
completely blind to this.
And I'm sure you are curious what the lower-right plot shows. It turns out
that a simple solution to this conundrum is to define the wavelet FWHM in
the frequency domain, not in the time domain. That produces vertical bars,
in other words, the wavelet width (and therefore also signal power) is
fixed over all frequencies for a given FWHM. This plot shows the signal
power, not the width, which is why it's a bit noisy. Also note that the
Gaussian in the code is defined differently: in the frequency domain, time
is inverted (e.g., Hz = 1/s), so the conversion from FWHM to Gaussian width
works differently than in the time domain. I mention this explicitly
because if you want to implement this and try using the "standard" Gaussian
function, it will give incorrect results.
So why don't we always define wavelets in the frequency domain? I actually
don't know. I guess it's just an accident of history, that things started
off in the time domain and we keep doing it that way just cuz. It does seem
more intuitive to describe wavelet widths in terms of Hz rather than in
number-of-cycles, but that's been the standard practice in the field. I use
the frequency-domain implementation when I want to specify wavelet width in
Hz (it's easier than starting from number-of-cycles), but for all of my
"normal" TF analyses, I still use the good old time-domain implementation
with the number-of-cycles parameter. Again, if you are applying a baseline
normalization, none of this matters at all (well, the width of the wavelet
definitely matters; the raw power values don't matter).
I hope this clears up some confusion. Let me know if you have more
questions.
Mike
<-- code starts here -->
% illustration of wavelet widths in time and frequency domains
% mikexcohen at gmail.com
%% visual example
% specify parameters
wfreq = 40; % in Hz
nCyc = 20/(2*pi*wfreq); % numerator is number of cycles; denominator is
a scaling factor
% create wavelet
wavet = -1:1/srate:1;
srate = 1000;
cmw = exp(1i*2*pi*wfreq*wavet) .* exp( -.5*(wavet/nCyc).^2 );
% create 10-second signal
signal = randn(1,10*srate);
% convolution parameters
nData = length(signal);
nKern = length(wavet);
nConv = nData+nKern-1;
% FFT of wavelet and signal
cmwX = fft(cmw,nConv);
sigX = fft(signal,nConv);
% (part of) convolution in frequency domain
asX = cmwX.*sigX;
% frequencies in Hz (valid only up to Nyquist)
hz = linspace(0,srate,nConv);
% plotting
figure(1), clf
subplot(211), plot(wavet,real(cmw),'k')
title('Time domain')
xlabel('Time (sec)'), ylabel('Amplitude (a.u.)')
subplot(212)
plot(hz,abs(cmwX)/max(abs(cmwX)),'k','linew',2), hold on
plot(hz,abs(sigX)/max(abs(sigX)),'r')
plot(hz,abs(asX)/max(abs(asX)),'b')
title('Frequency domain')
legend({'wavelet';'signal';'multiplication'})
xlabel('Frequency (Hz)'), ylabel('Amplitude (normalized)')
set(gca,'xlim',[0 200])
%% FWHM by frequency and number of cycles
% define parameters
frex = linspace(5,200,100); % vary frequency in Hz
nCyx = linspace(4,40,80); % vary number-of-cycles for time-domain
wavelet
fwhm = linspace(1,20,length(nCyx)); % vary FWHM of wavelet defined in the
frequency domain
% initializations
fwhms = zeros(length(frex),length(nCyx));
sigpow = zeros(3,length(frex),length(nCyx));
% double-loops
for fi=1:length(frex)
for ncyci=1:length(nCyx)
%% create wavelet in time domain
cmw = exp(1i*2*pi*frex(fi)*wavet) .* exp( -.5*(wavet/
(nCyx(ncyci)/(2*pi*frex(fi))) ).^2 );
cmwX = fft(cmw,nConv);
%% compute empirical FWHM of wavelet
cmwXp = abs(cmwX.^2)/max(abs(cmwX).^2); % must be [0 1] normalized
for this algorithm to work
idx = dsearchn(hz',frex(fi)); % peak frequency
fwhms(fi,ncyci) = hz(idx-1+dsearchn(cmwXp(idx:end)',.5)) -
hz(dsearchn(cmwXp(1:idx)',.5));
%% signal power
% non-normalized wavelet
as = ifft( cmwX.*sigX );
sigpow(1,fi,ncyci) = mean( as.*conj(as) );
% normalized wavelet
cmwX = cmwX./max(cmwX);
as = ifft( cmwX.*sigX );
sigpow(2,fi,ncyci) = mean( as.*conj(as) );
%% now define the Gaussian in the frequency domain
s = fwhm(ncyci)*(2*pi-1)/(4*pi); % normalized width
x = hz-frex(fi); % shifted frequencies
fx = exp(-.5*(x/s).^2); % gaussian
fx = fx./max(fx); % gain-normalize
sigpow(3,fi,ncyci) = mean(2*abs(sigX.*fx).^2);
end
end
% and plot these spectra
figure(2), clf
% FWHM of wavelets defined in the time domain
subplot(221)
contourf(nCyx,frex,log(fwhms),40,'linecolor','none'), axis square
xlabel('Number of cycles'), ylabel('Frequency (Hz)')
title('Wavelet log(FWHM) in frequency domain')
% signal power from normalized wavelet
subplot(222)
contourf(nCyx,frex,log(squeeze(sigpow(2,:,:))),40,'linecolor','none'), axis
square
xlabel('Number of cycles'), ylabel('Frequency (Hz)')
title('log(sigpow), norm.')
% signal power from non-normalized wavelet
subplot(223)
contourf(nCyx,frex,log(squeeze(sigpow(1,:,:))),40,'linecolor','none'), axis
square
xlabel('Number of cycles'), ylabel('Frequency (Hz)')
title('log(sigpow), nonnorm.')
% FWHM of wavelets defined in the frequency domain
subplot(224)
contourf(fwhm,frex,log(squeeze(sigpow(3,:,:))),40,'linecolor','none'), axis
square
xlabel('FWHM'), ylabel('Frequency (Hz)')
title('log(sigpow), F.D. Gaus.')
%% end.
On Fri, Aug 19, 2016 at 8:56 PM, Makoto Miyakoshi <mmiyakoshi at ucsd.edu>
wrote:
> Dear Mike,
>
> After reading the material Andreas referred to, now I'm not sure about
> this:
>
> > Particularly with a constant time-domain Gaussian width, the wavelet
> gets wider in the frequency domain with increasing frequency.
>
> > it is due to the increasing width of the wavelet in the frequency domain.
>
> This seems valid for the case of short-time Frourier Transform, for
> example. However, in our example code, the wavelet being used changes
> Gaussian width, right? Then the 'Heisenberg box' changes its shape
> accordingly, and you can't say that the wavelet lets through more signal
> because narrower time width... am I wrong? What do you mean by this example?
>
> Makoto
>
> On Wed, Aug 17, 2016 at 11:05 AM, Mike X Cohen <mikexcohen at gmail.com>
> wrote:
>
>> Hi everyone. I agree with Andreas that normalization is a tricky issue
>> and, to some extent, a philosophical one. In general, I recommend against
>> any interpretation of "absolute" values, because (1) they depend on a
>> variety of uninteresting factors like electrode montage, equipment, filter
>> characteristics, and so on, (2) they are entirely incomparable across
>> methods. You can compare dB or % change between EEG, MEG, and LFP, but it
>> is impossible to compare EEG microvolts with LFP microvolts, MEG teslas,
>> change in light fluorescence, etc.
>>
>> I point this out because I think we have here mainly an academic
>> discussion for the vast majority of neuroscience research, particularly for
>> any neuroscience researchers that hope to link their findings to other
>> pockets of neuroscience regardless of montage, species, decade, etc. That
>> said, if there's one thing academics love, it's an academic discussion ;)
>> so here are my two cents (the Dutch don't use pennies, so you'll have to
>> decide whether to round down to zero or up to .05 euros).
>>
>> From Andreas' code, you can add the following two lines after "signal,"
>> which will make a new signal, a chirp. You can then add colorbars to both
>> TF plots to see that the power is accurately reconstructed after max-val
>> normalization. The two numbers in variable f are for the start and end
>> frequencies of the linear chirp.
>>
>> f=[25 60];
>> signal = sin(2*pi.*linspace(f(1),f(2)*mean(f)/f(2),length(t)).*t)';
>>
>> The next point concerned the increase in power over frequency. This is a
>> feature, not a bug. First of all, it is highly dependent on the number of
>> cycles. For example, note that the power in the top-middle plot goes up to
>> just over .2. Now change the 'cycles' parameter to 30; the power now goes
>> up to around .05. In other words, the horrible linear increase was cut to a
>> quarter. A constant number of cycles over a large range of frequencies is a
>> poor choice of parameter, and it should come as no surprise that poor
>> parameter choices lead to poor results.
>>
>> So why does this even happen? Particularly with a constant time-domain
>> Gaussian width, the wavelet gets wider in the frequency domain with
>> increasing frequency. This means that more of the signal is being let
>> through the filter. More signal = more power. I do not see how this is an
>> artifact, or even a problem. The more of the spectrum you look at, the more
>> power you will see. If you want to maximize the power, then use the entire
>> spectrum. In fact, total FFT power is the same as total time-domain power,
>> so the most power you can get from the FFT will be sum(signal.^2), which is
>> a lot more than what you'd get from any wavelet.
>>
>> In other words, the increase in power with increasing frequency is *not*
>> due to increasing frequency; it is due to the increasing width of the
>> wavelet in the frequency domain. This seems worse for white noise because
>> of the flat spectrum, but it will be less noticeable for real brain
>> signals, which have 1/f^c shape (whether EEG is broadband and noisy depends
>> very much on the characteristics of the signal one is investigating). And
>> again, this also depends on the wavelet width parameter.
>>
>> I'll conclude by reiterating that interpreting any "absolute" voltage
>> value should be avoided whenever possible. Of course, there is always the
>> occasional exception, but I think we can all agree that we should focus
>> more on effect sizes rather than on arbitrary values. Some kind of baseline
>> normalization is almost always best, and really the best way to make sure
>> your findings can be compared across the growing span of brain imaging
>> techniques in neuroscience.
>>
>> Mike
>>
>> --
>> Mike X Cohen, PhD
>> mikexcohen.com
>>
>> _______________________________________________
>> Eeglablist page: http://sccn.ucsd.edu/eeglab/eeglabmail.html
>> To unsubscribe, send an empty email to eeglablist-unsubscribe at sccn.uc
>> sd.edu
>> For digest mode, send an email with the subject "set digest mime" to
>> eeglablist-request at sccn.ucsd.edu
>>
>
>
>
> --
> Makoto Miyakoshi
> Swartz Center for Computational Neuroscience
> Institute for Neural Computation, University of California San Diego
>
--
Mike X Cohen, PhD
mikexcohen.com
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160819/6a2921b7/attachment.html>
More information about the eeglablist
mailing list