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

Question

Joanne Hall 2022-10-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1824003-contourf-lower-edge-effects-how-to-fix-and-how-does-contourf-work

评论： Star Strider 2022-10-15

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2022-10-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1824003-contourf-lower-edge-effects-how-to-fix-and-how-does-contourf-work#answer_1072843

在 MATLAB Online 中打开

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')

.

4 个评论
显示 2更早的评论隐藏 2更早的评论

Star Strider 2022-10-12

在 MATLAB Online 中打开

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

%%% 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.

%

% mikexcohen@gmail.com

%% 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 2022-10-15

Thanks so much Star!!!

Star Strider 2022-10-15

As always, my pleasure!

请先登录，再进行评论。

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论
显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

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

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

4 个评论 显示 2更早的评论隐藏 2更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论