text:
function err = alignmentError(pcCalib, pcGroundTruth, param, camW, camH, verbose, figName)
%% Calculate alignment error between calibration and Intel RealSense F200 captured point cloud.
% Calculate and compare the 3D alignment errors between reconstructed
% point cloud and Intel RealSense F200 captured point cloud. Four
% calibration methods are compared.

%% License
% ACADEMIC OR NON-PROFIT ORGANIZATION NONCOMMERCIAL RESEARCH USE ONLY
% Copyright (c) 2018 Bingyao Huang
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:

% The above copyright notice and this permission notice shall be included in all
% copies or substantial portions of the Software.

% If you publish results obtained using this software, please cite our paper.

% THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
% IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
% FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
% AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
% LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
% OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
% SOFTWARE.

%%
if(nargin < 7)
    figName = '';
end

%%
% convert point cloud to XYZ map
% imPC = zeros(camH, camW, 3);
% imPC(:,:,1) = (reshape(pcGroundTruth.Location(:,1), camW, camH))';
% imPC(:,:,2) = (reshape(pcGroundTruth.Location(:,2), camW, camH))';
% imPC(:,:,3) = (reshape(pcGroundTruth.Location(:,3), camW, camH))';
% imPC = double(imPC);

% pts3dGT = double(pcGroundTruth.Location);
% pts3dGT = pts3dGT(pts3dGT(:,3)>0,:);
% pts2dGT = cv.projectPoints(pts3dGT, [0,0,0], [0,0,0], param.camK, 'DistCoeffs', param.camKc);
% imD = zeros(camH, camW);
% pts2dGT = round(pts2dGT);
% inlierIdx = pts2dGT(:,1)<=camW & pts2dGT(:,2)<=camH & pts2dGT(:,1)>0 & pts2dGT(:,2)>0;
% pts2dGT = pts2dGT(inlierIdx,:);
% pts3dGT = pts3dGT(inlierIdx,:);
% imD(sub2ind([camH, camW], pts2dGT(:,2), pts2dGT(:,1))) = pts3dGT(:,3);

pts3d = double(pcCalib.Location);
pts3d = pts3d(pts3d(:,3)>0,:);
pts2d = cv.projectPoints(pts3d, [0,0,0], [0,0,0], param.camK, 'DistCoeffs', param.camKc);
imD = zeros(camH, camW);
pts2d = round(pts2d);
inlierIdx = pts2d(:,1)<=camW & pts2d(:,2)<=camH & pts2d(:,1)>0 & pts2d(:,2)>0;
pts2d = pts2d(inlierIdx,:);
pts3d = pts3d(inlierIdx,:);
imD(sub2ind([camH, camW], pts2d(:,2), pts2d(:,1))) = pts3d(:,3);

% find depth image mask
imMask = false(size(imD));
imMask(imD > 0) = 1;
se = strel('disk', 5);
% imMask = imclose(imopen(imMask,  se), se);
imMask = imclose(imMask, se);
imMask = imfill(imMask, 'holes');
imMask = bwareafilt(imMask, 1);
imMask = imgaussfilt(double(imMask), 7) > 0.2; % smooth edges
% fs(imMask);

%% align reconstructed point cloud with ground truth point cloud using icp

% downsample and denoise to save time and memory
% gridSize = 0.1;
% fixed = pcdownsample(pcdenoise(pcGroundTruth, 'Threshold',5), 'gridAverage', gridSize);
% % moving = pcdenoise(pcdownsample(pointCloud(pts3dInterp), 'gridAverage', gridSize));
% moving = pcdownsample(pcdenoise(pcCalib,'Threshold',5), 'gridAverage', gridSize);
%
% % register pcCalib to fixed pcGroundTruth
% [tform, ~, ~] = pcregistericp(moving, fixed,'Extrapolate', true, 'MaxIterations', 50);
% pcCalibAligned = pctransform(pcCalib, tform);
% figure;pcshowpair(fixed, pcCalibAligned)

%% interpolate point cloud
rp = regionprops(imMask, 'boundingbox');
desity = 1;
npX = rp.BoundingBox(3)*desity;
npY = rp.BoundingBox(4)*desity;

% xRange contains lower and higher bounds
xRange = double(pcCalib.XLimits);
xDiff = xRange(1) - xRange(2);
xStep = -xDiff / npX;

yRange = double(pcCalib.YLimits);
yDiff = yRange(1) - yRange(2);
yStep = -yDiff/npY;

% use griddata to interpolate a mesh from given point cloud
% pts3d = pcCalibAligned.Location;
[xq,yq] = meshgrid(xRange(1):xStep:xRange(2), yRange(1):yStep:yRange(2));
zq = griddata(pts3d(:,1),pts3d(:,2),pts3d(:,3), xq, yq,'cubic');

% pcInterp is the interpolated denser point cloud
pcInterp = pointCloud([xq(:), yq(:), zq(:)]);

% finer registration
% moving = pcdownsample(pcdenoise(pcInterpAligned,'Threshold',5), 'gridAverage', gridSize);
% [tform, ~, ~] = pcregrigid(moving, fixed,'Extrapolate', true);
% pcInterpAligned = pctransform(pcInterpAligned, tform);

[~, ~, ~, pcInlierIdx] = Reconstruct.filterPointCloud(pcInterp, imMask, param, camW, camH);
pcOutlierIdx = setdiff(1:size(xq,1)*size(xq,2), pcInlierIdx);

if(verbose)
    figure('Name', figName);
    h = surf(xq, yq, zq, 'FaceColor', 'interp', 'FaceLighting', 'gouraud'); % square mesh
%     h = trisurf(delaunay(xq,yq),xq, yq, zq, 'FaceColor', 'interp', 'FaceLighting', 'gouraud'); % triangular mesh

    % remove interpolated meth that are out of ground truth range
    h.XData(pcOutlierIdx) = nan;
    h.YData(pcOutlierIdx) = nan;
    h.ZData(pcOutlierIdx) = nan;
    h.CData(pcOutlierIdx) = nan;

    h.EdgeColor = 'none';
%     h.EdgeColor = [0.3,0.3,0.3];
    h.LineStyle = ':';
    colormap parula

    % title
    title(figName)
    caxis([0 50])
%     colorbar

    % remove ticks
    set(gca,'xtick',[])
    set(gca,'xticklabel',[])
    set(gca,'ytick',[])
    set(gca,'yticklabel',[])
    set(gca,'ztick',[])
    set(gca,'zticklabel',[])

    % font
    set(gca, 'FontSize', 20);
    daspect([1 1 1]);
    axis vis3d tight

    rotate3d on
    box on
    set(gca, 'visible', 'off'); % ge trid of ugly box edges
end

light('Position',[0 1 -1]);
material dull % metal, shiny, default

%% calculte 3D alignment error using knn search
[~, err] = knnsearch(pcGroundTruth.Location, pcInterp.Location);

% visualize errors as pseudocolor
if(verbose)
    errOutlierIdx = isnan(err);
    errInlierIdx = ~errOutlierIdx;
    h.ZData(errOutlierIdx) = nan;
    h.CData(errInlierIdx) = err(errInlierIdx); % use err as color
%      h.CData(errInlierIdx) = abs(err(errInlierIdx)); % use abs err as color
end

% remove nan values
err = err(~isnan(err));
end

author: BingyaoHuang

repo: single-shot-pro-cam-calib


text:
function value = r8vec_in_ab ( n, x, a, b )

%*****************************************************************************80
%
%% R8VEC_IN_AB is TRUE if the entries of an R8VEC are in the range [A,B].
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    26 March 2012
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer N, the number of entries.
%
%    Input, real X(N), the vector.
%
%    Input, real A, B, the limits, with A <= B.
%
%    Output, logical VALUE, is TRUE if every entry is
%    between A and B.
%
  if ( any ( x(1:n) < a | b < x(1:n) ) )
    value = 0;
  else
    value = 1;
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function prepareDeepPerFrameFtrsCompiled(moviefilename,trackfilename,stationary,MCRloc,tempdir,maxcores)

% function prepareDeepPerFrameFtrsCompiled(moviefilename,trackfilename,stationary,maxcores)
setup;
usedeep = true;
method = 'deep-sup';
flowname = 'DS';

if ~exist(maxcores,'var')
  totcores = features('numcores');
  if totcores>10,
    maxcores = totcores-4;
  else
    maxcores = totcores-1;
  end
end
gparams = getParams;
[readfcn,nframes,fid,headerinfo] = get_readframe_fcn(moviename);

nblocks = ceil(nframes/gparams.blocksize);
expdir = fileparts(moviefilename);
[~,expname] = fileparts(expdir);
tempname = fullfile(tempdir,expname);
partcmd = sprintf('%s %s %s %s %d deep-sup ',...
  'computeFeaturesCompiled/for_testing/run_computeFeaturesCompiled.sh',...
  MCRloc,moviefilename,trackfilename,stationary);

author: kristinbranson

repo: JAABA


text:
function rr = PZX(rel_data, wins, up)
%PZX estimates the RR using a combination of zero-crossing detection, and
% peak detection.
%	            PZX(option, up)
%
%	Inputs:
%       rel_data    .t  vector of times
%                   .v  vector of resp Sig values
%                   .fs sampling freq
%       up              universal parameters structure
%       wins        .t  vector of start times
%
%	Outputs:
%       rr          .t  vector of times of estimated RRs
%                   .v  vector of estimated RRs
%

%% Cycle through windows
rr.t = mean([wins.t_start(:)' ; wins.t_end(:)']); rr.t = rr.t(:);
rr.v = nan(length(rr.t),1);

for win_no = 1 : length(wins.t_start)

    % extract relevant data
    rel_els = find(rel_data.t >= wins.t_start(win_no) & rel_data.t < wins.t_end(win_no));
    data.v = rel_data.v(rel_els);
    data.t = rel_data.t(rel_els);

    good_els = ~isnan(data.v);
    data.v = data.v(good_els); data.v = data.v(:);
    data.t = data.t(good_els); data.t = data.t(:);

    %% identify individual breaths from the raw signal using 3-pt peak detection:
    diffs_on_left_of_pt = diff(data.v); diffs_on_left_of_pt = diffs_on_left_of_pt(1:(end-1)); diffs_on_left_of_pt = logical(diffs_on_left_of_pt>0);
    diffs_on_right_of_pt = diff(data.v); diffs_on_right_of_pt = diffs_on_right_of_pt(2:end); diffs_on_right_of_pt = logical(diffs_on_right_of_pt<0);
    peaks = find(diffs_on_left_of_pt & diffs_on_right_of_pt)+1;

    %% identify individual breaths from the raw signal using 3-pt trough detection:
    diffs_on_left_of_pt = diff(-1*data.v); diffs_on_left_of_pt = diffs_on_left_of_pt(1:(end-1)); diffs_on_left_of_pt = logical(diffs_on_left_of_pt>0);
    diffs_on_right_of_pt = diff(-1*data.v); diffs_on_right_of_pt = diffs_on_right_of_pt(2:end); diffs_on_right_of_pt = logical(diffs_on_right_of_pt<0);
    troughs = find(diffs_on_left_of_pt & diffs_on_right_of_pt)+1;

    %% eliminate peaks or troughs which aren't in their half (above or below mean)
    mean_val = mean(data.v);
    temp_log = data.v(peaks) > mean_val;
    peaks = peaks(temp_log);
    temp_log = data.v(troughs) < mean_val;
    troughs = troughs(temp_log);

    %% Eliminate peaks or troughs which occur within 0.5s of the previous peak or trough
    thresh = 0.5;
    % peaks
    diffs = diff(data.t(peaks));
    rel_diffs = find(diffs > thresh); rel_diffs = rel_diffs(:);
    if isempty(rel_diffs)
        rr.v(win_no) = nan;
        continue
    end
    rel_diffs = [1; rel_diffs+1];
    peaks = peaks(rel_diffs);
    % troughs
    diffs = diff(data.t(troughs));
    rel_diffs = find(diffs > thresh); rel_diffs = rel_diffs(:);
    if isempty(rel_diffs)
        rr.v(win_no) = nan;
        continue
    end
    rel_diffs = [1; rel_diffs+1];
    troughs = troughs(rel_diffs);

    %% Eliminate peaks or troughs which aren't followed by their opposite
    comp.t = [peaks; troughs];
    comp.log = [ones(length(peaks),1); zeros(length(troughs),1)];
    [comp.t, orders] = sort(comp.t);
    comp.log = comp.log(orders);
    diffs = diff(comp.log); diffs = [1; diffs];
    comp.t = comp.t(diffs ~= 0);
    comp.log = comp.log(diffs ~= 0);
    peaks = comp.t(comp.log==1);
    troughs = comp.t(comp.log==0);

    %% Calc RR
    breaths = peaks;
    % Only using time period spanning a whole number of breaths
    win_length = data.t(breaths(end)) - data.t(breaths(1));
    rr.v(win_no) = 60*(length(breaths)-1)/win_length;

end

end

author: peterhcharlton

repo: RRest


