[FieldTrip] Simulate data to compare methods

Tom Holroyd tomh at kurage.nimh.nih.gov
Fri Jan 30 18:40:51 CET 2015 

This is more about the Subject than about filtering, but may I say
yay multitapering, and also yay Stockwell transforms. The latter are
somewhat easier to understand than wavelets, and the phase is easier to
extract. Also if you sum across time and inverse FFT the result is the
usual power specrtum. Enough about that.

Here is what I use to simulate MEG data. It's written in Python,
but it's pretty easy to translate. It creates a 1/f^2 noise and
then performs a fractional derivative to create a 1/f noise. The
noise demonstrates growth of variance over time but is nevertheless
normally distributed (mean is removed and s.d. = 1).

It makes good surrogate MEG data, properly scaled. Adding a couple ECDs
is beyond the scope of this post. :-)


from numpy import zeros, array, arange
from numpy.fft import fft, ifft
from numpy.random import normal

def meg_noise(l, n = .5):
    """Return l samples of 1/f noise."""

    l2 = l / 2

    d = zeros((l,), 'f')
    y = 0.
    for i in range(l):
        x = normal()    # white
        y += x          # brown
        d[i] = y

    # detrend
    d = d - arange(l) * (d[-1] - d[0]) / l

    # Fractional derivative of d. Regular derivative (n=1) adds 2 to the
    # exponent of the spectrum. Fractional derivative does a multiple
    # of that, so n = .5 adds 1 to the exponent. Thus for brown (-2)
    # you get pink (-1).

    w = array(range(l2) + range(-l2, 0))
    # now w = [ 0, 1, ..., l2 - 1, -l2, -l2 + 1, ..., -1 ]
    jwn = pow((1j) * w, n)
    D = fft(d)
    D = D * jwn
    dd = ifft(D).real / l

    dd -= dd.mean()
    dd /= dd.std()

    return dd



On Fri, 30 Jan 2015 13:34:16 +0100
"Jörn M. Horschig" <jorn at artinis.com> wrote:

