[Eeglablist] Units in output of timefreq - wavelet normalization
Makoto Miyakoshi
mmiyakoshi at ucsd.edu
Fri Aug 19 17:17:43 PDT 2016
Dear Mike,
Thanks for such a detailed explanation with a very nice demonstration. It
was fun to follow. I see it clearer now.
> 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.
Personally I find it is most reasonable. I'll talk about it with the
developers.
Makoto
On Fri, Aug 19, 2016 at 2:37 PM, Mike X Cohen <mikexcohen at gmail.com> wrote:
> 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
>
--
Makoto Miyakoshi
Swartz Center for Computational Neuroscience
Institute for Neural Computation, University of California San Diego
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20160819/237784ff/attachment.html>
More information about the eeglablist
mailing list