How to fit a line over a plot

Hi, I am trying to recreate this figure but I can't manage to figure out how to fit a smoothed line to summarise the frequencies. Also, the reference figure also includes the SEM which I'm struggling to recreate. My current plot is the grand averaged frequency spectra of simulated EEG data from 98 participants. Thanks!

Answer 1

[Mathieu NOE]

hello

seems to me you want to plot the envelope of your spectra

there is a matlab function (Signal Processing Toolbox required) for that (envelope.m)

If you don't have the toolbox , here's an alternative (adapt to your own data) :

t = 0:0.01:10;
x = exp(-(t-5).^2/2).*sin(2*pi*5*t);

% obtain the envelope data
%--------------------------------------------
y = abs(hilbert(x)); % option 1 with hilbert
[up,down] = envelope2(t,x,'linear'); % option 2 with envelope2 (see below)

plot(t,x,t,y,'k',t,up,'m',t,down,'c')
legend('signal','hilbert','envelope2 / upper env','envelope2 / lower env');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [up,down] = envelope2(x,y,interpMethod)

%ENVELOPE gets the data of upper and down envelope of the known input (x,y).
%
%   Input parameters:
%    x               the abscissa of the given data
%    y               the ordinate of the given data
%    interpMethod    the interpolation method
%
%   Output parameters:
%    up      the upper envelope, which has the same length as x.
%    down    the down envelope, which has the same length as x.
%
%   See also DIFF INTERP1

%   Designed by: Lei Wang, <[e-mail]>, 11-Mar-2003.
%   Last Revision: 21-Mar-2003.
%   Dept. Mechanical & Aerospace Engineering, NC State University.
% $Revision: 1.1 $  $Date: 3/21/2003 10:33 AM $

if length(x) ~= length(y)
    error('Two input data should have the same length.');
end

if (nargin < 2)|(nargin > 3),
 error('Please see help for INPUT DATA.');
elseif (nargin == 2)
    interpMethod = 'linear';
end

% Find the extreme maxim values 
% and the corresponding indexes
%----------------------------------------------------
extrMaxIndex = find(diff(sign(diff(y)))==-2)+1;
extrMaxValue = y(extrMaxIndex);

% Find the extreme minim values 
% and the corresponding indexes
%----------------------------------------------------
extrMinIndex = find(diff(sign(diff(y)))==+2)+1;
extrMinValue = y(extrMinIndex);

up = extrMaxValue;
up_x = x(extrMaxIndex);

down = extrMinValue;
down_x = x(extrMinIndex);

% Interpolation of the upper/down envelope data
%----------------------------------------------------
up = interp1(up_x,up,x,interpMethod); 
down = interp1(down_x,down,x,interpMethod); 
end

[Angela Yang]

Hi thanks for the reply! If I'm not wrong, this is the envelope in the time domain. When I try this function and plot it in the frequency domain it doesn't outline the frequency spectra, instead it hugs each peak if that makes sense..

[Mathieu NOE]

hello again

yes my demo can be interpreted as for a time domain signal but the envelope function does not care if it time or frequency domain... or whatever , it's data

If you're in trouble can you share your data (and code if possible ) ?

[Angela Yang]

Hi! Here's the code containing the simulated data for a single individual. Also, I think there might be a way to plot the outline of the frequency spectra without the envelope function because I saw some other student make a similar graph without it but I just can't figure out the code..

% sample size
conN=98;
sczN=95;

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

% ------- 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
    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));
for i=1:length(p2_hz)
    omega=2*pi*p2_hz; % radial frequency
    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];
for i=1:length(beta_hz)
    omega=2*pi*beta_hz; % radial frequency
    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];
for i=1:length(gamma_hz)
    omega=2*pi*gamma_hz; % radial frequency
    gammawave(i,:)=gamma_amp(i)*sin(omega(i)*t); 
end         

