contourf lower edge effects, how to fix and how does contourf work

Dear Matlab,

Using CONTOURF() to plot phase-amplitude coupling with an EEG signal, there appears what I call 'an edge effect' at the bottom of the plot (whereby the largest activity is always colored at the bottom of the plot) SEE ATTACHED IMAGE. I would like to capture the entire range of the 'hottest color blods' so I tried changing the y-axis range. BUT If I change the y-axis to range, no metter what - The 'hottest' colors always end up at the bottom of the plot. So I don't know how contourf() works then. Can you please help me? Thanks very much in advance!

Attached is sample data and pasted code below...

The variable that I change below to adjust/change the y-axis range is:

power_frex = linspace(15,60,45);

(and the resulting plot always has a blob at the bottom)

ATTACHED ARE EXAMPLE PLOTS FROM CHANGING THE Y-AXIS RANGE VIA THE 'power_frex' variable.

%% signal parameters

load signal

srate = 250;

lowerfreq = 6;
upperfreq = 40; % Hz

%%% Compute PAC

phase_frex = linspace(1,10,15);
power_frex = linspace(15,60,45);

pac = zeros(length(power_frex),length(phase_frex));

for lowfi=1:length(phase_frex)
    
    % get phase time series
    phasets = exp(1i*angle(hilbert(filterFGx(signal,srate,phase_frex(lowfi),1))));
    
    % loop over amplitude frequencies
    for hifi=1:length(power_frex)
        
        % get power time series
        powts = abs(hilbert(filterFGx(signal,srate,power_frex(hifi),lowerfreq*2))).^2;
        
        % PAC
        pac(hifi,lowfi) = abs(mean( powts.*phasets ));
    end
end


%%% visualize

figure(4), clf
contourf(phase_frex,power_frex,pac,40,'linecolor','none')
xlabel('Phase freq. (Hz)')
ylabel('Ampl. freq. (Hz)')

Joanne

Answer 1

[Star Strider]

We’re missing the filter function you’re using, so not able to run your code. However the pspectrum 'spectrogram' plot demonstrates that the pwer is concentrated at the lower frequencies (as expected). I’m not certain how you’re limiting the plots you’re creating, however setting the appropriate axis limits rather than restricting the values plotted may produce the result you want.

Try something like this —

LD = load(websave('TESTsignal','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1153968/TESTsignal.mat'));
signal = LD.signal;

srate = 250;
lowerfreq = 6;
upperfreq = 40; % Hz
%%% Compute PAC
phase_frex = linspace(1,10,15);
power_frex = linspace(15,60,45);

% figure
[p,f,t] = pspectrum(signal, srate, 'spectrogram');

figure
surfc(t,f,10*log10(p), 'EdgeColor','none')
colormap(turbo)
hcb = colorbar;
hcb.Label.String = 'Power (dB)';
zlim([-50 max(zlim)])
xlabel('Time')
ylabel('Frequency')
zlabel('Power (dB)')
title('Frequency: 0 - 125 Hz')

figure
surfc(t,f,10*log10(p), 'EdgeColor','none')
colormap(turbo)
hcb = colorbar;
hcb.Label.String = 'Power (dB)';
ylim([50 max(ylim)])
zlim([-50 max(zlim)])
xlabel('Time')
ylabel('Frequency')
zlabel('Power (dB)')
title('Frequency: 50 - 125 Hz')

.

[Joanne Hall]

Oh wait, ok sorry, here's the filter function to run the code (from Mike X Cohen), attached...

[Star Strider]

I set ‘power_frex’ to start at 1 here, then limited the axes (using ylim) to start at 15. The ‘lower edge effect’ appears to be gone with it.

Change the lower value in:

ylim([15 max(ylim)])

to be whatever you want to set the lower limit.

(For my own interest, I also plotted it as a surfc plot because I want to see wht the data are doing, and that tells me more than the contourf plot.)

Try this —

%% signal parameters

% load signal
LD = load(websave('TESTsignal','https://www.mathworks.com/matlabcentral/answers/uploaded_files/1153968/TESTsignal.mat'));
signal = LD.signal;

srate = 250;

lowerfreq = 6;
upperfreq = 40; % Hz

%%% Compute PAC

phase_frex = linspace(1,10,15);
power_frex = linspace(1,60,45);

pac = zeros(length(power_frex),length(phase_frex));

for lowfi=1:length(phase_frex)
    
    % get phase time series
    phasets = exp(1i*angle(hilbert(filterFGx(signal,srate,phase_frex(lowfi),1))));
    
    % loop over amplitude frequencies
    for hifi=1:length(power_frex)
        
        % get power time series
        powts = abs(hilbert(filterFGx(signal,srate,power_frex(hifi),lowerfreq*2))).^2;
        
        % PAC
        pac(hifi,lowfi) = abs(mean( powts.*phasets ));
    end
