Phalow_freqhigh

% method 'phalow_freqhigh' creates a frequency modulated signal
%   signal is build up of 3 components; s1, s2 and noise.
%     s1: represents the base signal that will be modulated
%     s2: signal that will be used for the frequency modulation
% and the output consists of the following channels:
%     1st channel: mixed signal = s1.ampl * cos(ins_pha) + noise
%     2nd channel: s1
%     3rd channel: s2
%     4th channel: noise
%     5th channel: inst_pha_base   instantaneous phase of the high (=base) frequency signal s1
%     6th channel: inst_pha_mod    low frequency phase modulation, this is equal to s2
%     7th channel: inst_pha        instantaneous phase, i.e. inst_pha_base + inst_pha_mod

cfg = [];
cfg.fsample   = 1000;
cfg.method    = 'phalow_freqhigh';
cfg.trllen    = 10;
cfg.numtrl    = 10;
cfg.s1.freq   = 20; %base frequency
cfg.s1.phase  = 0;
cfg.s1.ampl   = 1;
cfg.s2.freq   = 2;
cfg.s2.phase  = -0.5 * pi; %then base freq at t=0
cfg.s2.ampl   = pi; % should not be too high, then diff inst phase becomes negative or very close to zero and time stops or reverses
cfg.noise.ampl = 0;
cfg.output = 'all';

data = ft_freqsimulation(cfg);

figure;
sel = 1:1000;
subplot(3,3,1); plot(data.trial{1}(1,sel)); title(data.label{1})
subplot(3,3,2); plot(data.trial{1}(2,sel)); title(data.label{2})
subplot(3,3,3); plot(data.trial{1}(3,sel)); title(data.label{3})
subplot(3,3,4); plot(data.trial{1}(4,sel)); title(data.label{4})
subplot(3,3,5); plot(data.trial{1}(5,sel)); title(data.label{5})
subplot(3,3,6); plot(data.trial{1}(6,sel)); title(data.label{6})
subplot(3,3,7); plot(data.trial{1}(7,sel)); title(data.label{7})
print -dpng phalow_freqhigh_fig1.png

figure;plot(diff(data.trial{1}(7,:))); title('diff inst phase, should not be less than or close to zero')
print -dpng phalow_freqhigh_fig2.png

% show powerspectrum simulated data
cfg = [];
cfg.method    = 'mtmfft';
cfg.channel   = 'mix';
cfg.output    = 'pow';
cfg.taper     = 'hanning';
cfg.foilim    = [2 60];

fft_data = ft_freqanalysis(cfg,data);
figure; ft_singleplotER([],fft_data);
print -dpng phalow_freqhigh_fig3.png

% mtmconvol 1ex
cfg = [];
cfg.method    = 'mtmconvol';
cfg.channel   = 'mix';
cfg.output    = 'pow';
cfg.taper     = 'hanning';
cfg.foi       = 2:2:60;
cfg.toi       = data.time{1}(3001:7000);
cfg.t_ftimwin = 4./cfg.foi;
cfg.keeptrials = 'yes';

freq1 = ft_freqanalysis(cfg,data);
figure; imagesc(freq1.time(1:1000), freq1.freq, squeeze(freq1.powspctrm(1,1,:,1:1000)))
axis xy
print -dpng phalow_freqhigh_fig4.png

