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!
채택된 답변
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
Angela Yang
2023년 4월 12일
Mathieu NOE
2023년 4월 12일
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
2023년 4월 12일
Mathieu NOE
2023년 4월 12일
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
2023년 4월 12일
Mathieu NOE
2023년 4월 13일
As always, my pleasure !
추가 답변 (0개)
카테고리
도움말 센터 및 File Exchange에서 Descriptive Statistics에 대해 자세히 알아보기
제품
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
