<div dir="ltr">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. <div><br></div><div>Figure 1: <a href="http://mikexcohen.com/figure1.png">http://mikexcohen.com/figure1.png</a></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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. </div><div><br></div><div>Figure 2: <a href="http://mikexcohen.com/figure2.png">http://mikexcohen.com/figure2.png</a></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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.)</div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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). </div><div><div><br></div><div>I hope this clears up some confusion. Let me know if you have more questions.</div><div>Mike</div><div><br></div><div><br></div><div><br></div><div><-- code starts here --></div><div><br></div><div><div>% illustration of wavelet widths in time and frequency domains</div><div>% <a href="mailto:mikexcohen@gmail.com">mikexcohen@gmail.com</a></div><div><br></div><div>%% visual example</div><div><br></div><div>% specify parameters</div><div>wfreq = 40; % in Hz</div><div>nCyc = 20/(2*pi*wfreq); % numerator is number of cycles; denominator is a scaling factor</div><div><br></div><div>% create wavelet</div><div>wavet = -1:1/srate:1;</div><div>srate = 1000;</div><div>cmw = exp(1i*2*pi*wfreq*wavet) .* exp( -.5*(wavet/nCyc).^2 );</div><div><br></div><div>% create 10-second signal</div><div>signal = randn(1,10*srate);</div><div><br></div><div>% convolution parameters</div><div>nData = length(signal);</div><div>nKern = length(wavet);</div><div>nConv = nData+nKern-1;</div><div><br></div><div>% FFT of wavelet and signal</div><div>cmwX = fft(cmw,nConv);</div><div>sigX = fft(signal,nConv);</div><div><br></div><div>% (part of) convolution in frequency domain</div><div>asX = cmwX.*sigX;</div><div><br></div><div>% frequencies in Hz (valid only up to Nyquist)</div><div>hz = linspace(0,srate,nConv);</div><div><br></div><div>% plotting</div><div>figure(1), clf</div><div>subplot(211), plot(wavet,real(cmw),'k')</div><div>title('Time domain')</div><div>xlabel('Time (sec)'), ylabel('Amplitude (a.u.)')</div><div><br></div><div>subplot(212)</div><div>plot(hz,abs(cmwX)/max(abs(cmwX)),'k','linew',2), hold on</div><div>plot(hz,abs(sigX)/max(abs(sigX)),'r')</div><div>plot(hz,abs(asX)/max(abs(asX)),'b')</div><div>title('Frequency domain')</div><div>legend({'wavelet';'signal';'multiplication'})</div><div>xlabel('Frequency (Hz)'), ylabel('Amplitude (normalized)')</div><div><br></div><div>set(gca,'xlim',[0 200])</div><div><br></div><div>%% FWHM by frequency and number of cycles</div><div><br></div><div>% define parameters</div><div>frex = linspace(5,200,100); % vary frequency in Hz</div><div>nCyx = linspace(4,40,80); % vary number-of-cycles for time-domain wavelet</div><div>fwhm = linspace(1,20,length(nCyx)); % vary FWHM of wavelet defined in the frequency domain</div><div><br></div><div>% initializations</div><div>fwhms = zeros(length(frex),length(nCyx));</div><div>sigpow = zeros(3,length(frex),length(nCyx));</div><div><br></div><div>% double-loops</div><div>for fi=1:length(frex)</div><div> for ncyci=1:length(nCyx)</div><div> </div><div> %% create wavelet in time domain</div><div> </div><div> cmw = exp(1i*2*pi*frex(fi)*wavet) .* exp( -.5*(wavet/ (nCyx(ncyci)/(2*pi*frex(fi))) ).^2 );</div><div> cmwX = fft(cmw,nConv);</div><div> </div><div> %% compute empirical FWHM of wavelet</div><div> </div><div> cmwXp = abs(cmwX.^2)/max(abs(cmwX).^2); % must be [0 1] normalized for this algorithm to work</div><div> idx = dsearchn(hz',frex(fi)); % peak frequency</div><div> fwhms(fi,ncyci) = hz(idx-1+dsearchn(cmwXp(idx:end)',.5)) - hz(dsearchn(cmwXp(1:idx)',.5));</div><div> </div><div> %% signal power</div><div> </div><div> % non-normalized wavelet</div><div> as = ifft( cmwX.*sigX );</div><div> sigpow(1,fi,ncyci) = mean( as.*conj(as) );</div><div> </div><div> % normalized wavelet</div><div> cmwX = cmwX./max(cmwX);</div><div> as = ifft( cmwX.*sigX );</div><div> sigpow(2,fi,ncyci) = mean( as.*conj(as) );</div><div> </div><div> %% now define the Gaussian in the frequency domain</div><div> </div><div> s = fwhm(ncyci)*(2*pi-1)/(4*pi); % normalized width</div><div> x = hz-frex(fi); % shifted frequencies</div><div> fx = exp(-.5*(x/s).^2); % gaussian</div><div> fx = fx./max(fx); % gain-normalize</div><div> </div><div> sigpow(3,fi,ncyci) = mean(2*abs(sigX.*fx).^2);</div><div> </div><div> end</div><div>end</div><div><br></div><div>% and plot these spectra</div><div>figure(2), clf</div><div><br></div><div>% FWHM of wavelets defined in the time domain</div><div>subplot(221)</div><div>contourf(nCyx,frex,log(fwhms),40,'linecolor','none'), axis square</div><div>xlabel('Number of cycles'), ylabel('Frequency (Hz)')</div><div>title('Wavelet log(FWHM) in frequency domain')</div><div><br></div><div>% signal power from normalized wavelet</div><div>subplot(222)</div><div>contourf(nCyx,frex,log(squeeze(sigpow(2,:,:))),40,'linecolor','none'), axis square</div><div>xlabel('Number of cycles'), ylabel('Frequency (Hz)')</div><div>title('log(sigpow), norm.')</div><div><br></div><div>% signal power from non-normalized wavelet</div><div>subplot(223)</div><div>contourf(nCyx,frex,log(squeeze(sigpow(1,:,:))),40,'linecolor','none'), axis square</div><div>xlabel('Number of cycles'), ylabel('Frequency (Hz)')</div><div>title('log(sigpow), nonnorm.')</div><div><br></div><div>% FWHM of wavelets defined in the frequency domain</div><div>subplot(224)</div><div>contourf(fwhm,frex,log(squeeze(sigpow(3,:,:))),40,'linecolor','none'), axis square</div><div>xlabel('FWHM'), ylabel('Frequency (Hz)')</div><div>title('log(sigpow), F.D. Gaus.')</div><div><br></div><div>%% end.</div><div><br></div></div><div><br></div><div><br></div><div><br></div><div><br></div></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Fri, Aug 19, 2016 at 8:56 PM, Makoto Miyakoshi <span dir="ltr"><<a href="mailto:mmiyakoshi@ucsd.edu" target="_blank">mmiyakoshi@ucsd.edu</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr">Dear Mike,<div><br></div><div>After reading the material Andreas referred to, now I'm not sure about this:</div><span class=""><div><br></div><div>> Particularly with a constant time-domain Gaussian width, the wavelet gets wider in the frequency domain with increasing frequency.<br></div><div><br></div></span><span class=""><div>> it is due to the increasing width of the wavelet in the frequency domain.<br></div><div><br></div></span><div>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?</div><span class="HOEnZb"><font color="#888888"><div><br></div><div>Makoto</div></font></span><div class="gmail_extra"><br><div class="gmail_quote"><span class="">On Wed, Aug 17, 2016 at 11:05 AM, Mike X Cohen <span dir="ltr"><<a href="mailto:mikexcohen@gmail.com" target="_blank">mikexcohen@gmail.com</a>></span> wrote:<br></span><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left-width:1px;border-left-color:rgb(204,204,204);border-left-style:solid;padding-left:1ex"><div><div class="h5"><div dir="ltr"><div>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. <br></div><div><br></div><div>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).</div><div><br></div><div>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.</div><div><br></div><div>f=[25 60];</div><div>signal = sin(2*pi.*linspace(f(1),f(2)*m<wbr>ean(f)/f(2),length(t)).*t)';</div><div><br></div><div>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.<br></div><div><br></div><div>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.</div><div><br></div><div>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.</div><div><br></div><div>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.<span><font color="#888888"><br></font></span></div><span><font color="#888888"><div><br></div><div>Mike</div><div><br></div>-- <br><div data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><font color="#9900ff">Mike X Cohen, PhD<br><a href="http://mikexcohen.com" target="_blank">mikexcohen.com</a></font></div></div></div></div>
</font></span></div>
<br></div></div><span class="">______________________________<wbr>_________________<br>
Eeglablist page: <a href="http://sccn.ucsd.edu/eeglab/eeglabmail.html" rel="noreferrer" target="_blank">http://sccn.ucsd.edu/eeglab/ee<wbr>glabmail.html</a><br>
To unsubscribe, send an empty email to <a href="mailto:eeglablist-unsubscribe@sccn.ucsd.edu" target="_blank">eeglablist-unsubscribe@sccn.uc<wbr>sd.edu</a><br>
For digest mode, send an email with the subject "set digest mime" to <a href="mailto:eeglablist-request@sccn.ucsd.edu" target="_blank">eeglablist-request@sccn.ucsd.e<wbr>du</a><br></span></blockquote></div><span class=""><br><br clear="all"><div><br></div>-- <br><div data-smartmail="gmail_signature"><div dir="ltr">Makoto Miyakoshi<br>Swartz Center for Computational Neuroscience<br>Institute for Neural Computation, University of California San Diego<br></div></div>
</span></div></div>
</blockquote></div><br><br clear="all"><div><br></div>-- <br><div class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><font color="#9900ff">Mike X Cohen, PhD<br><a href="http://mikexcohen.com" target="_blank">mikexcohen.com</a></font></div></div></div></div>
</div>