Shading in the SEM region surrounding a line plot

Hi! I'm trying to shade in the SEM of 98 participants in my plot but I can't get it to work.. I just can't manage to get the SEM shading to show up on the figure.

I've included a figure of what I am trying to recreate.

Here's my code for simulating the data, I've only included lines for the control group:

% simulating control EEG
srate=2*100;
t=0:1/srate:500; 

% simulating grand average
% ------- part 1 theta + alpha wave parameters (3-10Hz)
theta_hz=3:10;
theta_amp=linspace(75,168.5,length(theta_hz)); % normalised -1:1 translated into 0:200
omega=2*pi*theta_hz; % radial frequency

thetawave=zeros(length(theta_hz),length(t)); % initialising
for i=1:length(theta_hz)
    % generate theta wave
    con_thetawave(i,:)=theta_amp(i)*sin(omega(i)*t); 
end

% ------- part 2 alpha wave (11-14Hz)
p2_hz=11:14;
p2_amp=linspace(160,107.5,length(p2_hz));
omega=2*pi*p2_hz; % radial frequency

p2wave=zeros(length(p2_hz),length(t)); % initialising
for i=1:length(p2_hz)
    con_p2wave(i,:)=p2_amp(i)*sin(omega(i)*t); 
end 

% ------- beta wave (15-30Hz)
beta_hz=15:30;
beta_amp=[110 linspace(112.5,115,4) linspace(115,105,5) 103 linspace(101,98,4) 100];
omega=2*pi*beta_hz; % radial frequency

betawave=zeros(length(beta_hz),length(t)); % initialising
for i=1:length(beta_hz)
    con_betawave(i,:)=beta_amp(i)*sin(omega(i)*t); 
end 

% ------- gamma wave (30-50Hz)
gamma_hz=31:50;
gamma_amp=[linspace(95,90,6) linspace(90,93,5) linspace(93,97.5,4) 96.5 linspace(96,100,3) 99];
omega=2*pi*gamma_hz; % radial frequency

gammawave=zeros(length(gamma_hz),length(t)); % initialising
for i=1:length(gamma_hz)
    con_gammawave(i,:)=gamma_amp(i)*sin(omega(i)*t); 
end         

% combine all frequencies
con_avgEEG=sum(con_thetawave)+sum(con_p2wave)+sum(con_betawave)+sum(con_gammawave);

% fourier transform
con_avgEEGX=fftshift(fft(con_avgEEG));
lx=length(con_avgEEGX);
hz=linspace(0,srate/2,floor(lx/2)+1); % scaling frequency axis into Hz
con_yX=2*abs(con_avgEEGX(ceil(lx/2):end));
con_norm_yX=normalize(con_yX,'range',[min(con_yX) max(con_yX)]); % normalise values to [-1 1]

% fitting line to envelope fourier transform
[con_locs,con_pks]=peakseek(con_norm_yX);
con_xpeaks=hz(con_locs);

% visualisation
figure(5), subplot(221)
plot(hz,con_norm_yX,con_xpeaks,con_pks,'b')
xlim([min(theta_hz) max(gamma_hz)])
xlabel('Frequency (Hz)'), ylabel('Normalised power (AU)')
title('Frequency Peaks With Envelope')

My code for generating the SEM

