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!