I'd appreciate some help from you guys.
I'm currently analyzing an experiment where I've entrained participants to these frequencies:
Condition 1: 4.16Hz-8.33Hz
Condition 2: 2.77Hz-8.33Hz
Condition 3: 2.4Hz
For each condition I passed 8 trials of 26 seconds at Fs = 500Hz. I'm analysing the data using FieldTrip and some toolboxes in Matlab. I calculated the Power and PSD using pwelch function after processing the data in Fieldtrip.
I want my data in dB units and get rid of the $1/f^\gamma$ trend where $0< \gamma<2$ (Demanuele et al., 2007)
I have 2 problems:
1) I don't know how to remove the 1/f trend correctly since I don't understand how the article that I mention does it (see below) in the frequency domain (I have tried fitting $1/f^\gamma$ with best $\gamma$, using the lsqcurvefit function. Using the medfilt1 with order 15 gives me the best results. Though I think there should be a better way, also that can be defended on an article).
"the spectrogram of one EEG channel was calculated. The median across all time windows was then found for every frequency point. Thus, a graph of the median PSD value for every frequency was obtained. The same procedure was repeated for all the EEG channels. Then, the overall median of the median PSD curves of all the channels was calculated. The same was done for each of the 'task' and 'rest' condition. The average was then calculated across all conditions and this was used as our normalisation curve."
EDIT I thought it might be helpful to provide the code I made. Sorry for the inefficiencies you may find since I am no developer:
% Choose subject
load matfiles/Subject05
% loading subject5 where (channel,samples,trial,condition) is a 4D matrix
samples = numel(subject5(1,:,1,1)); % number of samples
% Number of samples per window, the less windows the better frequency resolution and worst time resolution. I tried different windows.
nsc = floor(samples/4.5);
nov = floor(nsc*0.75); % overlap of 75%
fr = 0.5:0.05:15; % foi
fs = 500; % sampling freq
%initialize variable
medps = zeros(291,28,21,3);
% loop conditions
for cond = 1 : length(subject5(1,1,1,:))
%loop trials
for trial = 1 : length(subject5(1,1,:,cond))
% loop channels
for chan = 1 : length(subject5(:,1,trial,cond))
% calculate spectrogram for each channel, trial and condition
%gives t(time) and ps (PSD estimate for each freq and time in Amplitude^2/Hz units)
[~,~,t,ps] = spectrogram(subject5(chan,:,trial,cond),hamming(nsc),nov,fr,fs);
% calculate median over time points per channel, trial and condition
medps(:,chan,trial,cond) = median(ps,2);
end
end
end
%Calculate median over channels per trial and condition
medchan = median(medps,2); %this produces matrix freqx1channelxtrialsxcond
medchan = permute(medchan,[1 3 4 2]); %get 4D into 3D
%Calculate median over trials per condition
medtrial = median(medchan,2);
medtrial = permute(medtrial,[1 3 2]); %Get a 2D matrix freqxcond
%Mean over conditions is our Normalization Curve
normalizationCurve = mean(medtrial,2);
2) Once I transform the data to dB units (i.e. $10*log10$) it gets really noisy, making the peaks at the foi (frequencies of interest) as big as the background or almost. I've looked into SNR but really don't know how to apply it here. The only thing I've done and still doesn't remove all the noise is the medfilt1 mentioned in 1).
I hope any of you with experience in this can help me solve this issue!