% control: N=98 | schizophrenia: N=95 (sample size same as Adams, Pinotsis et al. (2022)
conN=98; % sample sizes
sczN=95;

% control group
% adding random noise for each participant
allhzs=[theta_hz p2_hz beta_hz gamma_hz];
allamps=[theta_amp p2_amp beta_amp gamma_amp]; %!! remove
allwavs=[con_thetawave;con_p2wave;con_betawave;con_gammawave]; % combine all 'averaged' frequency waves into one variable

% normally distributed noise amplitude
mu=.03; % mean
sigma=.5; % standard deviation
con_noiseamp=sigma.*randn(1,conN)+mu; % generating normal distribution

rng(69) % set seed
con_noisyEEG=zeros(conN,length(con_pks)); % intialising group variable
%noisyEEG=zeros(length(allhzs),length(t)); % intialising individual variable
for i=1:conN
    noise=randn(size(con_pks))*con_noiseamp(i); % creating noise
    con_noisyEEG(i,:)=con_pks+noise;

end 
% calculating SEM
con_sem=(std(con_noisyEEG)/sqrt(conN)); 

x_fill=[hz(con_locs),fliplr(hz(con_locs))];
y_upper=con_locs+con_sem;
y_lower=con_locs-con_sem;
y_fill=[y_upper,fliplr(y_lower)];

% scaling
scale_x_fill=rescale(x_fill,min(con_pks),max(con_pks));
scale_y_fill=rescale(y_fill,min(con_pks),max(con_pks));

% visualisation
subplot(222) % spectra envelope only
plot(con_xpeaks,con_pks,'b'), hold on
xline([7 14 31],':')

patch(x_fill,y_fill,'g')

xlim([min(theta_hz) max(gamma_hz)])
xlabel('Frequency (Hz)'), ylabel('Normalised power (AU)')
title('Spectra Envelope')

Answer 1

[Mathieu NOE]

hello Angela and welcome back

here is my suggestion based on surface

seems also that there is no grid in your reference picture so I removed it (grid off)

this is the visual result

% simulating control EEG
srate=2*100;
t=0:1/srate:500; 
% simulating grand average
% ------- part 1 theta + alpha wave parameters (3-10Hz)
theta_hz=3:10;
theta_amp=linspace(75,168.5,length(theta_hz)); % normalised -1:1 translated into 0:200
omega=2*pi*theta_hz; % radial frequency
thetawave=zeros(length(theta_hz),length(t)); % initialising
for i=1:length(theta_hz)
    % generate theta wave
    con_thetawave(i,:)=theta_amp(i)*sin(omega(i)*t); 
end
% ------- part 2 alpha wave (11-14Hz)
p2_hz=11:14;
p2_amp=linspace(160,107.5,length(p2_hz));
omega=2*pi*p2_hz; % radial frequency
p2wave=zeros(length(p2_hz),length(t)); % initialising
for i=1:length(p2_hz)
    con_p2wave(i,:)=p2_amp(i)*sin(omega(i)*t); 
end 
% ------- beta wave (15-30Hz)
beta_hz=15:30;
beta_amp=[110 linspace(112.5,115,4) linspace(115,105,5) 103 linspace(101,98,4) 100];
omega=2*pi*beta_hz; % radial frequency
betawave=zeros(length(beta_hz),length(t)); % initialising
for i=1:length(beta_hz)
    con_betawave(i,:)=beta_amp(i)*sin(omega(i)*t); 
end 
% ------- gamma wave (30-50Hz)
gamma_hz=31:50;
gamma_amp=[linspace(95,90,6) linspace(90,93,5) linspace(93,97.5,4) 96.5 linspace(96,100,3) 99];
omega=2*pi*gamma_hz; % radial frequency
gammawave=zeros(length(gamma_hz),length(t)); % initialising
for i=1:length(gamma_hz)
    con_gammawave(i,:)=gamma_amp(i)*sin(omega(i)*t); 
end         
% combine all frequencies
con_avgEEG=sum(con_thetawave)+sum(con_p2wave)+sum(con_betawave)+sum(con_gammawave);
% fourier transform
con_avgEEGX=fftshift(fft(con_avgEEG));
lx=length(con_avgEEGX);
hz=linspace(0,srate/2,floor(lx/2)+1); % scaling frequency axis into Hz
con_yX=2*abs(con_avgEEGX(ceil(lx/2):end));
con_norm_yX=normalize(con_yX,'range',[min(con_yX) max(con_yX)]); % normalise values to [-1 1]
% fitting line to envelope fourier transform
[con_locs,con_pks]=peakseek(con_norm_yX);
con_xpeaks=hz(con_locs);
% visualisation
figure(5), subplot(221)
plot(hz,con_norm_yX,con_xpeaks,con_pks,'b')
xlim([min(theta_hz) max(gamma_hz)])
xlabel('Frequency (Hz)'), ylabel('Normalised power (AU)')
title('Frequency Peaks With Envelope')
%%My code for generating the SEM
% control: N=98 | schizophrenia: N=95 (sample size same as Adams, Pinotsis et al. (2022)
conN=98; % sample sizes
sczN=95;
% control group
% adding random noise for each participant
allhzs=[theta_hz p2_hz beta_hz gamma_hz];
allamps=[theta_amp p2_amp beta_amp gamma_amp]; %!! remove
allwavs=[con_thetawave;con_p2wave;con_betawave;con_gammawave]; % combine all 'averaged' frequency waves into one variable
% normally distributed noise amplitude
mu=.03; % mean
sigma=.5; % standard deviation
con_noiseamp=sigma.*randn(1,conN)+mu; % generating normal distribution
rng(69) % set seed
con_noisyEEG=zeros(conN,length(con_pks)); % intialising group variable
%noisyEEG=zeros(length(allhzs),length(t)); % intialising individual variable
for i=1:conN
    noise=randn(size(con_pks))*con_noiseamp(i); % creating noise
    con_noisyEEG(i,:)=con_pks+noise;
end 
% calculating SEM
con_sem=(std(con_noisyEEG)/sqrt(conN)); 
y_upper=con_locs+con_sem;
y_lower=con_locs-con_sem;
% visualisation
subplot(222) % spectra envelope only
plot(con_xpeaks,con_pks,'b'), hold on, grid off
xline([7 14 31],':')


    % shading around the line made with surface 
    x = con_xpeaks;
    y = con_pks;
    z = zeros(size(x));
    surface([x;x],[y;y],[z;z], 'linewidth',7, 'edgecolor','g', 'edgealpha',0.2);
    hold off;
    
xlim([min(theta_hz) max(gamma_hz)])
xlabel('Frequency (Hz)'), ylabel('Normalised power (AU)')
title('Spectra Envelope')

Unrecognized function or variable 'peakseek'.

[Mathieu NOE]

hello Angela

did my suggestion fullfil your request ?

[Angela Yang]

completly forgot to answer but yes Mathieu, super helpful as always!

[Mathieu NOE]

thanks !

glad I could help