How to denoise and get rid of 1/f trend of a PSD in an EEG synchronization experiment 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&lt; \gamma&lt;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!
Apparently the paper you mention is simply evaluating the median in the EEG spetrum. Say that for all your channels you get a 3D matrix PSD, with dimension 1 the number of channel, dimension 2 the time windows, and dimension 3 the frequencies. In pseudocode the article computes M = median(median(PSD,2),1), then removes the trend with CorrectedPSD = PSD - M. Which is actually fairly simple.