text:
% compute abscissas and weight factors for Gaussian-Hermite quadrature
%
% CALL:  [x,w] = gauher(N)
%
%  x = base points (abscissas)
%  w = weight factors
%  N = number of base points (abscissas) (integrates an up to (2N-1)th order
%      polynomial exactly)
%
%  p(x)=exp(-x^2/2)/sqrt(2*pi), a =-Inf, b = Inf
%
%  The Gaussian Quadrature integrates a (2n-1)th order
%  polynomial exactly and the integral is of the form
%           b                         N
%          Int ( p(x)* F(x) ) dx  =  Sum ( w_j* F( x_j ) )
%           a                        j=1
%
%      this procedure uses the coefficients a(j), b(j) of the
%      recurrence relation
%
%           b p (x) = (x - a ) p   (x) - b   p   (x)
%            j j            j   j-1       j-1 j-2
%
%      for the various classical (normalized) orthogonal polynomials,
%      and the zero-th moment
%
%           1 = integral w(x) dx
%
%      of the given polynomial's weight function w(x).  Since the
%      polynomials are orthonormalized, the tridiagonal matrix is
%      guaranteed to be symmetric.

function [x,w]=gauher(N)
    if N==20 % return precalculated values
        x=[ -7.619048541679757;-6.510590157013656;-5.578738805893203;
            -4.734581334046057;-3.943967350657318;-3.18901481655339 ;
            -2.458663611172367;-1.745247320814127;-1.042945348802751;
            -0.346964157081356; 0.346964157081356; 1.042945348802751;
             1.745247320814127; 2.458663611172367; 3.18901481655339 ;
             3.943967350657316; 4.734581334046057; 5.578738805893202;
             6.510590157013653; 7.619048541679757];
        w=[  0.000000000000126; 0.000000000248206; 0.000000061274903;
             0.00000440212109 ; 0.000128826279962; 0.00183010313108 ;
             0.013997837447101; 0.061506372063977; 0.161739333984   ;
             0.260793063449555; 0.260793063449555; 0.161739333984   ;
             0.061506372063977; 0.013997837447101; 0.00183010313108 ;
             0.000128826279962; 0.00000440212109 ; 0.000000061274903;
             0.000000000248206; 0.000000000000126 ];
    else
        b = sqrt( (1:N-1)/2 )';
        [V,D] = eig( diag(b,1) + diag(b,-1) );
        w = V(1,:)'.^2;
        x = sqrt(2)*diag(D);
    end

author: benfulcher

repo: hctsa

save_path: github-repos/MATLAB/benfulcher-hctsa


text:
function OUT = NaN2Num(DATA)
% =======================================================================
% Substitute NaN values with 123456789
% =======================================================================
% OUT = Num2NaN(DATA)
% -----------------------------------------------------------------------
% INPUT
%	- DATA: a (m x n) matrix with NaNs
% -----------------------------------------------------------------------
% OUTPUT
%	- OUT: a (m x n) matrix with 123456789 instead of NaNs
% -----------------------------------------------------------------------
% EXAMPLE
%   x = [NaN 2; NaN 4; 5 6; 7 8; 9 10];
%   OUT = NaN2Num(x)
% =======================================================================
% VAR Toolbox 3.0
% Ambrogio Cesa-Bianchi
% ambrogiocesabianchi@gmail.com
% March 2012. Updated November 2020
% -----------------------------------------------------------------------

OUT=DATA;
for i=1:size(DATA,1)
    for j=1:size(DATA,2)
        if isnan(DATA(i,j)) == 1
            OUT(i,j) = 123456789;
        end
    end
end

author: ambropo

repo: VAR-Toolbox

save_path: github-repos/MATLAB/ambropo-VAR-Toolbox


text:
function []=ps_load_initial_gamma(endian)
%PS_LOAD_INITIAL_GAMMA Initial load of files into matlab workspace
%
%   Andy Hooper, Dec 2012
%
%   =======================================================================
%   01/2016 DB: Replace save with stamps_save which checks for var size when
%               saving
%   01/2016 AH: Save sensor
%   12/2016 AH: Set xy from lonlat
%   08/2017 AH: Allow for missing null ifg
%   =======================================================================

if nargin<1
   endian='b';
end

%NB IFGS assumed in ascending date order
logit;
logit('Loading data into matlab...')

phname=['./pscands.1.ph']; % for each PS candidate, a float complex value for each ifg
ijname=['./pscands.1.ij']; % ID# Azimuth# Range# 1 line per PS candidate
llname=['./pscands.1.ll']; % output from pt2geo
xyname=['./pscands.1.xy']; % output from pt2geo
hgtname=['./pscands.1.hgt']; %
daname=['./pscands.1.da']; %
rscname = ['../rsc.txt'];
pscname=['../pscphase.in'];

psver=1;
fid=fopen(rscname);
if fid<0
    error([rscname,' does not exist'])
end
rslcpar=textscan(fid,'%s');
rslcpar=rslcpar{1}{1};
fclose(fid);

fid=fopen(pscname);
if fid<0
    error([pscname,' does not exist'])
end
ifgs=textscan(fid,'%s');
fclose(fid);
ifgs=ifgs{1}(2:end);
nb=length(ifgs{1});
master_day=str2num(ifgs{1}(nb-21:nb-14));
n_ifg=length(ifgs);
n_image=n_ifg;
day=zeros(n_ifg,1);
for i=1:n_ifg
    day(i)=str2num(ifgs{i}(nb-12:nb-5));
end

year=floor(day/10000);
month=floor((day-year*10000)/100);
monthday=day-year*10000-month*100;
day=datenum(year,month,monthday);

master_day_yyyymmdd=master_day;
year=floor(master_day/10000);
month=floor((master_day-year*10000)/100);
monthday=master_day-year*10000-month*100;
master_day=datenum(year,month,monthday);

master_ix=sum(day<master_day)+1;
if day(master_ix)~=master_day
    master_master_flag='0' % no null master-master ifg provided
    day=[day(1:master_ix-1);master_day;day(master_ix:end)];
else
    master_master_flag='1' % yes, null master-master ifg provided
end

heading=readparm(rslcpar,'heading:');
setparm('heading',heading,1);

freq=readparm(rslcpar,'radar_frequency:');
lambda=299792458/freq;
setparm('lambda',lambda,1);

sensor=readparm(rslcpar,'sensor:');
if ~isempty(strfind('sensor','ASAR'))
    platform='ENVISAT';
else
    platform=sensor; % S1A for Sentinel-1A
end
setparm('platform',platform,1);

ij=load(ijname);
n_ps=size(ij,1);

rps=readparm(rslcpar,'range_pixel_spacing');
rgn=readparm(rslcpar,'near_range_slc');
se=readparm(rslcpar,'sar_to_earth_center');
re=readparm(rslcpar,'earth_radius_below_sensor');
rgc=readparm(rslcpar,'center_range_slc');
naz=readparm(rslcpar,'azimuth_lines');
prf=readparm(rslcpar,'prf');

mean_az=naz/2-0.5; % mean azimuth line

rg=rgn+ij(:,3)*rps;
look=acos((se^2+rg.^2-re^2)./(2*se*rg)); % Satellite look angles

bperp_mat=zeros(n_ps,n_image,'single');
for i=1:n_ifg
    basename=[ifgs{i}(1:nb-4),'base'];
    B_TCN=readparm(basename,'initial_baseline(TCN):',3);
    BR_TCN=readparm(basename,'initial_baseline_rate:',3);
    bc=B_TCN{2}+BR_TCN{2}*(ij(:,2)-mean_az)/prf;
    bn=B_TCN{3}+BR_TCN{3}*(ij(:,2)-mean_az)/prf;
    bperp_mat(:,i)=bc.*cos(look)-bn.*sin(look); % Convert baselines from (T)CN to perp-para coordinates
    %bpara=bc*sin(look)+bn*cos(look)
end

bperp=mean(bperp_mat)';
if master_master_flag==1
    bperp_mat=bperp_mat(:,[1:master_ix-1,master_ix+1:end]);
else
    bperp=[bperp(1:master_ix-1);0;bperp(master_ix:end)];