% plot power at different frequencies
pow_10Hz = squeeze(freq1.powspctrm(1,1,5,:));
pow_12Hz = squeeze(freq1.powspctrm(1,1,6,:));
pow_14Hz = squeeze(freq1.powspctrm(1,1,7,:));
pow_16Hz = squeeze(freq1.powspctrm(1,1,8,:));
pow_18Hz = squeeze(freq1.powspctrm(1,1,9,:));
pow_20Hz = squeeze(freq1.powspctrm(1,1,10,:));
pow_22Hz = squeeze(freq1.powspctrm(1,1,11,:));
pow_24Hz = squeeze(freq1.powspctrm(1,1,12,:));
pow_26Hz = squeeze(freq1.powspctrm(1,1,13,:));
pow_28Hz = squeeze(freq1.powspctrm(1,1,14,:));
pow_30Hz = squeeze(freq1.powspctrm(1,1,15,:));
figure;
subplot(3,4,1); plot(freq1.time(1:1000),pow_10Hz(1:1000));title('pow @ 10 Hz');ylim([0 0.4])
subplot(3,4,2); plot(freq1.time(1:1000),pow_12Hz(1:1000));title('pow @ 12 Hz');ylim([0 0.4])
subplot(3,4,3); plot(freq1.time(1:1000),pow_14Hz(1:1000));title('pow @ 14 Hz');ylim([0 0.4])
subplot(3,4,4); plot(freq1.time(1:1000),pow_16Hz(1:1000));title('pow @ 16 Hz');ylim([0 0.4])
subplot(3,4,5); plot(freq1.time(1:1000),pow_18Hz(1:1000));title('pow @ 18 Hz');ylim([0 0.4])
subplot(3,4,6); plot(freq1.time(1:1000),pow_20Hz(1:1000));title('pow @ 20 Hz');ylim([0 0.4])
subplot(3,4,7); plot(freq1.time(1:1000),pow_22Hz(1:1000));title('pow @ 22 Hz');ylim([0 0.4])
subplot(3,4,8); plot(freq1.time(1:1000),pow_24Hz(1:1000));title('pow @ 24 Hz');ylim([0 0.4])
subplot(3,4,9); plot(freq1.time(1:1000),pow_26Hz(1:1000));title('pow @ 26 Hz');ylim([0 0.4])
subplot(3,4,10); plot(freq1.time(1:1000),pow_28Hz(1:1000));title('pow @ 28 Hz');ylim([0 0.4])
subplot(3,4,11); plot(freq1.time(1:1000),pow_30Hz(1:1000));title('pow @ 30 Hz');ylim([0 0.4])
print -dpng phalow_freqhigh_fig5.png

% mtmfft on freq1 with output power
cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'pow';
cfg.taper     = 'hanning';
cfg.foilim    = [2 60];
cfg.keeptrials = 'no';

freq2 = ft_freqanalysis(cfg,freq1); %FieldTrip automatically converts the freq1 data to raw data. 
                                 %Every frequency in the powerspectrum is converted to a a channel labeled mix@xxHz

freq2.freq2 = [2:2:60];

%plot
figure; imagesc(freq2.freq, freq2.freq2, freq2.powspctrm)
axis xy
print -dpng phalow_freqhigh_fig6.png


%zoom in
figure; imagesc(freq2.freq(1:33), freq2.freq2(5:15), freq2.powspctrm(5:15,1:33))
axis xy
print -dpng phalow_freqhigh_fig7.png

In figure 7 you can see that most frequencies are modulated at 2 Hz, which was indeed the frequency of the modulation (s2 in data). Harmonics are seen at 4 and 6 Hz. The base frequency (20 Hz, s1 in data) is the only frequency that shows the strongest correlation with the modulation frequency (s2.freq = 2Hz) * 2 = 4Hz.

% bandpass
cfg = [];
cfg.channel = 'mix';
cfg.bpfilter = 'yes';
cfg.bpfreq = [10 40];
data_bp = ft_preprocessing(cfg,data);

% check
figure;
sel = 1:1000;
plot(data.trial{1}(1,sel))
hold on 
plot(data_bp.trial{1}(1,sel),'r')
legend('original data','bandpassed data','location','Best')
print -dpng phalow_freqhigh_fig8.png

% estimate instantaneous phase by angle Hilbert
cfg = [];
cfg.channel = 'mix';
cfg.hilbert = 'angle';
data_H = ft_preprocessing(cfg,data_bp);

% check
figure;
sel = 1:1000;
plot(data.trial{1}(7,sel));
hold on 
plot(data_H.trial{1}(1,sel),'r')
legend(data.label{7},'estimated ins phase','location','Best')
print -dpng phalow_freqhigh_fig9.png

% calculate derivative of instantaneous phase
cfg = [];
cfg.channel = 'mix';
cfg.absdiff = 'yes'
data_diff = ft_preprocessing(cfg,data_H);

% check
figure;
sel = 1:3000;
plot(data_diff.trial{1}(1,sel),'r')
print -dpng phalow_freqhigh_fig10.png

cfg = [];
cfg.method    = 'mtmfft';
cfg.output    = 'pow';
cfg.taper     = 'hanning';
cfg.foilim    = [1 20];
cfg.keeptrials = 'no';

fft_data_diff = ft_freqanalysis(cfg, data_diff);
figure; ft_singleplotER([],fft_data_diff);
print -dpng phalow_freqhigh_fig11.png