How to fit a line over a plot
Show older comments
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!
Accepted Answer
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, <WangLeiBox@hotmail.com>, 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
6 Comments
Angela Yang
on 12 Apr 2023
Mathieu NOE
on 12 Apr 2023
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
on 12 Apr 2023
Mathieu NOE
on 12 Apr 2023
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
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
on 12 Apr 2023
Mathieu NOE
on 13 Apr 2023
As always, my pleasure !
More Answers (0)
Categories
Find more on Measurements and Feature Extraction in Help Center and File Exchange
Products
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