allhzs=[theta_hz,p2_hz,beta_hz,gamma_hz];
allamps=[theta_amp p2_amp beta_amp gamma_amp];
individEEG=sum(thetawave)+sum(p2wave)+sum(betawave)+sum(gammawave); % combine all frequencies

% fourier transform
individEEGX=fftshift(fft(individEEG));
lx=length(individEEGX);
hz=linspace(0,srate/2,floor(lx/2)+1); % scaling frequency axis into Hz

% fitting line
figure, subplot(221)
plot(t,individEEG)
xlim([0 50])

subplot(223)
plot(hz,2*abs(individEEGX(ceil(lx/2):end)))
xlim([min(allhzs) max(allhzs)])

subplot(224)
plot(allhzs,allamps)
xlim([min(allhzs) max(allhzs)])

[Mathieu NOE]

I assume you wanted to add a line (red here) on top of the fft peaks to create that envelope.

You can do it simply with findpeaks or with the simpler (and faster) peakseek function provided at the end of the code

(available on the Fex here : PeakSeek - File Exchange - MATLAB Central (mathworks.com) )

hope it helps

% sample size
conN=98;
sczN=95;
% simulating control EEG
srate=2*100;
t=0:1/srate:500; 
% ------- 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
    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));
for i=1:length(p2_hz)
    omega=2*pi*p2_hz; % radial frequency
    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];
for i=1:length(beta_hz)
    omega=2*pi*beta_hz; % radial frequency
    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];
for i=1:length(gamma_hz)
    omega=2*pi*gamma_hz; % radial frequency
    gammawave(i,:)=gamma_amp(i)*sin(omega(i)*t); 
end         
allhzs=[theta_hz,p2_hz,beta_hz,gamma_hz];
allamps=[theta_amp p2_amp beta_amp gamma_amp];
individEEG=sum(thetawave)+sum(p2wave)+sum(betawave)+sum(gammawave); % combine all frequencies
% fourier transform
individEEGX=fftshift(fft(individEEG));
lx=length(individEEGX);
hz=linspace(0,srate/2,floor(lx/2)+1); % scaling frequency axis into Hz
% fitting line
figure, subplot(221)
plot(t,individEEG)
xlim([0 50])
subplot(223)
y = 2*abs(individEEGX(ceil(lx/2):end));
% find peaks (with findpeaks or with peakseek (simpler and faster)) and
% draw a "top line" on top the fft lines
[locs, pks]=peakseek(y);
xpeaks = hz(locs);
plot(hz,y,xpeaks,pks,'r')
xlim([min(allhzs) max(allhzs)])
subplot(224)
plot(allhzs,allamps)
xlim([min(allhzs) max(allhzs)])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [locs, pks]=peakseek(x,minpeakdist,minpeakh)
% x is a vector input (generally a timecourse)
% minpeakdist is the minimum desired distance between peaks (optional, defaults to 1)
% minpeakh is the minimum height of a peak (optional)
%
% (c) 2010
% Peter O'Connor
% peter<dot>ed<dot>oconnor .AT. gmail<dot>com
if size(x,2)==1, x=x'; end
% Find all maxima and ties
locs=find(x(2:end-1)>=x(1:end-2) & x(2:end-1)>=x(3:end))+1;
if nargin<2, minpeakdist=1; end % If no minpeakdist specified, default to 1.
if nargin>2 % If there's a minpeakheight
    locs(x(locs)<=minpeakh)=[];
end
if minpeakdist>1
    while 1
        del=diff(locs)<minpeakdist;
        if ~any(del), break; end
        pks=x(locs);
        [garb, mins]=min([pks(del) ; pks([false del])]); %#ok<ASGLU>
        deln=find(del);
        deln=[deln(mins==1) deln(mins==2)+1];
        locs(deln)=[];
    end
end
if nargout>1
    pks=x(locs);
end
end

[Angela Yang]

Ahh yes this is exactly what I'm looking for! Bless your soul, thank you so much! :)

[Mathieu NOE]

As always, my pleasure !