> Hi Todor,
> 
> maybe this matlab function helps illustrating what dpss multitapers
> are, and will thus clarify what makes them so powerful compared to
> other techniques:
> https://www.dropbox.com/s/0uifk9l8rb6m5vl/Tapering.m?dl=0 (go to
> example 5)
> 
> Best,
> Jörn
> 
> 
> 
> --
> 
> Jörn M. Horschig, Software Engineer
> Artinis Medical Systems  |  +31 481 350 980 
> 
> > -----Original Message-----
> > From: fieldtrip-bounces at science.ru.nl [mailto:fieldtrip-
> > bounces at science.ru.nl] On Behalf Of Eelke Spaak
> > Sent: Friday, January 30, 2015 11:52 AM
> > To: FieldTrip discussion list
> > Subject: Re: [FieldTrip] Simulate data to compare methods
> > 
> > Hi Todor,
> > 
> > Although your procedure would also yield smoothing in the frequency
> > domain which is independent from that in the time domain, it is not
> > at all equivalent to using multitapers!
> > 
> > The tapers in the discrete prolate spheroidal sequence (dpss, ==
> > multitaper in fieldtrip) are pairwise orthogonal, hence their
> > estimates are independent from one another. This will result in
> > there being more information extracted from the signal than if you
> > used a single taper and then apply Gaussian smoothing over
> > frequencies. You could have a look at
> > https://en.wikipedia.org/wiki/Multitaper which gives quite a decent
> > overview of multitapering. Or for the full details, refer to the
> > original paper by David Thompson:
> > http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=1456701
> > 
> > Best.
> > Eelke
> > 
> > On 30 January 2015 at 11:10, tjordanov at besa.de <tjordanov at besa.de>
> > wrote:
> > > Hi Eelke,
> > >
> > > I found your answer very interesting. If I understand you
> > > correctly, the
> > advantage of the multitaper method is that it smoothes in the
> > frequency domain independently of the smoothing in the time domain.
> > Then it should be equivalent (or similar) with the following
> > procedure:
> > > 1) Calculate single trial single taper decomposition of the
> > > signal. 2) Choose an appropriate 1D Gauss function (note that it
> > > is important to be 1D else it would smooth in both - time and
> > > frequency) 3) Apply the selected Gauss function on the decomposed
> > > signal only in the
> > frequency direction (not in time in order to avoid temporal
> > smearing). Do this for all trials and all time points.
> > > 4) Calculate the average over the trials.
> > > In this procedure the choice of the Gaussian would determine the
> > > amount
> > of smearing in the frequency domain.
> > >
> > > Is this so, or I misunderstood something?
> > >
> > > Best,
> > > Todor
> > >
> > >
> > > -----Original Message-----
> > > From: fieldtrip-bounces at science.ru.nl
> > > [mailto:fieldtrip-bounces at science.ru.nl] On Behalf Of Eelke Spaak
> > > Sent: Mittwoch, 28. Januar 2015 12:24
> > > To: FieldTrip discussion list
> > > Subject: Re: [FieldTrip] Simulate data to compare methods
> > >
> > > Hi Nico,
> > >
> > > As for question (2), you probably first need to think about what
> > > constitutes
> > a "better" result. Using more tapers with dpss will always result
> > in more frequency smoothing. If your source signal is primarily
> > composed of pure sinusoids, and you interpret a spectrum as "better"
> > > if it shows clearer peaks, then you will always get the "best"
> > > result for the
> > single-taper case.
> > >
> > > Multitapering allows optimal control over the amount of smoothing
> > > you
> > obtain in the frequency domain, which is more or less independent
> > of the amount of smoothing you obtain in the time domain (as
> > opposed to e.g. wavelets, where these are fundamentally linked).
> > When dealing with brain signals, you will often find that a certain
> > stimulus might induce e.g. a gamma response at 40-50 Hz in one
> > subject and one trial, while in another subject or another trial
> > the same stimulus might induce a 50-60 Hz response or so. Of
> > course, in the average over trials (and subjects), this
> > heterogeneity (i.e., noise) will wash out, but it will severely
> > damage your statistical sensitivity. Therefore, using multitapers
> > to add smoothing can produce a much more consistent result and
> > therefore be "better" in the sense of actually understanding the
> > brain.
> > >
> > > As for your simulation, perhaps using filtered noise would be
> > > better than
> > sinusoids. Also, since multitapering benefits you most strongly
> > when taking variation over observations into account, you could
> > consider simulating different observations, each consisting of
> > noise filtered in a slightly different randomly chosen bandwidth,
> > and inspecting the resulting variation over observations in the
> > spectra.
> > >
> > > Best,
> > > Eelke
> > >
> > > On 27 January 2015 at 18:36, Max Cantor <mcantor at umich.edu> wrote:
> > >> Hi Nico,
> > >>
> > >> I'm not sure about the second question, but as for the first
> > >> question, you can manually set the scales for ft_singleplotTFR
> > >> using
> > cfg.zlim.
> > >>
> > >> Hope that helps,
> > >>
> > >> Max
> > >>
> > >> On Tue, Jan 27, 2015 at 11:50 AM, Nico Weeger
> > >> <nico.weeger at googlemail.com>
> > >> wrote:
> > >>>
> > >>> Hello FieldTrip community,
> > >>>
> > >>>
> > >>>
> > >>> I am new to FieldTrip and I try to simulate data to compare the
> > >>> ft_frequanalysis methods Hanning, Multitaper and Wavelet.
> > >>>
> > >>> Therefore I simulate Data manually using different latency,
> > >>> amplitude and frequency combinations using the following
> > >>> equation:
> > >>>
> > >>> sig1 = exp(-(t-lat1).^2/(2*sigma1))*amp1.*sin(2*pi*f1*t);
> > >>>
> > >>> sig2 = exp(-(t-lat2).^2/(2*sigma2))*amp2.*sin(2*pi*f2*t);
> > >>>
> > >>> sig3 = exp(-(t-lat1).^2/(2*sigma1))*amp1.*sin(2*pi*f2*t);
> > >>>
> > >>> sig4 = exp(-(t-lat2).^2/(2*sigma2))*amp2.*sin(2*pi*f1*t);
> > >>>
> > >>> sig = sig1+sig2+sig3+sig4;
> > >>>
> > >>> where amp1=20; amp2 = 5; lat1= 1.7; lat2 = 2.3; f1 = 12; f2 =
> > >>> 60;
> > >>>
> > >>>
> > >>> After using ft_frequanalysis (see the following cfgs)
> > >>>
> > >>>
> > >>> Cfg Wavelet:
> > >>>
> > >>> cfg = [];
> > >>>
> > >>> cfg.output     = 'pow';
> > >>>
> > >>> cfg.channel    = labels;
> > >>>
> > >>> cfg.method     = 'wavelet';
> > >>>
> > >>> cfg.width      = 7;
> > >>>
> > >>> cfg.gwidth     = 3;
> > >>>
> > >>> cfg.foilim     = [1 70];
> > >>>
> > >>> cfg.toi        = 0:0.05:2;
> > >>>
> > >>> TFRwave = ft_freqanalysis(cfg, data_preproc);
> > >>>
> > >>>
> > >>>
> > >>> Cfg Hanning / Multitaper:
> > >>>
> > >>> cfg = [];
> > >>>
> > >>> cfg.output     = 'pow';
> > >>>
> > >>> cfg.channel    = labels;
> > >>>
> > >>> cfg.method     = 'mtmconvol'
> > >>>
> > >>> cfg.foi        = 1:1:70
> > >>>
> > >>> cfg.tapsmofrq  = 0.2*cfg.foi;
> > >>>
> > >>> cfg.taper      = 'dpss' / ‘hanning’;
> > >>>
> > >>> cfg.t_ftimwin  = 4./cfg.foi;
> > >>>
> > >>> cfg.toi        = 0:0.05:2;
> > >>>
> > >>> TFRmult1 = ft_freqanalysis(cfg, data_preproc);
> > >>>
> > >>>
> > >>>
> > >>>
> > >>> the data is plotted with ft_singleplotTFR (see cfg below)
> > >>>
> > >>>
> > >>> cfg singleplot:
> > >>>
> > >>> cfg = [];
> > >>>
> > >>> cfg.maskstyle    = 'saturation';
> > >>>
> > >>> cfg.colorbar     = 'yes';
> > >>>
> > >>> cfg.layout       = 'AC_Osc.lay';
> > >>>
> > >>> ft_singleplotTFR(cfg, TFRwave);
> > >>>
> > >>>
> > >>> Two problems occur. First, the power scale of wavelet and
> > >>> Multitaper/Hanning differs extremely (Multi 0-~100 and Wavelet
> > >>> 0-
> > ~15*10^4).
> > >>>
> > >>> 1.       How can I get the scale of all methods equal, or do I
> > >>> have to change the Wavelet settings to get the right scale of
> > >>> the values?
> > >>>
> > >>> Second, the best result of Multitaper analysis is performed
> > >>> using only one Taper. The goal was to get a result, where the
> > >>> advantages and disadvantages of Multitaper analysis compared to
> > >>> the other
> > methods can be seen.
> > >>>
> > >>> 2.       How can I change the simulation so that more tapers
> > >>> show better results than a single taper does?
> > >>>
> > >>>
> > >>> Thank you for your time and help.
> > >>>
> > >>>
> > >>> Regards,
> > >>>
> > >>>
> > >>>
> > >>> Nicolas Weeger
> > >>>
> > >>> Student of Master-Program Appied Research,
> > >>>
> > >>> University Ansbach,
> > >>>
> > >>> Germany
> > >>>
> > >>>
> > >>> _______________________________________________
> > >>> fieldtrip mailing list
> > >>> fieldtrip at donders.ru.nl
> > >>> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> > >>
> > >>
> > >>
> > >>
> > >> --
> > >> Max Cantor
> > >> Lab Manager
> > >> Computational Neurolinguistics Lab
> > >> University of Michigan
> > >
> > > _______________________________________________
> > > fieldtrip mailing list
> > > fieldtrip at donders.ru.nl
> > > http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> > >
> > >
> > > _______________________________________________
> > > fieldtrip mailing list
> > > fieldtrip at donders.ru.nl
> > > http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> > 
> > _______________________________________________
> > fieldtrip mailing list
> > fieldtrip at donders.ru.nl
> > http://mailman.science.ru.nl/mailman/listinfo/fieldtrip
> 
> 
> _______________________________________________
> fieldtrip mailing list
> fieldtrip at donders.ru.nl
> http://mailman.science.ru.nl/mailman/listinfo/fieldtrip



-- 
Dr. Tom
--
I would dance and be merry,
Life would be a ding-a-derry,
If I only had a brain.
        -- The Scarecrow
-------------- next part --------------
A non-text attachment was scrubbed...
Name: GrowthOfVariance.png
Type: image/png
Size: 28854 bytes
Desc: not available
URL: <http://mailman.science.ru.nl/pipermail/fieldtrip/attachments/20150130/83159027/attachment-0002.png>

More information about the fieldtrip mailing list