end
%bperp=[bperp(1:master_ix-1);0;bperp(master_ix:end)]; % insert master-master bperp (zero)
%bperp_mat=repmat(single(bperp([1:master_ix-1,master_ix+1:end]))',n_ps,1);

inci=acos((se^2-re^2-rg.^2)./(2*re*rg));
mean_incidence=mean(inci);
mean_range=rgc;

fid=fopen(phname,'r',endian);
ph=zeros(n_ps,n_ifg,'single');
byte_count=n_ps*2;
for i=1:n_ifg
    [ph_bit,byte_count]=fread(fid,[(n_ps)*2,1],'float=>single');
    ph(:,i)=complex(ph_bit(1:2:end),ph_bit(2:2:end));
end

zero_ph=sum(ph==0,2);
nonzero_ix=zero_ph<=1;       % if more than 1 phase is zero, drop node
if master_master_flag==1
    ph(:,master_ix)=1;
else
    ph=[ph(:,1:master_ix-1),ones(n_ps,1),ph(:,master_ix:end)];
    n_ifg=n_ifg+1;
    n_image=n_image+1;
end

if exist(llname,'file')
    fid=fopen(llname,'r');
    lonlat=fread(fid,[2,inf],'float',endian);
    lonlat=lonlat';
    fclose(fid);
else
    error([llname,' does not exist']);
end

ll0=(max(lonlat)+min(lonlat))/2;

xy=llh2local(lonlat',ll0)'*1000;

sort_x=sortrows(xy,1);
sort_y=sortrows(xy,2);
n_pc=round(n_ps*0.001);
bl=mean(sort_x(1:n_pc,:)); % bottom left corner
tr=mean(sort_x(end-n_pc:end,:)); % top right corner
br=mean(sort_y(1:n_pc,:)); % bottom right  corner
tl=mean(sort_y(end-n_pc:end,:)); % top left corner

theta=(180-heading)*pi/180;
if theta>pi
    theta=theta-2*pi;
end

rotm=[cos(theta),sin(theta); -sin(theta),cos(theta)];
xy=xy';
xynew=rotm*xy; % rotate so that scene axes approx align with x=0 and y=0
if max(xynew(1,:))-min(xynew(1,:))<max(xy(1,:))-min(xy(1,:)) &...
   max(xynew(2,:))-min(xynew(2,:))<max(xy(2,:))-min(xy(2,:))
    xy=xynew; % check that rotation is an improvement
    disp(['Rotating by ',num2str(theta*180/pi),' degrees']);
end

xy=single(xy');
[dummy,sort_ix]=sortrows(xy,[2,1]); % sort in ascending y order
xy=xy(sort_ix,:);
xy=[[1:n_ps]',xy];
xy(:,2:3)=round(xy(:,2:3)*1000)/1000; % round to mm

ph=ph(sort_ix,:);
ij=ij(sort_ix,:);
ij(:,1)=1:n_ps;
lonlat=lonlat(sort_ix,:);
bperp_mat=bperp_mat(sort_ix,:);

% removing Nan values if pressent in lonlat, phase
ix_nan = (sum(isnan(lonlat),2)>=1 | sum(isnan(ph),2)>=1);
lonlat(ix_nan,:)=[];
ij(ix_nan,:)=[];
xy(ix_nan,:)=[];
n_ps = size(lonlat,1);
% update the first column
ij(:,1)=[1:n_ps]';
xy(:,1)=[1:n_ps]';

savename=['ps',num2str(psver)];
stamps_save(savename,ij,lonlat,xy,bperp,day,master_day,master_ix,n_ifg,n_image,n_ps,sort_ix,ll0,master_ix,mean_incidence,mean_range);

save psver psver

phsavename=['ph',num2str(psver)];
ph(ix_nan,:)=[];
stamps_save(phsavename,ph);

bpsavename=['bp',num2str(psver)];
bperp_mat(ix_nan,:)=[];
stamps_save(bpsavename,bperp_mat);

lasavename=['la',num2str(psver)];
la=inci(sort_ix); % store incidence not look angle for gamma
la(ix_nan,:)=[];
stamps_save(lasavename,la);

if exist(daname,'file')
  D_A=load(daname);
  D_A=D_A(sort_ix);
  D_A(ix_nan,:)=[];
  dasavename=['da',num2str(psver)];
  stamps_save(dasavename,D_A);
end

if exist(hgtname,'file')
    fid=fopen(hgtname,'r');
    hgt=fread(fid,[1,inf],'float',endian);
    hgt=hgt';
    hgt=hgt(sort_ix);
    fclose(fid);
    hgt(ix_nan,:)=[];
    hgtsavename=['hgt',num2str(psver)];
    stamps_save(hgtsavename,hgt);
end

logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function test_issue840

% WALLTIME 00:10:00
% MEM 2gb
% DEPENDENCY ft_megplanar ft_combineplanar ft_appenddata

%%

if false
  % let's make three small data structures
  % the data is from https://openneuro.org/datasets/ds000117/versions/1.0.3

  filename = {
    'sub-01_ses-meg_task-facerecognition_run-01_meg.fif'
    'sub-01_ses-meg_task-facerecognition_run-02_meg.fif'
    'sub-01_ses-meg_task-facerecognition_run-03_meg.fif'
    };

  for i=1:numel(filename)
    cfg = [];
    cfg.trialfun = 'ft_trialfun_general';
    cfg.trialdef.triallength = 1;
    cfg.trialdef.ntrials = 10;
    cfg.dataset = filename{i};
    cfg = ft_definetrial(cfg);
    data{i} = ft_preprocessing(cfg);
  end

  data1 = data{1};
  data2 = data{2};
  data3 = data{3};
  clear data
  save data1 data1
  save data2 data2
  save data3 data3

else
  % load it from the mat files
  load(dccnpath('/home/common/matlab/fieldtrip/data/test/issue840/data1.mat'));
  load(dccnpath('/home/common/matlab/fieldtrip/data/test/issue840/data2.mat'));
  load(dccnpath('/home/common/matlab/fieldtrip/data/test/issue840/data3.mat'));

end

%%

cfg = [];
data1_combined = ft_combineplanar(cfg, data1);
data2_combined = ft_combineplanar(cfg, data2);
data3_combined = ft_combineplanar(cfg, data3);

cfg = [];
data_combined_append = ft_appenddata(cfg, data1_combined, data2_combined, data3_combined);
assert( isfield(data_combined_append, 'elec')); % this one is consistent
assert(~isfield(data_combined_append, 'grad')); % this one is not consistent

%%

cfg = [];
data_append = ft_appenddata(cfg, data1, data2, data3);
assert( isfield(data_append, 'elec')); % this one is consistent
assert(~isfield(data_append, 'grad')); % this one is not consistent

%%

cfg = [];
cfg.channel = 'MEG*3';
cfg.method = 'distance';
cfg.grad = ft_convert_units(data1.grad, 'mm');
cfg.neighbourdist = 40;
neighbours = ft_prepare_neighbours(cfg);

%%

cfg = [];
cfg.channel = 'MEG*3';
cfg.neighbours = neighbours;

% the FT_MEGPLANAR function refuses to work on neuromag306 data. Here we carefully
% selected only the magnetometers (both as cfg.channel and in cfg.neighbours) so in
% principle it should work fine. The hack I am using is to make FT_MEGPLANAR think it
% is an itab153 sensor array
data1.grad.type = 'itab153';
data2.grad.type = 'itab153';
data3.grad.type = 'itab153';

data1_planar = ft_megplanar(cfg, data1);
data2_planar = ft_megplanar(cfg, data2);
data3_planar = ft_megplanar(cfg, data3);

%%

cfg = [];
data_planar_append = ft_appenddata(cfg, data1_planar, data2_planar, data3_planar);
assert( isfield(data_planar_append, 'elec')); % this one is consistent
assert(~isfield(data_planar_append, 'grad')); % this one is not consistent

%%

cfg = [];
cfg.channel = {'*_dV', '*_dV'};
data1_planar_combined = ft_combineplanar(cfg, data1_planar);
data2_planar_combined = ft_combineplanar(cfg, data2_planar);
data3_planar_combined = ft_combineplanar(cfg, data3_planar);

%%

cfg = [];
data_planar_combined_append = ft_appenddata(cfg, data1_planar_combined, data2_planar_combined, data3_planar_combined);
assert( isfield(data_planar_combined_append, 'elec')); % this one is consistent
assert(~isfield(data_planar_combined_append, 'grad')); % this one is not consistent

%%

% this fails because it cannot determine the sensor type: the data contains elec, but no grad

% However, if it were to contain grad, it would probably also fail since I modified
% the grad.type above from neuromag306 into itab153. Overall, this is a pretty crappy
% test script. I'll keep it nevertheless, since there are some usefiel parts in it.

if false
  cfg = [];
  cfg.channel = {'*_dV', '*_dV'};
  data_planar_append_combined = ft_combineplanar(cfg, data_planar_append);
end

% Error using ft_checkdata (line 573)
% This function requires data with an 'ctf151_planar', 'ctf275_planar', 'neuromag122',
% 'neuromag306', 'bti248_planar', 'bti148_planar', 'itab153_planar', 'yokogawa160_planar',
% 'yokogawa64_planar' or 'yokogawa440_planar' sensor array.

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function [x,k,res,negCurv] = cg(A,b,optTol,maxIter,verbose,precFunc,precArgs,matrixVectFunc,matrixVectArgs)
% [x,k,res,negCurv] =
% cg(A,b,optTol,maxIter,verbose,precFunc,precArgs,matrixVectFunc,matrixVect
% Args)
% Linear Conjugate Gradient, where optionally we use
% - preconditioner on vector v with precFunc(v,precArgs{:})
% - matrix multipled by vector with matrixVectFunc(v,matrixVectArgs{:})

x = zeros(size(b));
r = -b;

% Apply preconditioner (if supplied)
if nargin >= 7 && ~isempty(precFunc)
    y = precFunc(r,precArgs{:});
else
    y = r;
end

ry = r'*y;
p = -y;
k = 0;

res = norm(r);
done = 0;
negCurv = [];
while res > optTol & k < maxIter & ~done
    % Compute Matrix-vector product
    if nargin >= 9
        Ap = matrixVectFunc(p,matrixVectArgs{:});
    else
        Ap = A*p;
    end
    pAp = p'*Ap;

    % Check for negative Curvature
    if pAp <= 1e-16
        if verbose
            fprintf('Negative Curvature Detected!\n');
        end

        if nargout == 4
           if pAp < 0
              negCurv = p;
              return
           end
        end

        if k == 0
            if verbose
                fprintf('First-Iter, Proceeding...\n');
            end
            done = 1;
        else
            if verbose
                fprintf('Stopping\n');
            end
            break;
        end
    end

    % Conjugate Gradient
    alpha = ry/(pAp);
    x = x + alpha*p;
    r = r + alpha*Ap;

    % If supplied, apply preconditioner
    if nargin >= 7 && ~isempty(precFunc)
        y = precFunc(r,precArgs{:});
    else
        y = r;
    end

    ry_new = r'*y;
    beta = ry_new/ry;
    p = -y + beta*p;
    k = k + 1;

    % Update variables
    ry = ry_new;
    res = norm(r);
end
end

author: huashiyiqike

repo: LSTM-MATLAB


text:
function [cMap]=bloodbone(varargin)

switch nargin
    case 0
        n=250;
    case 1
        n=varargin{1};
end

b=blood(8*20);
bb=linspacen([1 0.95 0.8],b(1,:),3*20)';
cMap=[bb(1:end-1,:); b];
[cMap]=resampleColormap(cMap,n);

end

%%
% _*GIBBON footer text*_
%
% License: <https://github.com/gibbonCode/GIBBON/blob/master/LICENSE>
%
% GIBBON: The Geometry and Image-based Bioengineering add-On. A toolbox for
% image segmentation, image-based modeling, meshing, and finite element
% analysis.
%
% Copyright (C) 2006-2022 Kevin Mattheus Moerman and the GIBBON contributors
%
% This program is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program.  If not, see <http://www.gnu.org/licenses/>.

author: gibbonCode

repo: GIBBON

save_path: github-repos/MATLAB/gibbonCode-GIBBON


text:
function y=str2strs(s)

% this function takes s (a string deliminated by tabs, commas, etc.)
%  and breaks it up into sub strings stored in y (a cell array of strings)
%
% spaces are eliminated by the function strcat so these are replaced in the
%  beginning by an underscore. space delimination can be handled, though s
%  should be checked for any pre-existing underscrores.
%
% recovery of the individual elements is by using the char command
%
% e.g., applying it to a tab deliminated string...
%   fields = str2strs('abc	d%*	34_EF	5	s p a c e s')     returns
%   'abc'    'd%*'    '34_EF'    '5'    's_p_a_c_e_s'
% and then
%   char(fields(3))     returns
%   34_EF
%
% David Malicky  malicky@umich.edu

% replace each space w/ a character that strcat won't miss
s=strrep(s,' ','_');
delim='	'; % that long space is a tab, though any character may be used

lens=length(s);
tabs=findstr(s,delim);
numtabs=length(tabs);	% this is 1 less than the number of sub-strings

i=0;
for j=1:numtabs+1
	y(j)=cellstr('');
	while ( (i<lens) & (s(i+1)~=delim)  )  % another tab
		y(j)=strcat(y(j),s(i+1));
		i=i+1;
	end;
	i=i+1;
end;

y;

author: Sable

repo: mcbench-benchmarks


text:
function [newnodes, len]=polylinesimplify(nodes, minangle)
%
% [newnodes, len]=polylinesimplify(nodes, minangle)
%
% Calculate a simplified polyline by removing nodes where two adjacent
% segment have an angle less than a specified limit
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    node: an N x 3 array defining each vertex of the polyline in
%          sequential order
%    minangle:(optional) minimum segment angle in radian, if not given, use
%          0.75*pi
%
% output:
%    newnodes: the updated node list; start/end will not be removed
%    len: the length of each segment between the start and the end points
%
%
% -- this function is part of brain2mesh toolbox (http://mcx.space/brain2mesh)
%    License: GPL v3 or later, see LICENSE.txt for details
%

if(nargin<2)
    minangle=0.75*pi;
end

v=segvec(nodes(1:end-1,:), nodes(2:end,:));
ang=acos(max(min(sum(-v(1:end-1,:).*(v(2:end,:)),2),1),-1));

newnodes=nodes;
newv=v;
newang=ang;

idx=find(newang<minangle);

while(~isempty(idx))
    newnodes(idx+1,:)=[];
    newv(idx+1,:)=[];
    newang(idx)=[];
    idx=unique(idx-(0:(length(idx)-1))');
    idx1=idx(idx<size(newnodes,1));
    newv(idx1,:)  =segvec(newnodes(idx1,:),newnodes(idx1+1,:));
    idx1=idx(idx<size(newv,1));
    newang(idx1)  =acos(sum(-newv(idx1,:).*(newv(idx1+1,:)),2));
    idx0=idx(idx>1);
    newang(idx0-1)=acos(sum(-newv(idx0-1,:).*(newv(idx0,:)),2));
    idx=find(newang<minangle);
end

if(nargout>1)
    len=newnodes(1:end-1,:) - newnodes(2:end,:);
    len=sqrt(sum(len.*len,2));
end

function v=segvec(n1, n2)

v=n2-n1;
normals=sqrt(sum(v.*v,2));
v=v./repmat(normals,1,size(v,2));

author: fangq


text:
function intlablogo(angle,screensaver)
%INTLABLOGO   Show INTLAB logo
%
%  intlablogo
%
%put cursor into figure and press
%
%  o  open
%  c  close
%  q  quit
%
%The call
%
%  intlablogo(angles)
%
%shows the dodecahedron in opening angles.
%

%The call
%   intlablogo([],1)
%generates a screensaver picture using
%  generate file: print -djpeg100 -zbuffer picturefile
%

% written  10/16/98     S.M. Rump
% modified 04/04/04     S.M. Rump  revision
% modified 10/13/05     S.M. Rump  name and fonts adapted
% modified 09/06/07     S.M. Rump  version number
% modified 06/24/08     S.M. Rump  screensaver
%

  if nargin<=1
    screensaver = 0;
  end

  psi1 = 10.8123169636;
  psi2 = 52.6226318594;

  for i=1:5
    P{i} = polar2rect(1,72*i,psi2);
    P{5+i} = polar2rect(1,72*i,psi1);
    P{10+i} = polar2rect(1,72*i-36,-psi1);
    P{15+i} = polar2rect(1,72*i-36,-psi2);
  end

  M = zeros(18,5,3);

  M(1,:,:) = [ P{ 1} P{ 5} P{10} .5*(P{10}+P{11}) .5*(P{1}+P{6}) ]';
  M(2,:,:) = [ P{ 1} P{ 2} P{ 3} P{ 4} P{ 5} ]';
  M(3,:,:) = [ P{ 5} P{ 4} P{ 9} P{15} P{10} ]';
  M(4,:,:) = [ P{ 3} P{ 8} P{14} P{ 9} P{ 4} ]';
  M(5,:,:) = [ P{10} P{15} .5*(P{15}+P{20}) .5*(P{10}+P{11}) P{10}]';
  M(6,:,:) = [ P{ 1} P{ 2} .5*(P{2}+P{7}) .5*(P{6}+P{1}) P{1}]';
  M(7,:,:) = [ P{15} P{ 9} P{14} .5*(P{14}+P{19}) .5*(P{20}+P{15}) ]';
  M(8,:,:) = [ P{2} P{3} P{8} .5*(P{8}+P{13}) .5*(P{2}+P{7}) ]';
  M(9,:,:) = [ P{8} P{14} .5*(P{14}+P{19}) .5*(P{8}+P{13}) P{8} ]';
  M(10:18,:,:) = -M(1:9,:,:);

  dil = .25*( P{14}+P{19} + P{20}+P{15} );
  for i=1:3
    M(:,:,i) = M(:,:,i) - dil(i);
  end
  phi = atan2(dil(2),dil(1));
  T = [ cos(phi) -sin(phi) 0 ; ...
      sin(phi)  cos(phi) 0 ; ...
      0          0     1 ];
  M = reshape( reshape(M,90,3)*T , 18,5,3 );
  Az = 23.75; El = 12.5;
  close

  % show a couple of angles
  if nargin~=0
    if ~isempty(angle)
      psi = angle(1);
      omegaold = 0;
      for omega=angle
        psi = omega-omegaold;
        M = showlogo(psi,M,screensaver);
        omegaold = omega;
        pause(0.05)
      end
      return
    end
  end

  % open/close by hand
  psi = 30;
  delta = 5;
  omega = psi;
  change = 1;

  while change
    clf
    M = showlogo(psi,M,screensaver);
    if screensaver
      print -djpeg100 -zbuffer picturefile
      close
      return
    end
    change = 0;

    if ~screensaver
      title( [ 'move cursor into window and press  "o"  or  "c"  or,  "q" for quit' ], ...
               'FontName','Sans Serif','Fontsize',10 );
    end
    shg
    waitforbuttonpress
    switch get(gcf,'CurrentCharacter')
      case 'o',
        if omega <= 180,
          psi = delta; omega = omega + delta;
        else
          psi = 0;
        end
        change = 1;
      case 'c',
        if omega >= delta
          psi = -delta; omega = omega - delta;
        else
          psi = 0;
        end
        change = 1;
    end
  end
  close

function M = showlogo(psi,M,screensaver)
% opens logo by angle psi from current setting
  clf
  axis off
  axis equal
  if screensaver
    M = 1.25*M;
    axis([-2.4 1.2 -1.5 2.1 -1 1.8])
  else
    axis([-1.5 1.5 -1.7 1.7 -.8 1.5])
  end
  axis manual
  hold on
  Az = 23.75; El = 12.5;
  view(Az,El)

  for i=10:18
    rand('state',0);
    fill3(M(i,:,1),M(i,:,2),M(i,:,3),rand(1,5,3));
  end

  psi_ = psi/180*pi;
  T = [ cos(psi_) 0 -sin(psi_) ; ...
      0     1      0      ; ...
      sin(psi_) 0  cos(psi_) ];
  M(1:9,:,:) = reshape( reshape(M(1:9,:,:),45,3)*T , 9,5,3 );

  for i=1:9
    rand('state',0);
    fill3(M(i,:,1),M(i,:,2),M(i,:,3),rand(1,5,3));
  end
  addtext(screensaver)
  hold off

function x = polar2rect(r,phi,psi)     % in degrees
  x = zeros(3,1);
  sigma = pi/180;
  phi = sigma*phi;
  psi = sigma*psi;
  x(1) = cos(phi).*cos(psi);
  x(2) = sin(phi).*cos(psi);
  x(3) = sin(psi);
  x = r*x;

function addtext(screensaver)
  if screensaver
    text(-1,-4  ,['INTLAB  -  INTerval LABoratory (Version ' intvalinit('version',0) ')'], ...
      'Fontsize',12, ...
      'FontName','MS Reference Sans Serif', ...
      'HorizontalAlignment','left');
    text(-.42,-5.2,'The Matlab toolbox for Reliable Computing  -  www.ti3.tuhh.de/rump', ...
      'Fontsize',10, ...
      'FontName','Arial', ...
      'HorizontalAlignment','left');
    text(-.36,-6.2,'Siegfried M. Rump, Institute for Reliable Computing, Hamburg University of Technology  ', ...
      'Fontsize',10, ...
      'FontName','Arial', ...
      'HorizontalAlignment','left');
  else
    text(-1,-4  ,['INTLAB  -  INTerval LABoratory (Version ' intvalinit('version',0) ')'], ...
      'Fontsize',17, ...
      'FontName','MS Reference Sans Serif', ...
      'HorizontalAlignment','left');
    text(-.1,-5.2,'The Matlab toolbox for Reliable Computing  -  www.ti3.tuhh.de/rump', ...
      'Fontsize',10, ...
      'FontName','Arial', ...
      'HorizontalAlignment','left');
    text(-.02,-6.2,'Siegfried M. Rump, Institute for Reliable Computing, Hamburg University of Technology  ', ...
      'Fontsize',10, ...
      'FontName','Arial', ...
      'HorizontalAlignment','left');
  end
  if ~screensaver
    str(1) =    {' x = hessianinit(xs);'};
    str(2) =     {' y = f(x);'};
    str(3) = {' xs = xs - y.hx\\y.dx'';'};
    text(-.5, 7.2,str, ...
      'HorizontalAlignment','right', ...
      'Fontsize',13, ...
      'FontName','Courier');
  end

repo: OpenMAS


text:
function [im]=mrgb2gray(im,met)
% MRGB2GRAY - convert RGB images to grayscale
%
% im = mrgb2gray(M,met);
%
% M - image matrix or filename. The size of the image should be
% MxNx3 when it is loaded by imread.
%
% MET - method for conversion. Should be one of
%
%       'default'   - uses a weight of [0.3 0.59 0.11] for the R, G and
%                     B respectively
%       'max'       - uses the maximum value of the R, G or B for
%                     each pixel
%       'mean'      - uses the mean of the RGB for each pixel. This
%                     corresponds to using a weight of
%                     [0.33 0.33 0.33]
%       'median'    - as above except it uses the median value
%       'min'       - as above except it uses the minimum value
%       'desat'     - computes (max(R,G,B)+min(R,G,B))/2 for each
%                     pixel. This produces a lower contrast image.
%
% MET can also be a 1x3 array with weights.
%
% examples:
%           im=mrgb2gray(im,'max');
%           im=mrgb2gray(im);
%           im=mrgb2gray(im,[0.3 0.5 0.2]);
%           im=mrgb2gray('filename.jpg');
%           im=mrgb2gray('filename.jpg',[0.25 0.60 0.15]);
%           im=mrgb2gray('filename.jpg','desat');
%
% See also: rgb2gray, isrgb, misrgb
%

% version 0.22 - bug fix: under 'mean' changed occurence of 'min' to 'mean'
% version 0.21 - cleanup of code
% 24.08.2006
%
% version 0.2 - first released version
% J.K.Sveen@damtp.cam.ac.uk, 2004, August 27.
% Distributed under the terms of the GNU GPL:
% http://www.gnu.org/copyleft/gpl.html

if nargin==1
  met='default';
end

if ischar(im)
  im=imread(im);
end

[sx,sy,sz]=size(im);

if sz~=3
  disp('Error. This image is not RGB. Aborting.'); return
end

wasuint=isa(im,'uint8');
im=double(im);
if ~ischar(met)
  T=met(:);
  im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
  im=im*T; % multiply the weights
  im=reshape(im,sx,sy); % reshape the image back to correct size
else
  switch lower(met)
   case {'default','gimp','standard'}

    T=[0.3 0.59 0.11]'; % weights for the individual colors. These are
			% apparently the same weights as used in GIMP:
			%
			% http://gimp-savvy.com/BOOK/index.html?node54.html

    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=im*T; % multiply the weights, 0.3*R, 0.59*G, 0.11*B
    im=reshape(im,sx,sy); % reshape the image back to correct size
   case {'max','Max'}
    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=max(im,[],2); %rotate image and take mean of all columns, then
                  %rotate back. We now have an M*N by 1 array
    im=reshape(im,sx,sy); % reshape the image back to size M by N
   case {'mean','average'}
    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=mean(im,2);
    im=reshape(im,sx,sy); % reshape the image back to correct size
   case {'median'}
    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=median(im,2);
    im=reshape(im,sx,sy); % reshape the image back to correct size
   case {'min'}
    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=min(im,[],2);
    im=reshape(im,sx,sy); % reshape the image back to correct size
   case {'desat','desaturate'}
    im=reshape(im,sx*sy,sz); %put image in an M*N by 3 array
    im=(max(im,[],2) + min(im,[],2))/2;
    im=reshape(im,sx,sy); % reshape the image back to correct size
   otherwise
    disp('Unknown method. Check your input.'); return
  end
end

if wasuint % change back to uint8 if that was the input
  im=uint8(im);
end


text:
function y = snakeindex(IDX)
% SNAKEINDEX  Create index for adpative interpolating the snake
%     y = snakeindex(IDX)
%

N = length(IDX);
y=1:0.5:N+0.5;
x=1:N;
y(2*x(IDX==0))=[];


text:
% Jiao Xianjun (putaoshu@msn.com; putaoshu@gmail.com)

% Convert input signal(with sampling period shorter than expectation by
% ppm_val PPM) to output signal(sampling period error corrected)

% A script of project: https://github.com/JiaoXianjun/rtl-sdr-LTE

function r = sampling_period_correction(s, ppm_val)
sampling_time_in = ( 0 : (length(s)-1) ).*(1-ppm_val*1e-6);
sampling_time_out = 0 : (length(s)-1);
if ppm_val > 0
    ep =  floor( (length(s)-1).*(1-ppm_val*1e-6) );
    sampling_time_out = sampling_time_out(1: ep);
end
r = interp1(sampling_time_in, s, sampling_time_out, 'nearest').';

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function R = gqr3 (A)
%GQR3 QR factorization, based on Givens rotations
%
% Example:
%   R = gqr3 (A)
% See also: testall

%   Copyright 2006-2007, Timothy A. Davis.
%   http://www.cise.ufl.edu/research/sparse

[m n] = size (A) ;

% parent = cs_etree (sparse (A), 'col') ;

for i = 2:m
    % i
    for k = 1:min(i-1,n)
    % k
        % Givens rotation to zero out A(i,k) using A(k,k)
        G = givens2 (A(k,k), A(i,k)) ;
        A ([k i],k:n) = G * A ([k i],k:n) ;
        A (i,k) = 0 ;
        % fprintf ('A(21,25)=%g\n', A(21,25)) ;
        % if (A(21,25) ~= 0)
            % pause
        % end
    end
end
R = A ;


text:
function vline(x, varargin)

% VLINE plot a vertical line in the current graph

abc = axis;
x = [x x];
y = abc([3 4]);
if length(varargin)==1
  varargin = {'color', varargin{1}};
end
h = line(x, y);
set(h, varargin{:});


text:
% HNORMALISE - Normalises array of homogeneous coordinates to a scale of 1
%
% Usage:  nx = hnormalise(x)
%
% Argument:
%         x  - an Nxnpts array of homogeneous coordinates.
%
% Returns:
%         nx - an Nxnpts array of homogeneous coordinates rescaled so
%              that the scale values nx(N,:) are all 1.
%
% Note that any homogeneous coordinates at infinity (having a scale value of
% 0) are left unchanged.

% Peter Kovesi
% School of Computer Science & Software Engineering
% The University of Western Australia
% http://www.csse.uwa.edu.au/~pk
%
% February 2004

function nx = hnormalise(x)

    [rows,npts] = size(x);
    nx = x;

    % Find the indices of the points that are not at infinity
    finiteind = find(abs(x(rows,:)) > eps);

    %if length(finiteind) ~= npts
    %    warning('Some points are at infinity');
    %end

    % Normalise points not at infinity
    for r = 1:rows-1
	nx(r,finiteind) = x(r,finiteind)./x(rows,finiteind);
    end
    nx(rows,finiteind) = 1;

author: jianxiongxiao

repo: ProfXkit


text:
  function ob = strum(data, methods, varargin)
%|function ob = strum(data, methods, [options])
%|
%| Construct a "strum" object.  This acts like an ordinary matlab "struct"
%| except that it also can contain user-defined "methods."
%| For example, if x is a fieldname of the structure, then ob.x will return
%| the corresponding value, just like an ordinary structure.
%| But x might instead be a handle to a function, say @xfun,
%| and invoking ob.x will return xfun(st).
%| Furthermore, those functions can accept arguments, i.e., ob.x(arg[s]),
%| which will invoke xfun(st, arg[s]).
%| The invoked functions are passed the structure (really, the entire object)
%| so that the invoked functions can access all the "not so private" data
%| in the structure.
%|
%| in
%|	data	struct	initial structure, i.e., 'data'
%|	methods	cell	{'name1', @handle1, 'name2', @handle2, ...}
%|		or	{'name1', @handle1, 'doc1; ...} (documented version)
%|
%| option
%|	'base'	strum	add data / methods to an existing strum
%|
%| out
%|	ob		strum object
%|
%| Copyright 2006-1-19, Jeff Fessler, University of Michigan

if nargin == 1 && streq(data, 'test'), run_mfile_local strum_test; return, end

ob.caller = caller_name(1);
ob.meth = struct;
ob.data = struct;
ob.docs = {};

if nargin == 0	% required by Mathworks
	ob = class(ob, mfilename);
return
end

if nargin < 2, ir_usage, end

if nargin > 2 % base object provided
	arg.base = struct([]);
	arg = vararg_pair(arg, varargin);
	if isempty(arg.base), keyboard, fail('base required'), end

	ob = struct(arg.base); % initialize with base data/methods
	ob.caller = [ob.caller '/' caller_name(1)]; % append caller

	% add new data to base object
	if ~isstruct(data), fail('data must be a struct'), end
	names = fieldnames(data);
	for ii=1:length(names)
		name = names{ii};
		if isfield(ob.data, name) || isfield(ob.meth, name)
			warn(['name reuse: ', name])
		end
		ob.data.(name) = data.(name); % set or over-ride with new
	end

	method = ob.meth;
	has_base = true;

else
	ob.data = data;
	method = struct;%([]); % no methods
	has_base = false;
end

%if isempty(methods)
%	warning 'empty methods: why use strum instead of struct?'
%end

if 3 == size(methods,2) % name,handle,doc triples

	if has_base && isempty(ob.docs)
		error 'docs must be consistent between base and new'
	end

	for ii=1:size(methods,1)
		name = methods{ii,1};
		handle = methods{ii,2};
		doc = methods{ii,3};
		if isfield(ob.data, name) || isfield(method, name)
			warn('name reuse: %s', name)
		end
		method.(name) = handle;
		ob.docs{end+1} = doc;

		strum_check_nargout(handle)
	end

elseif 0 == rem(length(methods),2) % name,handle pairs

	if has_base && ~isempty(ob.docs)
		error 'docs must be consistent between base and new'
	end

	for ii=1:2:length(methods)
		name = methods{ii};
		handle = methods{ii+1};
		if isfield(ob.data, name) || isfield(method, name)
			warn('name reuse: %s', name)
		end
		method.(name) = handle;
		strum_check_nargout(handle)
	end

else
	error 'need name,value pairs, or name,value,doc; triples'
end

ob.meth = method;

ob = class(ob, mfilename);

% strum_check_nargout()
% todo: warn if not exactly one output?
function strum_check_nargout(handle)
nout = ir_nargout(handle);
if nout == -1
	return % in matlab, -1 means "variable number of arguments"
end

if nout ~= 1
	pr handle
	keyboard
end

author: JeffFessler

repo: mirt

save_path: github-repos/MATLAB/JeffFessler-mirt


text:
function t = isdigit(c)
%ISDIGIT True for decimal digits.
%
%   For a string C, ISDIGIT(C) is 1 for decimal digits and 0 otherwise.
%
%   (In ASCII: characters [0-9].)
%
%   See also: ISALNUM, ISALPHA, ISASCII, ISLOWER, ISPRTCHR, ISUPPER,
%   ISXDIGIT.

%   Author:      Peter John Acklam
%   Time-stamp:  2002-03-03 13:51:04 +0100
%   E-mail:      pjacklam@online.no
%   URL:         http://home.online.no/~pjacklam

   if verLessThan('matlab','9.0')
       error(nargchk(1, 1, nargin));
   else
       narginchk(1, 1);
   end

   t = ischar(c) & ( '0' <= c ) & ( c <= '9' );

author: BUNPC


text:
function [ y ] = vl_nnselection( X,idx,dzdy )
% SelectionLayer
% X	: 1 x 2 x n x b
% y : 1 x 2 x nidx x b

if nargin<3
    %forward
    y = zeros(1,size(X,2),length(idx),size(X,4),'single');
    y(1,:,:,:) = X(1,:,idx,:);
else
    %backward
    y = zeros(size(X),'single');
    y(1,:,idx,:) = dzdy;
end

end

author: anilbas


text:
function varargout = getPrmDflt( prm, dfs, checkExtra )
% Helper to set default values (if not already set) of parameter struct.
%
% Takes input parameters and a list of 'name'/default pairs, and for each
% 'name' for which prm has no value (prm.(name) is not a field or 'name'
% does not appear in prm list), getPrmDflt assigns the given default
% value. If default value for variable 'name' is 'REQ', and value for
% 'name' is not given, an error is thrown. See below for usage details.
%
% USAGE (nargout==1)
%  prm = getPrmDflt( prm, dfs, [checkExtra] )
%
% USAGE (nargout>1)
%  [ param1 ... paramN ] = getPrmDflt( prm, dfs, [checkExtra] )
%
% INPUTS
%  prm          - param struct or cell of form {'name1' v1 'name2' v2 ...}
%  dfs          - cell of form {'name1' def1 'name2' def2 ...}
%  checkExtra   - [0] if 1 throw error if prm contains params not in dfs
%                 if -1 if prm contains params not in dfs adds them
%
% OUTPUTS (nargout==1)
%  prm    - parameter struct with fields 'name1' through 'nameN' assigned
%
% OUTPUTS (nargout>1)
%  param1 - value assigned to parameter with 'name1'
%   ...
%  paramN - value assigned to parameter with 'nameN'
%
% EXAMPLE
%  dfs = { 'x','REQ', 'y',0, 'z',[], 'eps',1e-3 };
%  prm = getPrmDflt( struct('x',1,'y',1), dfs )
%  [ x y z eps ] = getPrmDflt( {'x',2,'y',1}, dfs )
%
% See also INPUTPARSER
%
% Piotr's Computer Vision Matlab Toolbox      Version 2.60
% Copyright 2014 Piotr Dollar.  [pdollar-at-gmail.com]
% Licensed under the Simplified BSD License [see external/bsd.txt]

if( mod(length(dfs),2) ), error('odd number of default parameters'); end
if nargin<=2, checkExtra = 0; end

% get the input parameters as two cell arrays: prmVal and prmField
if iscell(prm) && length(prm)==1, prm=prm{1}; end
if iscell(prm)
  if(mod(length(prm),2)), error('odd number of parameters in prm'); end
  prmField = prm(1:2:end); prmVal = prm(2:2:end);
else
  if(~isstruct(prm)), error('prm must be a struct or a cell'); end
  prmVal = struct2cell(prm); prmField = fieldnames(prm);
end

% get and update default values using quick for loop
dfsField = dfs(1:2:end); dfsVal = dfs(2:2:end);
if checkExtra>0
  for i=1:length(prmField)
    j = find(strcmp(prmField{i},dfsField));
    if isempty(j), error('parameter %s is not valid', prmField{i}); end
    dfsVal(j) = prmVal(i);
  end
elseif checkExtra<0
  for i=1:length(prmField)
    j = find(strcmp(prmField{i},dfsField));
    if isempty(j), j=length(dfsVal)+1; dfsField{j}=prmField{i}; end
    dfsVal(j) = prmVal(i);
  end
else
  for i=1:length(prmField)
    dfsVal(strcmp(prmField{i},dfsField)) = prmVal(i);
  end
end

% check for missing values
if any(strcmp('REQ',dfsVal))
  cmpArray = find(strcmp('REQ',dfsVal));
  error(['Required field ''' dfsField{cmpArray(1)} ''' not specified.'] );
end

% set output
if nargout==1
  varargout{1} = cell2struct( dfsVal, dfsField, 2 );
else
  varargout = dfsVal;
end

author: pdollar

repo: toolbox

save_path: github-repos/MATLAB/pdollar-toolbox


text:
function [aa,z] = fitnesstra3f(pop,para)
[pops,numvar]=size(pop);
for i=1:pops
    [phi,conff]=invkin3(para(4),para(5),pop(i,4)) ;%final conficuration
    z(i)=phi;
    cond=[para(1:3),conff];
    chrom=[pop(i,1:3),pop(i,5:end)];
    kk=trajt3(cond,chrom);
    tt=torque3(kk);
    ft=ftorque3(tt);
    poo=kk(1:3,:);
    fq=sum(sum(abs(diff(poo'))));
    pos=forkin3(kk(1,:),kk(2,:),kk(3,:));xx=pos(1,:);yy=pos(2,:);
    x=diff(xx);
    y=diff(yy);
    dis=sqrt(x.^2+y.^2);
    fdis=sum(dis);
    time=pop(i,8)+pop(i,9);
    a1=1.7;
    a2=2;
    a3=2;
    a4=1;
    fit=a1*ft+a2*fq+a3*fdis+a4*time;
    aa(i)=1/fit;
end


text:
function net = initializeFaceCNN_simple_5(num_bins)

f=1/100 ;
net.layers = {} ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(9,9,1,10, 'single'), zeros(1, 10, 'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(7,7,10,10, 'single'), zeros(1,10,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'pool', ...
                           'method', 'max', ...
                           'pool', [2 2], ...
                           'stride', 2, ...
                           'pad', 0) ;
% This is basically an FC layer
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(10,10,10,50, 'single'), zeros(1,50,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'relu') ;
net.layers{end+1} = struct('type', 'conv', ...
                           'weights', {{f*randn(1,1,50,num_bins, 'single'), zeros(1,num_bins,'single')}}, ...
                           'stride', 1, ...
                           'pad', 0) ;
net.layers{end+1} = struct('type', 'softmaxloss') ;

net = vl_simplenn_tidy(net) ;

author: TadasBaltrusaitis

repo: OpenFace


text:
function r = ldivide(a,b)
% o .\ v
%
% Syntax
%   h = o .\ r
%
% Input
%  o - @orientation
%  r - @vector3d
%
% Output
%  h - @Miller indice
%
% See also
%

r = inv(a) .* b;

author: mtex-toolbox

repo: mtex

save_path: github-repos/MATLAB/mtex-toolbox-mtex


text:
function varargout = process_fooof(varargin)
% PROCESS_FOOOF: Applies the "Fitting Oscillations and One Over F" (specparam) algorithm on a Welch's PSD
%
% REFERENCE: Please cite the original algorithm:
%    Donoghue T, Haller M, Peterson E, Varma P, Sebastian P, Gao R, Noto T,
%    Lara AH, Wallis JD, Knight RT, Shestyuk A, Voytek B. Parameterizing
%    neural power spectra into periodic and aperiodic components.
%    Nature Neuroscience (2020)

% @=============================================================================
% This function is part of the Brainstorm software:
% https://neuroimage.usc.edu/brainstorm
%
% Copyright (c) University of Southern California & McGill University
% This software is distributed under the terms of the GNU General Public License
% as published by the Free Software Foundation. Further details on the GPLv3
% license can be found at http://www.gnu.org/copyleft/gpl.html.
%
% FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
% UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
% WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
% MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
% LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
%
% For more information type "brainstorm license" at command prompt.
% =============================================================================@
%
% Authors: Luc Wilson, Francois Tadel, 2020-2022

eval(macro_method);
end

%% ===== GET DESCRIPTION =====
function sProcess = GetDescription() %#ok<DEFNU>
    % Description the process
    sProcess.Comment     = 'specparam: Fitting oscillations and 1/f';
    sProcess.Category    = 'Custom';
    sProcess.SubGroup    = 'Frequency';
    sProcess.Index       = 490;
    sProcess.Description = 'https://neuroimage.usc.edu/brainstorm/Tutorials/Fooof';
    % Definition of the input accepted by this process
    sProcess.InputTypes  = {'timefreq'};
    sProcess.OutputTypes = {'timefreq'};
    sProcess.nInputs     = 1;
    sProcess.nMinFiles   = 1;
    % Definition of the options
    % === FOOOF TYPE
    sProcess.options.implementation.Comment = {'Matlab', 'Python 3 (3.7 recommended)', 'specparam implementation:'; 'matlab', 'python', ''};
    sProcess.options.implementation.Type    = 'radio_linelabel';
    sProcess.options.implementation.Value   = 'matlab';
    sProcess.options.implementation.Controller.matlab = 'Matlab';
    sProcess.options.implementation.Controller.python = 'Python';
    % === FREQUENCY RANGE
    sProcess.options.freqrange.Comment = 'Frequency range for analysis: ';
    sProcess.options.freqrange.Type    = 'freqrange_static';   % 'freqrange'
    sProcess.options.freqrange.Value   = {[1 40], 'Hz', 1};
    % === POWER LINE
    sProcess.options.powerline.Comment = {'None', '50 Hz', '60 Hz', 'Ignore power line frequencies:'; '-5', '50', '60', ''};
    sProcess.options.powerline.Type    = 'radio_linelabel';
    sProcess.options.powerline.Value   = '60';
    sProcess.options.powerline.Class   = 'Matlab';
    % === PEAK TYPE
    sProcess.options.peaktype.Comment = {'Gaussian', 'Cauchy*', 'Best of both* (* experimental)', 'Peak model:'; 'gaussian', 'cauchy', 'best', ''};
    sProcess.options.peaktype.Type    = 'radio_linelabel';
    sProcess.options.peaktype.Value   = 'gaussian';
    sProcess.options.peaktype.Class   = 'Matlab';
    % === PEAK WIDTH LIMITS
    sProcess.options.peakwidth.Comment = 'Peak width limits (default=[0.5-12]): ';
    sProcess.options.peakwidth.Type    = 'freqrange_static';
    sProcess.options.peakwidth.Value   = {[0.5 12], 'Hz', 1};
    % === MAX PEAKS
    sProcess.options.maxpeaks.Comment = 'Maximum number of peaks (default=3): ';
    sProcess.options.maxpeaks.Type    = 'value';
    sProcess.options.maxpeaks.Value   = {3, '', 0};
    % === MEAN PEAK HEIGHT
    sProcess.options.minpeakheight.Comment = 'Minimum peak height (default=3): ';
    sProcess.options.minpeakheight.Type    = 'value';
    sProcess.options.minpeakheight.Value   = {3, 'dB', 1};
    % === PROXIMITY THRESHOLD
    sProcess.options.proxthresh.Comment = 'Proximity threshold (default=2): ';
    sProcess.options.proxthresh.Type    = 'value';
    sProcess.options.proxthresh.Value   = {2, 'stdev of peak model', 1};
    sProcess.options.proxthresh.Class   = 'Matlab';
    % === APERIODIC MODE
    sProcess.options.apermode.Comment = {'Fixed', 'Knee', 'Aperiodic mode (default=fixed):'; 'fixed', 'knee', ''};
    sProcess.options.apermode.Type    = 'radio_linelabel';
    sProcess.options.apermode.Value   = 'fixed';
    % === GUESS WEIGHT
    sProcess.options.guessweight.Comment = {'None', 'Weak', 'Strong', 'Guess weight (default=none):'; 'none', 'weak', 'strong', ''};
    sProcess.options.guessweight.Type    = 'radio_linelabel';
    sProcess.options.guessweight.Value   = 'none';
    sProcess.options.guessweight.Class   = 'Matlab';

    % === SORT PEAKS TYPE
    sProcess.options.sorttype.Comment = {'Peak parameters', 'Frequency bands', 'Sort peaks using:'; 'param', 'band', ''};
    sProcess.options.sorttype.Type    = 'radio_linelabel';
    sProcess.options.sorttype.Value   = 'param';
    sProcess.options.sorttype.Controller.param = 'Param';
    sProcess.options.sorttype.Controller.band = 'Band';
    sProcess.options.sorttype.Group   = 'output';
    % === SORT PEAKS PARAM
    sProcess.options.sortparam.Comment = {'Frequency', 'Amplitude', 'Std dev.', 'Sort by peak...'; 'frequency', 'amplitude', 'std', ''};
    sProcess.options.sortparam.Type    = 'radio_linelabel';
    sProcess.options.sortparam.Value   = 'frequency';
    sProcess.options.sortparam.Class   = 'Param';
    sProcess.options.sortparam.Group   = 'output';
    % === SORT FREQ BANDS
    DefaultFreqBands = bst_get('DefaultFreqBands');
    sProcess.options.sortbands.Comment = '';
    sProcess.options.sortbands.Type    = 'groupbands';
    sProcess.options.sortbands.Value   = DefaultFreqBands(:,1:2);
    sProcess.options.sortbands.Class   = 'Band';
    sProcess.options.sortbands.Group   = 'output';
end

%% ===== FORMAT COMMENT =====
function Comment = FormatComment(sProcess) %#ok<DEFNU>
    Comment = sProcess.Comment;
end

%% ===== RUN =====
function OutputFile = Run(sProcess, sInputs) %#ok<DEFNU>
    % Initialize returned list of files
    OutputFile = {};

    % Fetch user settings
    implementation = sProcess.options.implementation.Value;
    opt.freq_range          = sProcess.options.freqrange.Value{1};
    opt.peak_width_limits   = sProcess.options.peakwidth.Value{1};
    opt.max_peaks           = sProcess.options.maxpeaks.Value{1};
    opt.min_peak_height     = sProcess.options.minpeakheight.Value{1} / 10; % convert from dB to B
    opt.aperiodic_mode      = sProcess.options.apermode.Value;
    opt.peak_threshold      = 2;   % 2 std dev: parameter for interface simplification
    opt.border_threshold    = 1;   % 1 std dev: proximity to edge of spectrum, static in Python
    opt.return_spectrum     = 0;   % SPM/FT: set to 1
    % Matlab-only options
    opt.power_line          = sProcess.options.powerline.Value;
    opt.peak_type           = sProcess.options.peaktype.Value;
    opt.proximity_threshold = sProcess.options.proxthresh.Value{1};
    opt.guess_weight        = sProcess.options.guessweight.Value;
    opt.thresh_after        = true;   % Threshold after fitting always selected for Matlab (mirrors the Python FOOOF closest by removing peaks that do not satisfy a user's predetermined conditions)
    % Python-only options
    opt.verbose    = false;
    % Output options
    opt.sort_type  = sProcess.options.sorttype.Value;
    opt.sort_param = sProcess.options.sortparam.Value;
	opt.sort_bands = sProcess.options.sortbands.Value;

    % Check input frequency bounds
    if (any(opt.freq_range < 0) || opt.freq_range(1) >= opt.freq_range(2))
        bst_report('error','Invalid Frequency range');
        return
    end

    hasOptimTools = 0;
    if exist('fmincon') == 2 && strcmp(implementation,'matlab')
        hasOptimTools = 1;
        disp('Using constrained optimization, Guess Weight ignored.')
    end

    % Initialize returned list of files
    OutputFile = {};
    for iFile = 1:length(sInputs)
        bst_progress('text',['Standby: FOOOFing spectrum ' num2str(iFile) ' of ' num2str(length(sInputs))]);
        % Load input file
        PsdMat = in_bst_timefreq(sInputs(iFile).FileName);
        % Exclude 0Hz from the computation
        if (opt.freq_range(1) == 0) && (PsdMat.Freqs(1) == 0) && (length(PsdMat.Freqs) >= 2)
            opt.freq_range(1) = PsdMat.Freqs(2);
        end

        % === COMPUTE FOOOF MODEL ===
        % Switch between implementations
        switch (implementation)
            case 'matlab'   % Matlab standalone FOOOF
                [FOOOF_freqs, FOOOF_data] = FOOOF_matlab(PsdMat.TF, PsdMat.Freqs, opt, hasOptimTools);
            case 'python'
                opt.peak_type = 'gaussian';
                [FOOOF_freqs, FOOOF_data] = process_fooof_py('FOOOF_python', PsdMat.TF, PsdMat.Freqs, opt);
                % Remove unnecessary structure level, allowing easy concatenation across channels, e.g. for display.
                FOOOF_data = FOOOF_data.FOOOF;
            otherwise
                error('Invalid implentation.');
        end

        % === FOOOF ANALYSIS ===
        TFfooof = PsdMat.TF(:,1,ismember(PsdMat.Freqs,FOOOF_freqs));
        [ePeaks, eAperiodics, eStats] = FOOOF_analysis(FOOOF_data, PsdMat.RowNames, TFfooof, opt.max_peaks, opt.sort_type, opt.sort_param, opt.sort_bands);

        % === PREPARE OUTPUT STRUCTURE ===
        % Create file structure
        PsdMat.Options.FOOOF = struct(...
            'options',    opt, ...
            'freqs',      FOOOF_freqs, ...
            'data',       FOOOF_data, ...
            'peaks',      ePeaks, ...
            'aperiodics', eAperiodics, ...
            'stats',      eStats);
        % Comment: Add FOOOF
        if ~isempty(strfind(PsdMat.Comment, 'PSD:'))
            PsdMat.Comment = strrep(PsdMat.Comment, 'PSD:', 'specparam:');
        else
            PsdMat.Comment = strcat(PsdMat.Comment, ' | specparam');
        end
        % History: Computation
        PsdMat = bst_history('add', PsdMat, 'compute', 'specparam');

        % === SAVE FILE ===
        % Filename: add _fooof tag
        [fPath, fName, fExt] = bst_fileparts(file_fullpath(sInputs(iFile).FileName));
        NewFile = file_unique(bst_fullfile(fPath, [fName, '_specparam', fExt]));
        % Save file
        bst_save(NewFile, PsdMat, 'v6');
        % Add file to database structure
        db_add_data(sInputs(iFile).iStudy, NewFile, PsdMat);
        % Return new file
        OutputFile{end+1} = NewFile;
    end
end

%% ===================================================================================
%  ===== MATLAB FOOOF ================================================================
%  ===================================================================================

%% ===== MATLAB STANDALONE FOOOF =====
function [fs, fg] = FOOOF_matlab(TF, Freqs, opt, hOT)
    % Find all frequency values within user limits
    fMask = (round(Freqs.*10)./10 >= opt.freq_range(1)) & (round(Freqs.*10)./10 <= opt.freq_range(2)) & ~mod(sum(abs(round(Freqs.*10)./10-[1;2;3].*str2double(opt.power_line)) >= 2),3);
    fs = Freqs(fMask);
    spec = log10(squeeze(TF(:,1,fMask))); % extract log spectra
    nChan = size(TF,1);
    if nChan == 1, spec = spec'; end
    % Initalize FOOOF structs
    fg(nChan) = struct(...
            'aperiodic_params', [],...
            'peak_params',      [],...
            'peak_types',       '',...
            'ap_fit',           [],...
            'fooofed_spectrum', [],...
            'peak_fit',         [],...
            'error',            [],...
            'r_squared',        []);
    % Iterate across channels
    for chan = 1:nChan
        bst_progress('set', round(chan./nChan.*100));
        % Fit aperiodic
        aperiodic_pars = robust_ap_fit(fs, spec(chan,:), opt.aperiodic_mode);
        % Remove aperiodic
        flat_spec = flatten_spectrum(fs, spec(chan,:), aperiodic_pars, opt.aperiodic_mode);
        % Fit peaks
        [peak_pars, peak_function] = fit_peaks(fs, flat_spec, opt.max_peaks, opt.peak_threshold, opt.min_peak_height, ...
            opt.peak_width_limits/2, opt.proximity_threshold, opt.border_threshold, opt.peak_type, opt.guess_weight,hOT);
        if opt.thresh_after && ~hOT  % Check thresholding requirements are met for unbounded optimization
            peak_pars(peak_pars(:,2) < opt.min_peak_height,:)     = []; % remove peaks shorter than limit
            peak_pars(peak_pars(:,3) < opt.peak_width_limits(1)/2,:)  = []; % remove peaks narrower than limit
            peak_pars(peak_pars(:,3) > opt.peak_width_limits(2)/2,:)  = []; % remove peaks broader than limit
            peak_pars = drop_peak_cf(peak_pars, opt.border_threshold, opt.freq_range); % remove peaks outside frequency limits
            peak_pars(peak_pars(:,1) < 0,:) = []; % remove peaks with a centre frequency less than zero (bypass drop_peak_cf)
            peak_pars = drop_peak_overlap(peak_pars, opt.proximity_threshold); % remove smallest of two peaks fit too closely
        end
        % Refit aperiodic
        aperiodic = spec(chan,:);
        for peak = 1:size(peak_pars,1)
            aperiodic = aperiodic - peak_function(fs,peak_pars(peak,1), peak_pars(peak,2), peak_pars(peak,3));
        end
        aperiodic_pars = simple_ap_fit(fs, aperiodic, opt.aperiodic_mode);
        % Generate model fit
        ap_fit = gen_aperiodic(fs, aperiodic_pars, opt.aperiodic_mode);
        model_fit = ap_fit;
        for peak = 1:size(peak_pars,1)
            model_fit = model_fit + peak_function(fs,peak_pars(peak,1),...
                peak_pars(peak,2),peak_pars(peak,3));
        end
        % Calculate model error
        MSE = sum((spec(chan,:) - model_fit).^2)/length(model_fit);
        rsq_tmp = corrcoef(spec(chan,:),model_fit).^2;
        % Return FOOOF results
        aperiodic_pars(2) = abs(aperiodic_pars(2));
        fg(chan).aperiodic_params = aperiodic_pars;
        fg(chan).peak_params      = peak_pars;
        fg(chan).peak_types       = func2str(peak_function);
        fg(chan).ap_fit           = 10.^ap_fit;
        fg(chan).fooofed_spectrum = 10.^model_fit;
        fg(chan).peak_fit         = 10.^(model_fit-ap_fit);
        fg(chan).error            = MSE;
        fg(chan).r_squared        = rsq_tmp(2);
        if opt.return_spectrum
            fg(chan).power_spectrum = spec(chan,:);
        end
        %plot(fs', [fg(chan).ap_fit', fg(chan).peak_fit', fg(chan).fooofed_spectrum'])
    end
end

%% ===== GENERATE APERIODIC =====
function ap_vals = gen_aperiodic(freqs,aperiodic_params,aperiodic_mode)
%       Generate aperiodic values, from parameter definition.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%       	Frequency vector to create aperiodic component for.
%       aperiodic_params : 1x3 array
%           Parameters that define the aperiodic component.
%       aperiodic_mode : {'fixed', 'knee'}
%           Defines absence or presence of knee in aperiodic component.
%
%       Returns
%       -------
%       ap_vals : 1d array
%           Generated aperiodic values.

    switch aperiodic_mode
        case 'fixed'  % no knee
            ap_vals = expo_nk_function(freqs,aperiodic_params);
        case 'knee'
            ap_vals = expo_function(freqs,aperiodic_params);
        case 'floor'
            ap_vals = expo_fl_function(freqs,aperiodic_params);
    end
end

%% ===== CORE MODELS =====
function ys = gaussian(freqs, mu, hgt, sigma)
%       Gaussian function to use for fitting.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency vector to create gaussian fit for.
%       mu, hgt, sigma : doubles
%           Parameters that define gaussian function (centre frequency,
%           height, and standard deviation).
%
%       Returns
%       -------
%       ys :    1xn array
%       Output values for gaussian function.

    ys = hgt*exp(-(((freqs-mu)./sigma).^2) /2);

end

function ys = cauchy(freqs, ctr, hgt, gam)
%       Cauchy function to use for fitting.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency vector to create cauchy fit for.
%       ctr, hgt, gam : doubles
%           Parameters that define cauchy function (centre frequency,
%           height, and "standard deviation" [gamma]).
%
%       Returns
%       -------
%       ys :    1xn array
%       Output values for cauchy function.

    ys = hgt./(1+((freqs-ctr)/gam).^2);

end

function ys = expo_function(freqs,params)
%       Exponential function to use for fitting 1/f, with a 'knee' (maximum at low frequencies).
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Input x-axis values.
%       params : 1x3 array (offset, knee, exp)
%           Parameters (offset, knee, exp) that define Lorentzian function:
%           y = 10^offset * (1/(knee + x^exp))
%
%       Returns
%       -------
%       ys :    1xn array
%           Output values for exponential function.

    ys = params(1) - log10(abs(params(2)) +freqs.^params(3));

end

function ys = expo_nk_function(freqs, params)
%       Exponential function to use for fitting 1/f, without a 'knee'.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Input x-axis values.
%       params : 1x2 array (offset, exp)
%           Parameters (offset, exp) that define Lorentzian function:
%           y = 10^offset * (1/(x^exp))
%
%       Returns
%       -------
%       ys :    1xn array
%           Output values for exponential (no-knee) function.

    ys = params(1) - log10(freqs.^params(2));

end

function ys = expo_fl_function(freqs, params)

    ys = log10(f.^(params(1)) * 10^(params(2)) + params(3));

end

%% ===== FITTING ALGORITHM =====
function aperiodic_params = simple_ap_fit(freqs, power_spectrum, aperiodic_mode)
%       Fit the aperiodic component of the power spectrum.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency values for the power spectrum, in linear scale.
%       power_spectrum : 1xn array
%           Power values, in log10 scale.
%       aperiodic_mode : {'fixed','knee'}
%           Defines absence or presence of knee in aperiodic component.
%
%       Returns
%       -------
%       aperiodic_params : 1xn array
%           Parameter estimates for aperiodic fit.

%       Set guess params for lorentzian aperiodic fit, guess params set at init
    options = optimset('Display', 'off', 'TolX', 1e-4, 'TolFun', 1e-6, ...
        'MaxFunEvals', 5000, 'MaxIter', 5000);

    switch (aperiodic_mode)
        case 'fixed'  % no knee
            exp_guess = -(power_spectrum(end)-power_spectrum(1))./log10(freqs(end)./freqs(1));
            guess_vec = [power_spectrum(1), exp_guess];
            aperiodic_params = fminsearch(@error_expo_nk_function, guess_vec, options, freqs, power_spectrum);
        case 'knee'
            exp_guess = -(power_spectrum(end)-power_spectrum(1))./log10(freqs(end)./freqs(1));
            guess_vec = [power_spectrum(1),0, exp_guess];
            aperiodic_params = fminsearch(@error_expo_function, guess_vec, options, freqs, power_spectrum);
    end

end

function aperiodic_params = robust_ap_fit(freqs, power_spectrum, aperiodic_mode)
%       Fit the aperiodic component of the power spectrum robustly, ignoring outliers.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency values for the power spectrum, in linear scale.
%       power_spectrum : 1xn array
%           Power values, in log10 scale.
%       aperiodic_mode : {'fixed','knee'}
%           Defines absence or presence of knee in aperiodic component.
%
%       Returns
%       -------
%       aperiodic_params : 1xn array
%           Parameter estimates for aperiodic fit.

    % Do a quick, initial aperiodic fit
    popt = simple_ap_fit(freqs, power_spectrum, aperiodic_mode);
    initial_fit = gen_aperiodic(freqs, popt, aperiodic_mode);

    % Flatten power_spectrum based on initial aperiodic fit
    flatspec = power_spectrum - initial_fit;

    % Flatten outliers - any points that drop below 0
    flatspec(flatspec(:) < 0) = 0;

    % Use percential threshold, in terms of # of points, to extract and re-fit
    perc_thresh = bst_prctile(flatspec, 0.025);
    perc_mask = flatspec <= perc_thresh;
    freqs_ignore = freqs(perc_mask);
    spectrum_ignore = power_spectrum(perc_mask);

    % Second aperiodic fit - using results of first fit as guess parameters

    options = optimset('Display', 'off', 'TolX', 1e-4, 'TolFun', 1e-6, ...
        'MaxFunEvals', 5000, 'MaxIter', 5000);
    guess_vec = popt;

    switch (aperiodic_mode)
        case 'fixed'  % no knee
            aperiodic_params = fminsearch(@error_expo_nk_function, guess_vec, options, freqs_ignore, spectrum_ignore);
        case 'knee'
            aperiodic_params = fminsearch(@error_expo_function, guess_vec, options, freqs_ignore, spectrum_ignore);
    end
end

function spectrum_flat = flatten_spectrum(freqs, power_spectrum, robust_aperiodic_params, aperiodic_mode)
%       Flatten the power spectrum by removing the aperiodic component.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency values for the power spectrum, in linear scale.
%       power_spectrum : 1xn array
%           Power values, in log10 scale.
%       robust_aperiodic_params : 1x2 or 1x3 array (see aperiodic_mode)
%           Parameter estimates for aperiodic fit.
%       aperiodic_mode : 1 or 2
%           Defines absence or presence of knee in aperiodic component.
%
%       Returns
%       -------
%       spectrum_flat : 1xn array
%           Flattened (aperiodic removed) power spectrum.

spectrum_flat = power_spectrum - gen_aperiodic(freqs,robust_aperiodic_params,aperiodic_mode);

end

function [model_params,peak_function] = fit_peaks(freqs, flat_iter, max_n_peaks, peak_threshold, min_peak_height, gauss_std_limits, proxThresh, bordThresh, peakType, guess_weight,hOT)
%       Iteratively fit peaks to flattened spectrum.
%
%       Parameters
%       ----------
%       freqs : 1xn array
%           Frequency values for the power spectrum, in linear scale.
%       flat_iter : 1xn array
%           Flattened (aperiodic removed) power spectrum.
%       max_n_peaks : double
%           Maximum number of gaussians to fit within the spectrum.
%       peak_threshold : double
%           Threshold (in standard deviations of noise floor) to detect a peak.
%       min_peak_height : double
%           Minimum height of a peak (in log10).
%       gauss_std_limits : 1x2 double
%           Limits to gaussian (cauchy) standard deviation (gamma) when detecting a peak.
%       proxThresh : double
%           Minimum distance between two peaks, in st. dev. (gamma) of peaks.
%       peakType : {'gaussian', 'cauchy', 'both'}
%           Which types of peaks are being fitted
%       guess_weight : {'none', 'weak', 'strong'}
%           Parameter to weigh initial estimates during optimization (None, Weak, or Strong)
%       hOT : 0 or 1
%           Defines whether to use constrained optimization, fmincon, or
%           basic simplex, fminsearch.
%
%       Returns
%       -------
%       gaussian_params : mx3 array, where m = No. of peaks.
%           Parameters that define the peak fit(s). Each row is a peak, as [mean, height, st. dev. (gamma)].
    switch peakType
        case 'gaussian' % gaussian only
            peak_function = @gaussian; % Identify peaks as gaussian
            % Initialize matrix of guess parameters for gaussian fitting.
            guess_params = zeros(max_n_peaks, 3);
            % Save intact flat_spectrum
            flat_spec = flat_iter;
            % Find peak: Loop through, finding a candidate peak, and fitting with a guess gaussian.
            % Stopping procedure based on either the limit on # of peaks,
            % or the relative or absolute height thresholds.
            for guess = 1:max_n_peaks
                % Find candidate peak - the maximum point of the flattened spectrum.
                max_ind = find(flat_iter == max(flat_iter));
                max_height = flat_iter(max_ind);

                % Stop searching for peaks once max_height drops below height threshold.
                if max_height <= peak_threshold * std(flat_iter)
                    break
                end

                % Set the guess parameters for gaussian fitting - mean and height.
                guess_freq = freqs(max_ind);
                guess_height = max_height;

                % Halt fitting process if candidate peak drops below minimum height.
                if guess_height <= min_peak_height
                    break
                end

                % Data-driven first guess at standard deviation
                % Find half height index on each side of the center frequency.
                half_height = 0.5 * max_height;

                le_ind = sum(flat_iter(1:max_ind) <= half_height);
                ri_ind = length(flat_iter) - sum(flat_iter(max_ind:end) <= half_height)+1;

                % Keep bandwidth estimation from the shortest side.
                % We grab shortest to avoid estimating very large std from overalapping peaks.
                % Grab the shortest side, ignoring a side if the half max was not found.
                % Note: will fail if both le & ri ind's end up as None (probably shouldn't happen).
                short_side = min(abs([le_ind,ri_ind]-max_ind));

                % Estimate std from FWHM. Calculate FWHM, converting to Hz, get guess std from FWHM
                fwhm = short_side * 2 * (freqs(2)-freqs(1));
                guess_std = fwhm / (2 * sqrt(2 * log(2)));

                % Check that guess std isn't outside preset std limits; restrict if so.
                % Note: without this, curve_fitting fails if given guess > or < bounds.
                if guess_std < gauss_std_limits(1)
                    guess_std = gauss_std_limits(1);
                end
                if guess_std > gauss_std_limits(2)
                    guess_std = gauss_std_limits(2);
                end

                % Collect guess parameters.
                guess_params(guess,:) = [guess_freq, guess_height, guess_std];

                % Subtract best-guess gaussian.
                peak_gauss = gaussian(freqs, guess_freq, guess_height, guess_std);
                flat_iter = flat_iter - peak_gauss;

            end
            % Remove unused guesses
            guess_params(guess_params(:,1) == 0,:) = [];

            % Check peaks based on edges, and on overlap
            % Drop any that violate requirements.
            guess_params = drop_peak_cf(guess_params, bordThresh, [min(freqs) max(freqs)]);
            guess_params = drop_peak_overlap(guess_params, proxThresh);

            % If there are peak guesses, fit the peaks, and sort results.
            if ~isempty(guess_params)
                model_params = fit_peak_guess(guess_params, freqs, flat_spec, 1, guess_weight, gauss_std_limits,hOT);
            else
                model_params = zeros(1, 3);
            end

        case 'cauchy' % cauchy only
            peak_function = @cauchy; % Identify peaks as cauchy
            guess_params = zeros(max_n_peaks, 3);
            flat_spec = flat_iter;
            for guess = 1:max_n_peaks
                max_ind = find(flat_iter == max(flat_iter));
                max_height = flat_iter(max_ind);
                if max_height <= peak_threshold * std(flat_iter)
                    break
                end
                guess_freq = freqs(max_ind);
                guess_height = max_height;
                if guess_height <= min_peak_height
                    break
                end
                half_height = 0.5 * max_height;
                le_ind = sum(flat_iter(1:max_ind) <= half_height);
                ri_ind = length(flat_iter) - sum(flat_iter(max_ind:end) <= half_height);
                short_side = min(abs([le_ind,ri_ind]-max_ind));

                % Estimate gamma from FWHM. Calculate FWHM, converting to Hz, get guess gamma from FWHM
                fwhm = short_side * 2 * (freqs(2)-freqs(1));
                guess_gamma = fwhm/2;
                % Check that guess gamma isn't outside preset limits; restrict if so.
                % Note: without this, curve_fitting fails if given guess > or < bounds.
                if guess_gamma < gauss_std_limits(1)
                    guess_gamma = gauss_std_limits(1);
                end
                if guess_gamma > gauss_std_limits(2)
                    guess_gamma = gauss_std_limits(2);
                end

                % Collect guess parameters.
                guess_params(guess,:) = [guess_freq(1), guess_height, guess_gamma];

                % Subtract best-guess cauchy.
                peak_cauchy = cauchy(freqs, guess_freq(1), guess_height, guess_gamma);
                flat_iter = flat_iter - peak_cauchy;

            end
            guess_params(guess_params(:,1) == 0,:) = [];
            guess_params = drop_peak_cf(guess_params, bordThresh, [min(freqs) max(freqs)]);
            guess_params = drop_peak_overlap(guess_params, proxThresh);

            % If there are peak guesses, fit the peaks, and sort results.
            if ~isempty(guess_params)
                model_params = fit_peak_guess(guess_params, freqs, flat_spec, 2, guess_weight, gauss_std_limits,hOT);
            else
                model_params = zeros(1, 3);
            end
        case 'best' % best of both: model both fits and compare error, save best
            % Gaussian Fit
            guess_params = zeros(max_n_peaks, 3);
            flat_spec = flat_iter;
            for guess = 1:max_n_peaks
                max_ind = find(flat_iter == max(flat_iter));
                max_height = flat_iter(max_ind);
                if max_height <= peak_threshold * std(flat_iter)
                    break
                end
                guess_freq = freqs(max_ind);
                guess_height = max_height;
                if guess_height <= min_peak_height
                    break
                end
                half_height = 0.5 * max_height;
                le_ind = sum(flat_iter(1:max_ind) <= half_height);
                ri_ind = length(flat_iter) - sum(flat_iter(max_ind:end) <= half_height)+1;
                short_side = min(abs([le_ind,ri_ind]-max_ind));
                fwhm = short_side * 2 * (freqs(2)-freqs(1));
                guess_std = fwhm / (2 * sqrt(2 * log(2)));
                if guess_std < gauss_std_limits(1)
                    guess_std = gauss_std_limits(1);
                end
                if guess_std > gauss_std_limits(2)
                    guess_std = gauss_std_limits(2);
                end
                guess_params(guess,:) = [guess_freq, guess_height, guess_std];
                peak_gauss = gaussian(freqs, guess_freq, guess_height, guess_std);
                flat_iter = flat_iter - peak_gauss;
            end
            guess_params(guess_params(:,1) == 0,:) = [];
            guess_params = drop_peak_cf(guess_params, bordThresh, [min(freqs) max(freqs)]);
            guess_params = drop_peak_overlap(guess_params, proxThresh);
            if ~isempty(guess_params)
                gauss_params = fit_peak_guess(guess_params, freqs, flat_spec, 1, guess_weight, gauss_std_limits,hOT);
                flat_gauss = zeros(size(freqs));
                for peak = 1:size(gauss_params,1)
                    flat_gauss =  flat_gauss + gaussian(freqs,gauss_params(peak,1),...
                        gauss_params(peak,2),gauss_params(peak,3));
                end
                error_gauss = sum((flat_gauss-flat_spec).^2);
            else
                gauss_params = zeros(1, 3); error_gauss = 1E10;
            end

            % Cauchy Fit
            guess_params = zeros(max_n_peaks, 3);
            flat_iter = flat_spec;
            for guess = 1:max_n_peaks
                max_ind = find(flat_iter == max(flat_iter));
                max_height = flat_iter(max_ind);
                if max_height <= peak_threshold * std(flat_iter)
                    break
                end
                guess_freq = freqs(max_ind);
                guess_height = max_height;
                if guess_height <= min_peak_height
                    break
                end
                half_height = 0.5 * max_height;
                le_ind = sum(flat_iter(1:max_ind) <= half_height);
                ri_ind = length(flat_iter) - sum(flat_iter(max_ind:end) <= half_height)+1;
                short_side = min(abs([le_ind,ri_ind]-max_ind));
                fwhm = short_side * 2 * (freqs(2)-freqs(1));
                guess_gamma = fwhm/2;
                if guess_gamma < gauss_std_limits(1)
                    guess_gamma = gauss_std_limits(1);
                end
                if guess_gamma > gauss_std_limits(2)
                    guess_gamma = gauss_std_limits(2);
                end
                guess_params(guess,:) = [guess_freq(1), guess_height, guess_gamma];
                peak_cauchy = cauchy(freqs, guess_freq(1), guess_height, guess_gamma);
                flat_iter = flat_iter - peak_cauchy;
            end
            guess_params(guess_params(:,1) == 0,:) = [];
            guess_params = drop_peak_cf(guess_params, bordThresh, [min(freqs) max(freqs)]);
            guess_params = drop_peak_overlap(guess_params, proxThresh);
            if ~isempty(guess_params)
                cauchy_params = fit_peak_guess(guess_params, freqs, flat_spec, 2, guess_weight, gauss_std_limits,hOT);
                flat_cauchy = zeros(size(freqs));
                for peak = 1:size(cauchy_params,1)
                    flat_cauchy =  flat_cauchy + cauchy(freqs,cauchy_params(peak,1),...
                        cauchy_params(peak,2),cauchy_params(peak,3));
                end
                error_cauchy = sum((flat_cauchy-flat_spec).^2);
            else
                cauchy_params = zeros(1, 3); error_cauchy = 1E10;
            end
            % Save least-error model
                if min([error_gauss,error_cauchy]) == error_gauss
                    model_params = gauss_params;
                    peak_function = @gaussian;
                else
                    model_params = cauchy_params;
                    peak_function = @cauchy;
                end
    end

end

function guess = drop_peak_cf(guess, bw_std_edge, freq_range)
%       Check whether to drop peaks based on center's proximity to the edge of the spectrum.
%
%       Parameters
%       ----------
%       guess : mx3 array, where m = No. of peaks.
%           Guess parameters for peak fits.
%
%       Returns
%       -------
%       guess : qx3 where q <= m No. of peaks.
%           Guess parameters for peak fits.

    cf_params = guess(:,1)';
    bw_params = guess(:,3)' * bw_std_edge;

    % Check if peaks within drop threshold from the edge of the frequency range.

    keep_peak = abs(cf_params-freq_range(1)) > bw_params & ...
        abs(cf_params-freq_range(2)) > bw_params;

    % Drop peaks that fail the center frequency edge criterion
    guess = guess(keep_peak,:);

end

function guess = drop_peak_overlap(guess, proxThresh)
%       Checks whether to drop gaussians based on amount of overlap.
%
%       Parameters
%       ----------
%       guess : mx3 array, where m = No. of peaks.
%           Guess parameters for peak fits.
%       proxThresh: double
%           Proximity threshold (in st. dev. or gamma) between two peaks.
%
%       Returns
%       -------
%       guess : qx3 where q <= m No. of peaks.
%           Guess parameters for peak fits.
%
%       Note
%       -----
%       For any gaussians with an overlap that crosses the threshold,
%       the lowest height guess guassian is dropped.

    % Sort the peak guesses, so can check overlap of adjacent peaks
    guess = sortrows(guess);

    % Calculate standard deviation bounds for checking amount of overlap

    bounds = [guess(:,1) - guess(:,3) * proxThresh, ...
        guess(:,1), guess(:,1) + guess(:,3) * proxThresh];

    % Loop through peak bounds, comparing current bound to that of next peak
    drop_inds =  [];

    for ind = 1:size(bounds,1)-1

        b_0 = bounds(ind,:);
        b_1 = bounds(ind + 1,:);

        % Check if bound of current peak extends into next peak
        if b_0(2) > b_1(1)
            % If so, get the index of the gaussian with the lowest height (to drop)
            drop_inds = [drop_inds (ind - 1 + find(guess(ind:ind+1,2) == ...
                min(guess(ind,2),guess(ind+1,2))))];
        end
    end
    % Drop any peaks guesses that overlap too much, based on threshold.
    guess(drop_inds,:) = [];
end

function peak_params = fit_peak_guess(guess, freqs, flat_spec, peak_type, guess_weight, std_limits, hOT)
%     Fits a group of peak guesses with a fit function.
%
%     Parameters
%     ----------
%       guess : mx3 array, where m = No. of peaks.
%           Guess parameters for peak fits.
%       freqs : 1xn array
%           Frequency values for the power spectrum, in linear scale.
%       flat_iter : 1xn array
%           Flattened (aperiodic removed) power spectrum.
%       peakType : {'gaussian', 'cauchy', 'best'}
%           Which types of peaks are being fitted.
%       guess_weight : 'none', 'weak', 'strong'
%           Parameter to weigh initial estimates during optimization.
%       std_limits: 1x2 array
%           Minimum and maximum standard deviations for distribution.
%       hOT : 0 or 1
%           Defines whether to use constrained optimization, fmincon, or
%           basic simplex, fminsearch.
%
%       Returns
%       -------
%       peak_params : mx3, where m =  No. of peaks.
%           Peak parameters post-optimization.

    if hOT % Use OptimToolbox for fmincon
        options = optimset('Display', 'off', 'TolX', 1e-3, 'TolFun', 1e-5, ...
        'MaxFunEvals', 3000, 'MaxIter', 3000); % Tuned options
        lb = [guess(:,1)-guess(:,3)*2,zeros(size(guess(:,2))),ones(size(guess(:,3)))*std_limits(1)];
        ub = [guess(:,1)+guess(:,3)*2,inf(size(guess(:,2))),ones(size(guess(:,3)))*std_limits(2)];
        peak_params = fmincon(@error_model_constr,guess,[],[],[],[], ...
            lb,ub,[],options,freqs,flat_spec, peak_type);
    else % Use basic simplex approach, fminsearch, with guess_weight
        options = optimset('Display', 'off', 'TolX', 1e-4, 'TolFun', 1e-5, ...
        'MaxFunEvals', 5000, 'MaxIter', 5000);
        peak_params = fminsearch(@error_model,...
            guess, options, freqs, flat_spec, peak_type, guess, guess_weight);
    end
end

%% ===== ERROR FUNCTIONS =====
function err = error_expo_nk_function(params,xs,ys)
    ym = -log10(xs.^params(2)) + params(1);
    err = sum((ys - ym).^2);
end

function err = error_expo_function(params,xs,ys)
    ym = expo_function(xs,params);
    err = sum((ys - ym).^2);
end

function err = error_model(params, xVals, yVals, peak_type, guess, guess_weight)
    fitted_vals = 0;
    weak = 1E2;
    strong = 1E7;
    for set = 1:size(params,1)
        switch (peak_type)
            case 1 % Gaussian
                fitted_vals = fitted_vals + gaussian(xVals, params(set,1), params(set,2), params(set,3));
            case 2 % Cauchy
                fitted_vals = fitted_vals + cauchy(xVals, params(set,1), params(set,2), params(set,3));
        end
    end
    switch guess_weight
        case 'none'
            err = sum((yVals - fitted_vals).^2);
        case 'weak' % Add small weight to deviations from guess m and amp
            err = sum((yVals - fitted_vals).^2) + ...
                 weak*sum((params(:,1)-guess(:,1)).^2) + ...
                 weak*sum((params(:,2)-guess(:,2)).^2);
        case 'strong' % Add large weight to deviations from guess m and amp
            err = sum((yVals - fitted_vals).^2) + ...
                 strong*sum((params(:,1)-guess(:,1)).^2) + ...
                 strong*sum((params(:,2)-guess(:,2)).^2);
    end
end

function err = error_model_constr(params, xVals, yVals, peak_type)
    fitted_vals = 0;
    for set = 1:size(params,1)
        switch (peak_type)
            case 1 % Gaussian
                fitted_vals = fitted_vals + gaussian(xVals, params(set,1), params(set,2), params(set,3));
            case 2 % Cauchy
                fitted_vals = fitted_vals + cauchy(xVals, params(set,1), params(set,2), params(set,3));
        end
    end
    err = sum((yVals - fitted_vals).^2);
end

%% ===================================================================================
%  ===== FOOOF STATS =================================================================
%  ===================================================================================
function [ePeaks, eAper, eStats] = FOOOF_analysis(FOOOF_data, ChanNames, TF, max_peaks, sort_type, sort_param, sort_bands)
    % ===== EXTRACT PEAKS =====
    % Organize/extract peak components from FOOOF models
    nChan = numel(ChanNames);
    maxEnt = nChan * max_peaks;
    switch sort_type
        case 'param'
            % Initialize output struct
            ePeaks = struct('channel', [], 'center_frequency', [],...
                'amplitude', [], 'std_dev', []);
            % Collect data from all peaks
            i = 0;
            for chan = 1:nChan
                if ~isempty(FOOOF_data(chan).peak_params)
                    for p = 1:size(FOOOF_data(chan).peak_params,1)
                        i = i +1;
                        ePeaks(i).channel = ChanNames(chan);
                        ePeaks(i).center_frequency = FOOOF_data(chan).peak_params(p,1);
                        ePeaks(i).amplitude = FOOOF_data(chan).peak_params(p,2);
                        ePeaks(i).std_dev = FOOOF_data(chan).peak_params(p,3);
                    end
                end
            end
            % Apply specified sort
            switch sort_param
                case 'frequency'
                    [tmp,iSort] = sort([ePeaks.center_frequency]);
                    ePeaks = ePeaks(iSort);
                case 'amplitude'
                    [tmp,iSort] = sort([ePeaks.amplitude]);
                    ePeaks = ePeaks(iSort(end:-1:1));
                case 'std'
                    [tmp,iSort] = sort([ePeaks.std_dev]);
                    ePeaks = ePeaks(iSort);
            end
        case 'band'
            % Initialize output struct
            ePeaks = struct('channel', [], 'center_frequency', [],...
                'amplitude', [], 'std_dev', [], 'band', []);
            % Generate bands from input
            bands = process_tf_bands('Eval', sort_bands);
            % Collect data from all peaks
            i = 0;
            for chan = 1:nChan
                if ~isempty(FOOOF_data(chan).peak_params)
                    for p = 1:size(FOOOF_data(chan).peak_params,1)
                        i = i +1;
                        ePeaks(i).channel = ChanNames(chan);
                        ePeaks(i).center_frequency = FOOOF_data(chan).peak_params(p,1);
                        ePeaks(i).amplitude = FOOOF_data(chan).peak_params(p,2);
                        ePeaks(i).std_dev = FOOOF_data(chan).peak_params(p,3);
                        % Find name of frequency band from user definitions
                        bandRanges = cell2mat(bands(:,2));
                        iBand = find(ePeaks(i).center_frequency >= bandRanges(:,1) & ePeaks(i).center_frequency <= bandRanges(:,2));
                        if ~isempty(iBand)
                            ePeaks(i).band = bands{iBand,1};
                        else
                            ePeaks(i).band = 'None';
                        end
                    end
                end
            end
    end

    % ===== EXTRACT APERIODIC =====
    % Organize/extract aperiodic components from FOOOF models
    hasKnee = length(FOOOF_data(1).aperiodic_params) - 2;
    % Initialize output struct
    eAper = struct('channel', [], 'offset', [], 'exponent', []);
    for chan = 1:nChan
            eAper(chan).channel = ChanNames(chan);
            eAper(chan).offset = FOOOF_data(chan).aperiodic_params(1);
        if hasKnee % Legacy FOOOF alters order of parameters
            eAper(chan).exponent = FOOOF_data(chan).aperiodic_params(3);
            eAper(chan).knee_frequency = FOOOF_data(chan).aperiodic_params(2);
        else
            eAper(chan).exponent = FOOOF_data(chan).aperiodic_params(2);
        end
    end

    % ===== EXTRACT STAT =====
    % Organize/extract stats from FOOOF models
    % Initialize output struct
    eStats = struct('channel', ChanNames);
    for chan = 1:nChan
        eStats(chan).MSE = FOOOF_data(chan).error;
        eStats(chan).r_squared = FOOOF_data(chan).r_squared;
        spec = squeeze(log10(TF(chan,:,:)));
        fspec = squeeze(log10(FOOOF_data(chan).fooofed_spectrum))';
        eStats(chan).frequency_wise_error = abs(spec-fspec);
    end
end

author: brainstorm-tools