end


%%% visualize

figure(4), clf
contourf(phase_frex,power_frex,pac,40,'linecolor','none')
xlabel('Phase freq. (Hz)')
ylabel('Ampl. freq. (Hz)')
ylim([15 max(ylim)])

figure
hsc = surfc(phase_frex,power_frex,pac, 'EdgeColor','none');
xlabel('Phase freq. (Hz)')
ylabel('Ampl. freq. (Hz)')
title('Surface Plot')
% ylim([15 max(ylim)])
hsc(2).LevelList = linspace(min(pac(:)), max(pac(:)), 40);
title('Original Y-Limts (1 - 60)')

figure
hsc = surfc(phase_frex,power_frex,pac, 'EdgeColor','none');
xlabel('Phase freq. (Hz)')
ylabel('Ampl. freq. (Hz)')
title('Surface Plot')
ylim([15 max(ylim)])
hsc(2).LevelList = linspace(min(pac(:)), max(pac(:)), 40);
title('Reduced Y-Limts (15 - 60)')

function [filtdat,empVals,fx] = filterFGx(data,srate,f,fwhm,showplot)
% filterFGx   Narrow-band filter via frequency-domain Gaussian
%  [filtdat,empVals] = filterFGx(data,srate,f,fwhm,showplot)
% 
% 
%    INPUTS
%       data : 1 X time or chans X time
%      srate : sampling rate in Hz
%          f : peak frequency of filter
%       fhwm : standard deviation of filter, 
%              defined as full-width at half-maximum in Hz
%   showplot : set to true to show the frequency-domain filter shape
% 
%    OUTPUTS
%    filtdat : filtered data
%    empVals : the empirical frequency and FWHM (in Hz and in ms)
% 
% Empirical frequency and FWHM depend on the sampling rate and the
% number of time points, and may thus be slightly different from
% the requested values.
% 
% [e-mail]

%% input check

if size(data,1)>size(data,2)
%     help filterFGx
%     error('Check data size')
end

if (f-fwhm)<0
%     help filterFGx
%     error('increase frequency or decrease FWHM')
end

if nargin<4
    help filterFGx
    error('Not enough inputs')
end

if fwhm<=0
    error('FWHM must be greater than 0')
end

if nargin<5
    showplot=false;
end

%% compute and apply filter

% frequencies
hz = linspace(0,srate,size(data,2));

% create Gaussian
s  = fwhm*(2*pi-1)/(4*pi); % normalized width
x  = hz-f;                 % shifted frequencies
fx = exp(-.5*(x/s).^2);    % gaussian
fx = fx./abs(max(fx));     % gain-normalized

%% filter

% filtdat = 2*real( ifft( bsxfun(@times,fft(data,[],2),fx) ,[],2) );
filtdat = 2*real( ifft( fft(data',[],1).*fx') )';

%% compute empirical frequency and standard deviation

idx = dsearchn(hz',f);
empVals(1) = hz(idx);

% find values closest to .5 after MINUS before the peak
empVals(2) = hz(idx-1+dsearchn(fx(idx:end)',.5)) - hz(dsearchn(fx(1:idx)',.5));

% also temporal FWHM
tmp = abs(hilbert(real(fftshift(ifft(fx)))));
tmp = tmp./max(tmp);
tx = (0:length(data)-1)/srate;
[~,idxt] = max(tmp);
empVals(3) = (tx(idxt-1+dsearchn(tmp(idxt:end)',.5)) - tx(dsearchn(tmp(1:idxt)',.5)))*1000;

%% inspect the Gaussian (turned off by default)

if showplot
    figure(10001+showplot),clf
    subplot(211)
    plot(hz,fx,'o-')
    hold on
    plot([hz(dsearchn(fx(1:idx)',.5)) hz(idx-1+dsearchn(fx(idx:end)',.5))],[fx(dsearchn(fx(1:idx)',.5)) fx(idx-1+dsearchn(fx(idx:end)',.5))],'k--')
    set(gca,'xlim',[max(f-10,0) f+10]);
    
    title([ 'Requested: ' num2str(f) ', ' num2str(fwhm) ' Hz; Empirical: ' num2str(empVals(1)) ', ' num2str(empVals(2)) ' Hz' ])
    xlabel('Frequency (Hz)'), ylabel('Amplitude gain')
    
    subplot(212)
    tmp1 = real(fftshift(ifft(fx))); tmp1 = tmp1./max(tmp1);
    tmp2 = abs(hilbert(tmp1));
    plot(tx,tmp1, tx,tmp2), zoom on
    xlabel('Time (s)'), ylabel('Amplitude gain')
end

%% done.
end

Also, to be used in this context, the ‘filterFGx’ function needs an end statement. I provided it here.

.

[Joanne Hall]

Thanks so much Star!!!

[Star Strider]

As always, my pleasure!