Multitaper<\/h4>\n
We can sacrifice time or frequency resolution for SNR. I.e. one can average estimates at neighbouring timewindows or estimates at neighbouring frequencies. A smart way to do this is to use multitapers, I will not go into details here. In the end, they allow us to average smartly over time or frequency.<\/p>\n
Parameter one usually specifies<\/h4>\n\n- foi: The frequencies of interest. Note that you can use arbitrary high resolution but if you have $\\delta foi < \\delta f$ then information will be displayed redundantly (but it looks smoother).<\/li>\n
- tapsmofrq: How much frequency smoothing should happen at each foi<\/em>.<\/li>\n
- t_ftimwin: The timewindow for each foi<\/em><\/li>\n<\/ul>\n
Plotting the parameters<\/h3>\n
With the matlab-function attached to the end of the post one can see the effective frequency-resolution of the parameters in a single plot. For this example I use the frequency parameters by Hipp 2011<\/a>. They logarithmically scale up the frequencybandwith. In order to do so for the low frequencies (1-16Hz) they change the size of the timewindow, for the higher frequencies they use multitaper with a constant window size of 250ms.<\/p>\n![\"hipp_theoretical\"](\"http:\/\/benediktehinger.de\/blog\/science\/upload\/sites\/2\/2016\/10\/hipp_theoretical-e1478469252720.jpg\")
<\/a>
\nThis is the theoretical<\/em> resolution (interactive graph<\/a>), i.e. a plot of the parameters you put in. foi +- tapsmofrq scaled by each timewindow size. Note: The individual bar size depicts the window-size of the respective frequency-bin.<\/p>\n![\"hipp_numeric\"](\"http:\/\/benediktehinger.de\/blog\/science\/upload\/sites\/2\/2016\/10\/hipp_numeric.jpg\")
<\/a>
\nThis is the numeric<\/em> resolution (interactice graph<\/a>). I used the frequency response of 100 repetitions of white noise (flat spectrum over repetitions) and calculated the frequency response of the multitaper filters. This is scaled in the x-axis by the time-window used.<\/p>\nCheck out the interactive graphs to see how the time-windows change as intendet with lower frequencies.<\/p>\n
The colorbar depict the number of tapers used. As intended, the number of tapers do not change from 1-16Hz.<\/p>\n
Code<\/h3>\nfunction [] = plot_taper(inpStruct,cfg)\n%% plot_taper\n% inpStruct: All fields need to have the same length\n% foi : The frequencies of interest. \n% tapsmofrq : How much frequency smoothing should happen at each foi (+-tapsmofrq)\n% t_ftimwin : The timewindow for each foi\n% You can also directly input a fieldtrip time-freq object\n%\n% cfg:\n% logy : plot the y-axis as log10\n% type : \"theoretical\" : This plots the foi+-tapsmofrq\n% \"numerical\" : Simulation of white noise to approx the\n% effective resolution\n\nif nargin == 1\n cfg = struct();\nend\n\nif isfield(cfg,'logy')\n assert(ismember(cfg.type,[0,1]),'logy has to be 0 or 1')\n logy = cfg.logy;\nelse \n logy = 0;\nend\nif isfield(cfg,'type')\n assert(ismember(cfg.type,{'numerical','theoretical'}),'type has to be ''numerical'' or ''theoretical''')\n type = cfg.type;\n\nelse\n type = 'numerical'; \nend\nnRep = 100; % repetitions of the numeric condition\nxscaling = 0.5;\nif ~isempty(inpStruct) && isfield(inpStruct,'t_ftimwin')\n fprintf('cfg detected\\n')\n cfg = inpStruct;\nelseif ~isempty(inpStruct) && isfield(inpStruct,'cfg')\n fprintf('fieldtrip structure detected\\n')\n cfg = inpStruct.cfg;\nelse\n fprintf('using parameters from Hipp 2010 Neuron \\n')\n cfg =[];\n\n % we use logspaced frequencies from 4 to 181 Hz\n cfg.foi = logspace(log10(4),log10(181),23);\n\n % The windows should have 3\/4 octave smoothing in frequency domain\n cfg.tapsmofrq = (cfg.foi*2^(3\/4\/2) - cfg.foi*2^(-3\/4\/2)) \/2; % \/2 because fieldtrip takes +- tapsmofrq\n\n % The timewindow should be so, that for freqs below 16, it results in n=1\n % Taper used, but for frequencies higher, it should be a constant 250ms.\n % To get the number of tapers we use: round(cfg.tapsmofrq*2.*cfg.t_ftimwin-1)\n cfg.t_ftimwin = [2.\/(cfg.tapsmofrq(cfg.foi<16)*2),repmat(0.25,1,length(cfg.foi(cfg.foi>=16)))];\nend\n\nmax_tapers = ceil(max(cfg.t_ftimwin.*cfg.tapsmofrq*2));\ncolor_tapers = parula(max_tapers);\n%%\nfigure\nfor fIdx = 1:length(cfg.t_ftimwin)\n timwin = cfg.t_ftimwin(fIdx);\n tapsmofrq = cfg.tapsmofrq(fIdx); % this is the fieldtrip standard, +- the taper frequency\n freqoi = cfg.foi(fIdx);\n switch type\n case 'theoretical'\n % This part simply plots the parameters as requested by the\n % user.\n y = ([freqoi-tapsmofrq freqoi+tapsmofrq freqoi+tapsmofrq freqoi-tapsmofrq]);\n\n ds = 0.3;\n ys = ([freqoi-ds freqoi+ds freqoi+ds freqoi-ds]);\n if logy\n y = log10(y);\n ys = log10(ys);\n end\n x = [-timwin\/2 -timwin\/2 +timwin\/2 +timwin\/2]+fIdx*xscaling;\n\n patch(x,y,color_tapers(round(2*timwin*tapsmofrq),:),'FaceAlpha',.7)\n patch(x,ys,[0.3 1 0.3],'FaceAlpha',.5)\n case 'numerical'\n fsample = 1000; % sample Fs \n % This part is copied together from ft_specest_mtmconvol\n timwinsample = round(timwin .* fsample); % the number of samples in the timewin\n\n % get the tapers\n tap = dpss(timwinsample, timwinsample .* (tapsmofrq .\/ fsample))';\n tap = tap(1:(end-1), :); % the kill the last one because 'it is weird'\n\n % don't fully understand, something about keeping the phase\n anglein = (-(timwinsample-1)\/2 : (timwinsample-1)\/2)' .* ((2.*pi.\/fsample) .* freqoi);\n\n % The general idea (I think) is that instead of windowing the\n % time-signal and shifting the wavelet per timestep, we can multiply the frequency representation of\n % wavelet and the frequency representation of the data. this\n % masks all non-interesting frequencies and we can go back to\n % timespace and select the portions of interest. In formula:\n % ifft(fft(signal).*fft(wavelet))\n %\n % here we generate the wavelets by windowing sin\/cos with the\n % tapers. We then calculate the spectrum\n wltspctrm = [];\n for itap = 1:size(tap,1)\n coswav = tap(itap,:) .* cos(anglein)';\n sinwav = tap(itap,:) .* sin(anglein)';\n wavelet = complex(coswav, sinwav);\n % store the fft of the complex wavelet\n wltspctrm(itap,:) = fft(wavelet,[],2);\n end\n\n\n% % We do this nRep times because a single random vector does not\n% have a flat spectrum.\n datAll = nan(nRep,timwinsample);\n for rep = 1:nRep\n sig = randi([0 1],timwinsample,1)*2-1; %http:\/\/dsp.stackexchange.com\/questions\/13194\/matlab-white-noise-with-flat-constant-power-spectrum\n spec = bsxfun(@times,fft(sig),wltspctrm'); % multiply the spectra \n datAll(rep,:) = sum(abs(spec),2); %save only the power\n end\n dat = mean(datAll); %average over reps\n\n\n % normalize the power\n dat = dat- min(dat);\n dat = dat.\/max(dat);\n\n t0 = 0 + fIdx*xscaling;\n t = t0 + dat*timwin\/2;\n tNeg = t0 - dat*timwin\/2;\n\n f = linspace(0,fsample-1\/timwin,timwinsample);%0:1\/timwin:fsample-1;\n if logy\n f = log10(f);\n end\n % this is a bit hacky violin plot, but looks good.\n patch([t tNeg(end:-1:1)],[f f(end:-1:1)],color_tapers(2*round(timwin*tapsmofrq),:),'FaceAlpha',.7)\n end\n\n\nend\nc= colorbar;colormap('parula')\nc.Label.String = '# Tapers Used';\ncaxis([1 size(color_tapers,1)])\n\nset(gca,'XTicks',[0,1])\nxlabel('width of plot depicts window size in [s]')\nif logy\n set(gca,'ylim',[0 log(max(1.2*(cfg.foi+cfg.tapsmofrq)))]);\n ylabel('frequency [log(Hz)]')\nelse\n set(gca,'ylim',[0 max(1.2*(cfg.foi+cfg.tapsmofrq))]);\n ylabel('frequency [Hz]')\nend\nend\n\n<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"Goal Give a graphical representation of multitaper time frequency parameters. Background Frequency Resolution and sFFT In M\/EEG analysis we are often interested in oscillations. One tool to analyse oscillations is by using time-frequency methods. In principle they work by segmenting the signal of interest in small windows and calculating a spectrum for each window. A tutorial can be found on the fieldtrip page. There exists a time-frequency tradeof: $\\delta f = \\frac{1}{T}$ with T the timewindow (in s) and $\\delta f$ the frequency resolution. The percentage-difference of two neighbouring frequency resolution gets smaller the higher the frequency goes, e.g. a…<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5],"tags":[],"class_list":["post-159","post","type-post","status-publish","format-standard","hentry","category-blog"],"_links":{"self":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/159","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/comments?post=159"}],"version-history":[{"count":0,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/159\/revisions"}],"wp:attachment":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/media?parent=159"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/categories?post=159"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/tags?post=159"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}
Plotting the parameters<\/h3>\n
With the matlab-function attached to the end of the post one can see the effective frequency-resolution of the parameters in a single plot. For this example I use the frequency parameters by Hipp 2011<\/a>. They logarithmically scale up the frequencybandwith. In order to do so for the low frequencies (1-16Hz) they change the size of the timewindow, for the higher frequencies they use multitaper with a constant window size of 250ms.<\/p>\n Check out the interactive graphs to see how the time-windows change as intendet with lower frequencies.<\/p>\n The colorbar depict the number of tapers used. As intended, the number of tapers do not change from 1-16Hz.<\/p>\n Goal Give a graphical representation of multitaper time frequency parameters. Background Frequency Resolution and sFFT In M\/EEG analysis we are often interested in oscillations. One tool to analyse oscillations is by using time-frequency methods. In principle they work by segmenting the signal of interest in small windows and calculating a spectrum for each window. A tutorial can be found on the fieldtrip page. There exists a time-frequency tradeof: $\\delta f = \\frac{1}{T}$ with T the timewindow (in s) and $\\delta f$ the frequency resolution. The percentage-difference of two neighbouring frequency resolution gets smaller the higher the frequency goes, e.g. a…<\/p>\n","protected":false},"author":2,"featured_media":0,"comment_status":"open","ping_status":"open","sticky":false,"template":"","format":"standard","meta":{"footnotes":""},"categories":[5],"tags":[],"class_list":["post-159","post","type-post","status-publish","format-standard","hentry","category-blog"],"_links":{"self":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/159","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/users\/2"}],"replies":[{"embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/comments?post=159"}],"version-history":[{"count":0,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/posts\/159\/revisions"}],"wp:attachment":[{"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/media?parent=159"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/categories?post=159"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/benediktehinger.de\/blog\/science\/wp-json\/wp\/v2\/tags?post=159"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}<\/a>
\nThis is the theoretical<\/em> resolution (interactive graph<\/a>), i.e. a plot of the parameters you put in. foi +- tapsmofrq scaled by each timewindow size. Note: The individual bar size depicts the window-size of the respective frequency-bin.<\/p>\n<\/a>
\nThis is the numeric<\/em> resolution (interactice graph<\/a>). I used the frequency response of 100 repetitions of white noise (flat spectrum over repetitions) and calculated the frequency response of the multitaper filters. This is scaled in the x-axis by the time-window used.<\/p>\nCode<\/h3>\n
function [] = plot_taper(inpStruct,cfg)\n%% plot_taper\n% inpStruct: All fields need to have the same length\n% foi : The frequencies of interest. \n% tapsmofrq : How much frequency smoothing should happen at each foi (+-tapsmofrq)\n% t_ftimwin : The timewindow for each foi\n% You can also directly input a fieldtrip time-freq object\n%\n% cfg:\n% logy : plot the y-axis as log10\n% type : \"theoretical\" : This plots the foi+-tapsmofrq\n% \"numerical\" : Simulation of white noise to approx the\n% effective resolution\n\nif nargin == 1\n cfg = struct();\nend\n\nif isfield(cfg,'logy')\n assert(ismember(cfg.type,[0,1]),'logy has to be 0 or 1')\n logy = cfg.logy;\nelse \n logy = 0;\nend\nif isfield(cfg,'type')\n assert(ismember(cfg.type,{'numerical','theoretical'}),'type has to be ''numerical'' or ''theoretical''')\n type = cfg.type;\n\nelse\n type = 'numerical'; \nend\nnRep = 100; % repetitions of the numeric condition\nxscaling = 0.5;\nif ~isempty(inpStruct) && isfield(inpStruct,'t_ftimwin')\n fprintf('cfg detected\\n')\n cfg = inpStruct;\nelseif ~isempty(inpStruct) && isfield(inpStruct,'cfg')\n fprintf('fieldtrip structure detected\\n')\n cfg = inpStruct.cfg;\nelse\n fprintf('using parameters from Hipp 2010 Neuron \\n')\n cfg =[];\n\n % we use logspaced frequencies from 4 to 181 Hz\n cfg.foi = logspace(log10(4),log10(181),23);\n\n % The windows should have 3\/4 octave smoothing in frequency domain\n cfg.tapsmofrq = (cfg.foi*2^(3\/4\/2) - cfg.foi*2^(-3\/4\/2)) \/2; % \/2 because fieldtrip takes +- tapsmofrq\n\n % The timewindow should be so, that for freqs below 16, it results in n=1\n % Taper used, but for frequencies higher, it should be a constant 250ms.\n % To get the number of tapers we use: round(cfg.tapsmofrq*2.*cfg.t_ftimwin-1)\n cfg.t_ftimwin = [2.\/(cfg.tapsmofrq(cfg.foi<16)*2),repmat(0.25,1,length(cfg.foi(cfg.foi>=16)))];\nend\n\nmax_tapers = ceil(max(cfg.t_ftimwin.*cfg.tapsmofrq*2));\ncolor_tapers = parula(max_tapers);\n%%\nfigure\nfor fIdx = 1:length(cfg.t_ftimwin)\n timwin = cfg.t_ftimwin(fIdx);\n tapsmofrq = cfg.tapsmofrq(fIdx); % this is the fieldtrip standard, +- the taper frequency\n freqoi = cfg.foi(fIdx);\n switch type\n case 'theoretical'\n % This part simply plots the parameters as requested by the\n % user.\n y = ([freqoi-tapsmofrq freqoi+tapsmofrq freqoi+tapsmofrq freqoi-tapsmofrq]);\n\n ds = 0.3;\n ys = ([freqoi-ds freqoi+ds freqoi+ds freqoi-ds]);\n if logy\n y = log10(y);\n ys = log10(ys);\n end\n x = [-timwin\/2 -timwin\/2 +timwin\/2 +timwin\/2]+fIdx*xscaling;\n\n patch(x,y,color_tapers(round(2*timwin*tapsmofrq),:),'FaceAlpha',.7)\n patch(x,ys,[0.3 1 0.3],'FaceAlpha',.5)\n case 'numerical'\n fsample = 1000; % sample Fs \n % This part is copied together from ft_specest_mtmconvol\n timwinsample = round(timwin .* fsample); % the number of samples in the timewin\n\n % get the tapers\n tap = dpss(timwinsample, timwinsample .* (tapsmofrq .\/ fsample))';\n tap = tap(1:(end-1), :); % the kill the last one because 'it is weird'\n\n % don't fully understand, something about keeping the phase\n anglein = (-(timwinsample-1)\/2 : (timwinsample-1)\/2)' .* ((2.*pi.\/fsample) .* freqoi);\n\n % The general idea (I think) is that instead of windowing the\n % time-signal and shifting the wavelet per timestep, we can multiply the frequency representation of\n % wavelet and the frequency representation of the data. this\n % masks all non-interesting frequencies and we can go back to\n % timespace and select the portions of interest. In formula:\n % ifft(fft(signal).*fft(wavelet))\n %\n % here we generate the wavelets by windowing sin\/cos with the\n % tapers. We then calculate the spectrum\n wltspctrm = [];\n for itap = 1:size(tap,1)\n coswav = tap(itap,:) .* cos(anglein)';\n sinwav = tap(itap,:) .* sin(anglein)';\n wavelet = complex(coswav, sinwav);\n % store the fft of the complex wavelet\n wltspctrm(itap,:) = fft(wavelet,[],2);\n end\n\n\n% % We do this nRep times because a single random vector does not\n% have a flat spectrum.\n datAll = nan(nRep,timwinsample);\n for rep = 1:nRep\n sig = randi([0 1],timwinsample,1)*2-1; %http:\/\/dsp.stackexchange.com\/questions\/13194\/matlab-white-noise-with-flat-constant-power-spectrum\n spec = bsxfun(@times,fft(sig),wltspctrm'); % multiply the spectra \n datAll(rep,:) = sum(abs(spec),2); %save only the power\n end\n dat = mean(datAll); %average over reps\n\n\n % normalize the power\n dat = dat- min(dat);\n dat = dat.\/max(dat);\n\n t0 = 0 + fIdx*xscaling;\n t = t0 + dat*timwin\/2;\n tNeg = t0 - dat*timwin\/2;\n\n f = linspace(0,fsample-1\/timwin,timwinsample);%0:1\/timwin:fsample-1;\n if logy\n f = log10(f);\n end\n % this is a bit hacky violin plot, but looks good.\n patch([t tNeg(end:-1:1)],[f f(end:-1:1)],color_tapers(2*round(timwin*tapsmofrq),:),'FaceAlpha',.7)\n end\n\n\nend\nc= colorbar;colormap('parula')\nc.Label.String = '# Tapers Used';\ncaxis([1 size(color_tapers,1)])\n\nset(gca,'XTicks',[0,1])\nxlabel('width of plot depicts window size in [s]')\nif logy\n set(gca,'ylim',[0 log(max(1.2*(cfg.foi+cfg.tapsmofrq)))]);\n ylabel('frequency [log(Hz)]')\nelse\n set(gca,'ylim',[0 max(1.2*(cfg.foi+cfg.tapsmofrq))]);\n ylabel('frequency [Hz]')\nend\nend\n\n<\/code><\/pre>\n","protected":false},"excerpt":{"rendered":"