text:
function [E, wE] = nfgCompareVolumeTrueError(fgGS, fgSS, g_radius,vR,phantomDir)
%Compare volume of two fiber groups using trueError program
%
%   [E, wE] = nfgCompareVolumeTrueError(fgGS, fgSS, g_radius,vR,phantomDir)
%
% Assumes that fgGS is the gold standard group and fgSS is the selected
% fibers that match the gold standard from the test group.  g_radius is the
% radius of the gold standard fibers and vR is a vector of radii that might
% be the thickness of the test fibers.
%
%
% NOTES:

% Let's see how trueError compares these fibers
tempFiberDir = nfgGetName('tempFiberDir',phantomDir);
ctrparamsFile = nfgGetName('ctrparamsFile',phantomDir);
noisyImg = nfgGetName('noisyImg',phantomDir);
bvalsFile = nfgGetName('bvalsFile',phantomDir);
bvecsFile = nfgGetName('bvecsFile',phantomDir);
b0File = nfgGetName('b0File',phantomDir);
tensorsFile = nfgGetName('tensorsFile',phantomDir);
wmROIFile = nfgGetName('wmROIFile',phantomDir);
gmROIFile = nfgGetName('gmROIFile',phantomDir);

% XXX special check to find fibers less than 5 points in gold standard
fiberLen = cellfun('size',fgGS.fibers,2);
lt5 = fiberLen<5;
if sum(lt5)>0
    disp(['Warning there are ' num2str(sum(lt5)) ' fibers with less than 5 pts.']);
    fgGS.fibers = fgGS.fibers(~lt5);
end

mkdir(tempFiberDir);
mtrExportFibers(fgSS,fullfile(tempFiberDir,'ss.pdb'));
mtrExportFibers(fgGS,fullfile(tempFiberDir,'gs.pdb'));

% Convert to SBfloat
argParams = [' -i ' ctrparamsFile];
argD_SS = [' -p ' fullfile(tempFiberDir,'ss_0.SBfloat')];
argInput = [' ' fullfile(tempFiberDir,'ss.pdb')];
argThresh = [' --thresh ' num2str(length(fgSS.fibers)) ' --seq '];
cmd = ['contrack_score' argParams argD_SS argThresh argInput];
%disp(cmd);
[s,r] = system(cmd);
argD_GS = [' -p ' fullfile(tempFiberDir,'gs_0.SBfloat')];
argInput = [' ' fullfile(tempFiberDir,'gs.pdb')];
argThresh = [' --thresh ' num2str(length(fgGS.fibers)) ' --seq '];
cmd = ['contrack_score' argParams argD_GS argThresh argInput];
%disp(cmd);
[s,r] = system(cmd);
% Run true error on gold
argD_GS = [' -d ' fullfile(tempFiberDir,'gs_0.SBfloat')];
argR = [' -r ' noisyImg];
argVal = [' --val ' bvalsFile];
argVec = [' --vec ' bvecsFile];
arg0 = [' -0 ' b0File];
argG = [' -g ' gmROIFile];
argM = [' -m ' wmROIFile];
argT = [' --ten ' tensorsFile];
argGroupSize = ' -v 2';
b0 = niftiRead(b0File);
argSubSize = [' -s 0,' num2str(b0.dim(1)-1) ',0,' num2str(b0.dim(2)-1) ',0,' num2str(b0.dim(3)-1)];
argW = ' -w 0';
argDiameter = [' --diameter ' num2str(max(g_radius)*2)];
fracGoldFile = fullfile(tempFiberDir,'fgs.nii.gz');
argFractionFile = [' --fraction ' fracGoldFile];
cmd = ['trueError' argR argD_GS argVal argVec arg0 argG argM argT argGroupSize argSubSize argW argDiameter argFractionFile];
%disp(cmd);
[s,r] = system(cmd);
fracG = niftiRead(fracGoldFile);

% Search for the right radius to compare the selection to the gold
argD_SS = [' -d ' fullfile(tempFiberDir,'ss_0.SBfloat')];
fracImgs = zeros(size(repmat(b0.data,[1 1 1 length(vR)])));
E = zeros(size(vR));
for rr=1:length(vR)
    strD = num2str(vR(rr)*2);
    disp(['Comparing trueError volume using test diameter ' strD ' ...']);
    argDiameter = [' --diameter ' strD];
    fracFile = fullfile(tempFiberDir,['fss' strD '.nii.gz']);
    argFractionFile = [' --fraction ' fracFile];
    cmd = ['trueError' argR argD_SS argVal argVec arg0 argG argM argT argGroupSize argSubSize argW argDiameter argFractionFile];
    %disp(cmd);
    [s,r] = system(cmd);
    ni = niftiRead(fracFile);
    fracImgs(:,:,:,rr) = ni.data;
    %E(rr) = sum(abs(ni.data(:)-fracG.data(:)))/sum(fracG.data(:)) * 100;
    E(rr) = abs((sum(ni.data(:))-sum(fracG.data(:)))/sum(fracG.data(:)) * 100);
end
wE = sum(fracG.data(:));
% Cleanup the temp space
rmdir(tempFiberDir,'s');

return;


text:
function [pd, labels] = imPerimeterDensity(img, varargin)
% Perimeter density of a 2D binary structure, using Crofton formula.
%
%   Pv = imPerimeterDensity(IMG)
%
%   Example
%   imPerimeterDensity
%
%   See also
%     imPerimeter

% ------
% Author: David Legland
% e-mail: david.legland@inrae.fr
% Created: 2010-01-21,    using Matlab 7.9.0.529 (R2009b)
% Copyright 2010 INRA - Cepia Software Platform.

%% Pre-processing

% check image dimension
if ndims(img) ~= 2 %#ok<ISMAT>
    error('first argument should be a 2D image');
end

% in case of a label image, return a vector with a set of results
if ~islogical(img)
    % extract labels (considers 0 as background)
    labels = unique(img);
    labels(labels==0) = [];

    % allocate result array
    nLabels = length(labels);
    pd = zeros(nLabels, 1);

    props = regionprops(img, 'BoundingBox');

    % compute perimeter of each label considered as binary image
    for i = 1:nLabels
        label = labels(i);
        bin = imcrop(img, props(label).BoundingBox) == label;
        pd(i) = imPerimeterDensity(bin, varargin{:});
    end

    return;
end

%% Extract input arguments

% in case of binary image, compute only one label
labels = 1;

% default number of directions
nDirs = 4;

% default image resolution
delta = [1 1];

% parse parameter name-value pairs
while ~isempty(varargin)
    var = varargin{1};

    if isnumeric(var)
        % option is either number of directions or resolution
        if isscalar(var)
            nDirs = var;
        else
            delta = var;
        end
        varargin(1) = [];

    elseif ischar(var)
        if length(varargin) < 2
            error('Parameter name must be followed by parameter value');
        end

        if strcmpi(var, 'ndirs')
            nDirs = varargin{2};
        elseif strcmpi(var, 'resolution')
            delta = varargin{2};
        else
            error(['Unknown parameter name: ' var]);
        end

        varargin(1:2) = [];
    end
end

%% Compute perimeter within image, and normalize by area

% component area in image
p = imPerimeterEstimate(img, nDirs, delta);

% total area of image, without borders
totalArea = prod(size(img)-1) * prod(delta);

% compute perimeter density
pd = p / totalArea;

author: mattools

repo: matImage

save_path: github-repos/MATLAB/mattools-matImage


text:
function plotVR=SCARA_VR_PLOT(T1,T2,T4,d4b,dist)
SCARA30 = vrworld('SCARA30.wrl');
open(SCARA30)
    radian=T1*pi/180;
    SCARA30.a1.rotation = [1, 0, 0, radian];

    radian=-T2*pi/180;
    SCARA30.a2.rotation = [1, 0, 0, radian];

    radian=T4*pi/180;
    SCARA30.d4b.rotation = [0, 1, 0, radian];

    SCARA30.d4b.translation = [0, d4b, 0];

EndV1=SCARA30.EndV1.translation;
EndV2=SCARA30.EndV2.translation;
x1=EndV1(1);
x2=EndV2(1);
% y1=EndV1(2);
% y2=EndV2(2);
% sighn=sign([y1 y2]);
% y1f=.5*sighn(1)*dist
% y2f=.5*sighn(2)*dist
y1=-.5*dist;
y2=.5*dist;

SCARA30.EndV1.translation = [x1, y1, 0];
SCARA30.EndV2.translation = [x2, y2, 0];

author: Sable

repo: mcbench-benchmarks


text:
function [cnmfAnalysisOutput] = computeCnmfSignalExtractionPatch(inputMovie,numExpectedComponents,varargin)
	% Brapper function for CNMF, update for most recent versions.
	% Building off of demo_script.m in CNMF github repo
	% Most recent commit tested on: https://github.com/epnev/ca_source_extraction/commit/187bbdbe66bca466b83b81861b5601891a95b8d1
	% https://github.com/epnev/ca_source_extraction/blob/master/demo_script_class.m
	% Biafra Ahanonu
	% started: 2019.03.11
	% inputs
		% inputMovie - a string or a cell array of strings pointing to the movies to be analyzed (recommended). Else, [x y t] matrix where t = frames.
		% numExpectedComponents - number of expected components
	% outputs
		% cnmfAnalysisOutput - structure containing extractedImages and extractedSignals along with input parameters to the algorithm
	% READ BEFORE RUNNING
		% Get CVX from http://cvxr.com/cvx/doc/install.html
		% Run the below commands in Matlab after unzipping
		% cvx_setup
		% cvx_save_prefs (permanently stores settings)

	[cnmfAnalysisOutput] = ciapkg.signal_extraction.computeCnmfSignalExtractionPatch(inputMovie,numExpectedComponents,'passArgs', varargin);
end

author: bahanonu

repo: ciatah

save_path: github-repos/MATLAB/bahanonu-ciatah


text:
% example_radar.m
% ====================================================>
% The following seem to be descent parameters:
%   - ProbOfDetection = 0.8
%   - ProbOfDeath = 0.2
%   - VelocityErrVariance = 2
%   - ObsErrVariance = 50
%   - ProbOfConfirm = 0.9
%   - PHD.BirthScheme = {'Expansion', 5000}
%

%% Radar specific settings
RadarName = 'staddon';
switch(RadarName)
    case('staddon')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    -212.289393440959;... % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
        RadarCoords.lat = 50.33933333;
        RadarCoords.lon = -4.12527778;
%         RadarCoords.lat = 50.346069;
%         RadarCoords.lon = -4.113670;
    case('longroom')
        % Surveillance region parameters
        V_bounds = [-8154.72944624983;... % X-min |
                    3000;...              % X-max | Surveillance region bounding
                    -7548.44272179096;... % Y-min | box coordinates (m)
                    4355.32645897434]';   % Y-max |
        RadarCoords.lat = 50.36286670714617;
        RadarCoords.lon = -4.156833300366998;
end

% Load dataset
load(strcat(RadarName,'_old.mat'));
N = size(DataList,2); % Simulation length

%% Plot & Recording settings
% Plot settings
ShowPlots = 1;              % Set to 0 to prevent showing any plots
ShowUpdate = 1;             % Set to 0 to skip showing update plots
ShowTrackInfo = 0;
NumPersistFrames = 50;

% Recording settings
Record = 1;                 % Set to (0|1) to turn video recording (off|on)
FrameRate = 0.5;            % Number of frames per second
VideoQuality = 100;         % Set to desired quality percentage
VideoPathName = strcat(RadarName,'.avi'); % Set to the desired path and name of produced recording
%% Instantiation of necessary components

% Instantiate a Dynamic model (CV model with q = 1 m/s^2)
dyn = ConstantVelocityModelX_2D('VelocityErrVariance',0.3);
dyn.TimeVariant = 2;

% Instantiate an Observation model (Variance of 50m^2 on each coordinate)
obs = LinGaussObsModelX_2D('NumStateDims',4,'ObsErrVariance',50,'Mapping',[1 3]);

% Compile the State-Space model
ssm = StateSpaceModelX(dyn,obs);

V = (abs(V_bounds(2)-V_bounds(1))*abs(V_bounds(4)-V_bounds(3))); % Total area of surveillance region

% Assign PHD parameter values
config.NumParticles = 100000;
config.Model = ssm;
config.BirthIntFcn = @(Np)[abs(V_bounds(2)-V_bounds(1))*rand(1,Np)+V_bounds(1);... % Uniform position across the surveillance
                           5*rand(1,Np);...                                      % region.
                           abs(V_bounds(4)-V_bounds(3))*rand(1,Np)+V_bounds(3);... % Uniform speed between 0-9 m/s in both
                           5*rand(1,Np)];                                        % X and Y.
config.PriorDistFcn = @ (Np) deal(config.BirthIntFcn(Np), repmat(1/Np, Np, 1)');   % Uniform position and weights.
config.BirthScheme = {'Expansion', 0.1*config.NumParticles};
config.ProbOfDeath = 0.005;                                                        % Probability of death = 0.5%
config.ProbOfDetection = 0.6;                                                      % Probability of detection = 70%
config.ResamplingScheme = 'Multinomial';                                           % Use Multinomial Resampling

% Instantiate PHD filter
myphd = SMC_PHDFilterX(config);

% Instantiate PF filter
mypf = ParticleFilterX(ssm);

% Initiate PDAF parameters
Params_pdaf.Clusterer = NaiveClustererX();
Params_pdaf.Gater = EllipsoidalGaterX(2,'ProbOfGating',0.99)';
Params_pdaf.ProbOfDetect = 0.6;
Params_pdaf.Hypothesiser = LoopyBeliefPropagationX('ConvergeThreshold',10^(-3));
mypdaf = JointIntegratedProbabilisticDataAssocX(Params_pdaf);

% Initiate Track Initiator
config_ti.Filter = mypf;
config_ti.PHDFilter = myphd;
config_ti.ProbOfGating = 0.99;
config_ti.ProbOfConfirm = 0.9;
myti = PhdExistProbTrackInitiatorX(config_ti);

TrackList = [];

%% Create plot windows
if(ShowPlots)

    % Map plot
    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
    plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
    %axis(ax(1),V_bounds)
    axis(ax(1),[-4.195 -4.11 50.31 50.375])

    % PHD Intensity plot
    figure('units','normalized','outerposition',[.5 0 .5 1])
    ax(2) = gca;

    plots = [];
end

%% START OF SIMULATION
% ===================>
for k=1:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    % Extract DataList at time k
    tempDataList = DataList{k}(1:2,:);
    tempDataList( :, ~any(tempDataList,1) ) = [];

    % Process JPDAF
    mypdaf.MeasurementList = tempDataList;
    mypdaf.TrackList = TrackList;
    for j=1:numel(TrackList)
        mypdaf.TrackList{j}.Filter.predict();
    end
    mypdaf.associate();
    mypdaf.updateTracks();

    tic;

    % Append state to target trajectories
    for j=1:numel(TrackList)
        if(isempty(TrackList{j}.Trajectory))
            TrackList{j}.Trajectory = TrackList{j}.Filter.StateMean;
        else
            TrackList{j}.Trajectory = [TrackList{j}.Trajectory, TrackList{j}.Filter.StateMean];
        end
    end

    % Plot update step results
    if(ShowPlots && ShowUpdate)

        % Delete all plots (other than map)
        for i = 1:numel(plots)
            delete(plots(i))
        end
        plots = [];
        hold on;

        if(size(DataList{k},2)>0)
            % Convert measurements to LLA and plot them
            [lat,lon,~] = ned2geodetic(DataList{k}(2,:),...
                                       DataList{k}(1,:),...
                                       0,...
                                       RadarCoords.lat,...
                                       RadarCoords.lon,...
                                       0,...
                                       referenceEllipsoid('wgs84'));
    %         lat = DataList{k}(3,:);
    %         lon = DataList{k}(4,:);
            plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
            plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
                               '-s','MarkerSize',20,...
                               'MarkerEdgeColor','red',...
                               'MarkerFaceColor',[1 .6 .6]);
        end

        % Plot all existing tracks
        for j=1:numel(TrackList)

            % Convert track trajectory to LLA and plot it
            [lat,lon,~] = ned2geodetic(TrackList{j}.Trajectory(3,:),...
                                       TrackList{j}.Trajectory(1,:),...
                                       0,...
                                       RadarCoords.lat,...
                                       RadarCoords.lon,...
                                       0,...
                                       referenceEllipsoid('wgs84'));
            traj_length = size(lon,2);
            if(traj_length>NumPersistFrames)
                start = traj_length-NumPersistFrames;
            else
                start = 1;
            end
            plots(end+1) = plot(ax(1), lon(:,start:end),lat(:,start:end),'-.w','LineWidth',2);
            plots(end+1) = plot(ax(1), lon(:,end),lat(:,end),'ws','MarkerSize',15);

            % Convert track velocity to LLA and plot it
            [lat_vel,lon_vel,~] = ned2geodetic(TrackList{j}.Trajectory(4,end),...
                                       TrackList{j}.Trajectory(2,end),...
                                       0,...
                                       lat(:,end),...
                                       lon(:,end),...
                                       0,...
                                       referenceEllipsoid('wgs84'));
            lat_vel = lat_vel-lat(:,end);
            lon_vel = lon_vel-lon(:,end);

            plots(end+1) = quiver(ax(1), lon(:,end),lat(:,end),20*lon_vel,20*lat_vel,'r','LineWidth',1.5);
            speed_kmph = sqrt(TrackList{j}.Trajectory(4,end)^2+TrackList{j}.Trajectory(2,end)^2)*3.6;
            speed_knot = speed_kmph/1.852;

            if(ShowTrackInfo)
                plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)+0.00027,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','w');
                plots(end+1) = text(ax(1), lon(:,end)+0.001,lat(:,end)-0.00027,strcat("PoE:",num2str(TrackList{j}.ProbOfExist*100,3)," %"),'FontSize',8,'Color','w');
            end
            % TODO: Convert heading and state covariance to LLA and plot them
            % [lat,lon,h] = ned2geodetic(North,East,0,50.346069,-4.113670,0,referenceEllipsoid('wgs84'));
            % h2 = plot_gaussian_ellipsoid(TrackList{j}.Filter.StateMean([1 3]), TrackList{j}.Filter.StateCovar([1 3],[1 3]),1,20,ax(1));
            % plots(end+1) = h2;
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-5,int2str(TrackList{j}.TrackID));
            % plots(end+1) = text(ax(1),TrackList{j}.Filter.StateMean(1)+20,TrackList{j}.Filter.StateMean(3)-130,num2str(TrackList{j}.ProbOfExist,2));
        end

        % Add axis labels
        xlabel('Longitude')
        ylabel('Latitude')

        % Plot PHD Intensity
%         cla(ax(2), 'reset');
%         [bandwidth,density,X,Y]=kde2d(myphd.Particles([1,3],:)');
%         h = surf(ax(2),X,Y,density);
%         shading interp
%         colormap(ax(2), jet(3000))
%         hold on;
%         plot(ax(2), myphd.Particles(1,:), myphd.Particles(3,:), '.',...
%                     myphd.MeasurementList(1,:), myphd.MeasurementList(2,:), 'y*');
%         axis(ax(2), [V_bounds]);
%         str = sprintf('PHD intensity (Update)');
%         xlabel(ax(2),'X position (m)')
%         ylabel(ax(2),'Y position (m)')
%         zlabel(ax(2),'Intensity')
%         title(ax(2),str)
%         pause(0.01)

        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end

    % Perform Track initiation
    myti.MeasurementList = tempDataList; % New observations
    myti.PHDFilter.ClutterRate = size(myti.MeasurementList,2)/V;
    myti.TrackList = TrackList;
    myti.AssocWeightsMatrix = mypdaf.AssocWeightsMatrix;
    TrackList = myti.initiateTracks();
end

% Create video file and write to it
if(Record)
    %F = F(2:end);
    vidObj = VideoWriter(char(VideoPathName));
    vidObj.Quality = VideoQuality;
    vidObj.FrameRate = FrameRate;
    open(vidObj);
    writeVideo(vidObj, F);
    close(vidObj);
end

author: sglvladi

repo: TrackingX

save_path: github-repos/MATLAB/sglvladi-TrackingX


text:
function table = r8mat_latinize ( m, n, table )

%*****************************************************************************80
%
%% R8MAT_LATINIZE "Latinizes" an R8MAT.
%
%  Discussion:
%
%    It is assumed, though not necessary, that the input dataset
%    has points that lie in the unit hypercube.
%
%    In any case, the output dataset will have this property.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 March 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, integer M, the spatial dimension.
%
%    Input, integer N, the number of cells.
%
%    Input, real TABLE(M,N), the dataset to be "Latinized".
%
%    Output, real TABLE(M,N), the Latinized dataset.
%
  for i = 1 : m
    indx = r8vec_sort_heap_index_a ( n, table(i,1:n) );
    for j = 1 : n
      table(i,indx(j)) = ( 2 * j - 1 ) / ( 2 * n );
    end
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
function gp_save_loudspeakers(file,x0)
%GP_SAVE_LOUDSPEAKERS saves x0 as a text file in a Gnuplot compatible format
%
%   Usage: gp_save_loudspeakers(file,x0)
%
%   Input parameters:
%       file        - filename of the data file
%       x0          - secondary sources [nx7]
%
%   GP_SAVE_LOUDSPEAKERS(file,x0) saves x0(:,1:2) as positions of the
%   loudspeakers, an orientation value calculated from x0(:,4:6), and the
%   activity x0(:,7) of the loudspeakers in a text file useable by gnuplot.
%
%   See also: gp_save, gp_save_matrix

%*****************************************************************************
% The MIT License (MIT)                                                      *
%                                                                            *
% Copyright (c) 2010-2019 SFS Toolbox Developers                             *
%                                                                            *
% Permission is hereby granted,  free of charge,  to any person  obtaining a *
% copy of this software and associated documentation files (the "Software"), *
% to deal in the Software without  restriction, including without limitation *
% the rights  to use, copy, modify, merge,  publish, distribute, sublicense, *
% and/or  sell copies of  the Software,  and to permit  persons to whom  the *
% Software is furnished to do so, subject to the following conditions:       *
%                                                                            *
% The above copyright notice and this permission notice shall be included in *
% all copies or substantial portions of the Software.                        *
%                                                                            *
% 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.                                                  *
%                                                                            *
% The SFS Toolbox  allows to simulate and  investigate sound field synthesis *
% methods like wave field synthesis or higher order ambisonics.              *
%                                                                            *
% https://sfs.readthedocs.io                            sfstoolbox@gmail.com *
%*****************************************************************************

%% ===== Checking of input  parameters ==================================
nargmin = 2;
nargmax = 2;
narginchk(nargmin,nargmax);
isargchar(file);
isargsecondarysource(x0);

%% ===== Main ============================================================
% Write header to the file
fid = fopen(file,'w');
fprintf(fid,'# Loudspeaker file generated by gp_save_loudspeakers.m\n');
fprintf(fid,'# x0 y0 phi ls_activity\n');
fclose(fid);

% Calculate phi
loudspeaker(:,1:2) = x0(:,1:2);
[loudspeaker(:,3),~,~] = cart2sph(x0(:,4),x0(:,5),x0(:,6));
loudspeaker(:,4) = x0(:,7);

% Append the data to the file using tabulator as a delimiter between the data
if isoctave
    dlmwrite(file,loudspeaker,'\t','-append');
else
    dlmwrite(file,loudspeaker,'delimiter','\t','-append');
end

author: sfstoolbox

repo: sfs-matlab


text:
function [FF,II] = remesh_at_points(V,F,P)
  % REMESH_AT_POINTS Given a surface mesh (V,F) and a list of points on/near the
  % surface P, subdivide triangles of (V,F) so that P are now contained in the
  % vertex set.
  %
  % Inputs:
  %   V  #V by 3 list of vertex positions
  %   F  #F by 3 list of face indices into V
  %   P  #P by 3 list of point positions
  % Outputs:
  %   FF  #FF by 3 list of face indices into [V;P]
  %   II  #FF list of indices into F revealing birth face
  %
  % Example:
  %   P = random_points_on_mesh(V,F,size(F,1));
  %   [FF,II] = remesh_at_points(V,F,P);
  %   tsurf(FF,[V;P],'CData',II)
  %

  VV = V;
  FF = F;
  II = (1:size(F,1))';
  JJ = zeros(size(V,1),1);
  L = (1:size(P,1))';
  PP = P;
  while true
    %clf;hold on;tsurf(FF,VV,'CData',II);scatter3(PP(:,1),PP(:,2),PP(:,3));hold off;pause
    % Find closest points/faces on mesh (This is a bit redundant. Really we
    % should only look at faces that _were_ affected last round, since those
    % must have had multiple points.)
    [~,IC,C] = point_mesh_squared_distance(PP,VV,FF);
    % Only keep one closest point per face
    [I,J] = unique(IC);
    % append keepers
    n = size(VV,1);
    JJ = [JJ;L(J)];
    VV = [VV;C(J,:)];
    nCI = n+(1:numel(J))';
    unaffected = setdiff(1:size(FF,1),I);
    II = [II(unaffected);repmat(II(I),3,1)];
    FF = [ ...
      FF(unaffected,:);
      FF(I,[1 2]) nCI; ...
      FF(I,[2 3]) nCI; ...
      FF(I,[3 1]) nCI];
    % Leftovers
    L = L(setdiff(1:end,J));
    PP = P(L,:);
    if isempty(PP)
      break;
    end
  end

  %% Re-order so that output is simply VV = [V;C]
  %C = zeros(size(P));
  %C(JJ(size(V,1)+1:end),:) = VV(size(V,1)+1:end,:);
  I = JJ+size(V,1);
  I(1:size(V,1)) = 1:size(V,1);
  FF = I(FF);
  %VV = [V;C];
end

author: alecjacobson

repo: gptoolbox


text:
function variance = power_variance ( a, b )

%*****************************************************************************80
%
%% POWER_VARIANCE returns the variance of the Power PDF.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 September 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real A, B, the parameters of the PDF.
%    0.0 < A, 0.0 < B.
%
%    Output, real VARIANCE, the variance of the PDF.
%
  variance = b * b * a / ( ( a + 1.0 )^2 * ( a + 2.0 ) );

  return
end


text:
function stamps(start_step,end_step,patches_flag,est_gamma_parm,patch_list_file,stamps_PART_limitation)
%STAMPS Stanford Method for Persistent Scatterers
%   STAMPS(START_STEP,END_STEP,PATCHES_FLAG,EST_GAMMA_FLAG) Default is to run all steps.
%   A subset of steps may be selected with START_STEP and/or END_STEP
%   STEP 1 = Initial load of data
%   STEP 2 = Estimate gamma
%   STEP 3 = Select PS pixels
%   STEP 4 = Weed out adjacent pixels
%   STEP 5 = Correct wrapped phase for spatially-uncorrelated look angle error and merge patches
%   STEP 6 = Unwrap phase
%   STEP 7 = Calculate spatially correlated look angle (DEM) error
%   STEP 8 = Filter spatially correlated noise
%   STEP 0 = Continue from the last known stage till the end-stage selected
%
%   PATCHES_FLAG Default 'y'. Set to 'n' to process all data as one patch
%
%   EST_GAMMA_PARM is an optional parameter passed to PS_EST_GAMMA_QUICK
%
%   PATCH_LIST_FILE is an optional argument specifying the file list of
%   patches to be processed. Note that from step 5 and above one should use
%   all patches to merge results.
%
%   If current directory is a single patch, stamps only operates in the
%   current directory, but if current directory contains many patches,
%   stamps operates on them all.
%
%   Andy Hooper, June 2006

%   =================================================================
%   07/2006 AH: END_STEP added
%   09/2006 AH: ps_load removed (obsolete)
%   09/2006 AH: small baselines added
%   11/2006 AH: patches added
%   01/2007 AH: calculate spatially correlated look angle error added
%   03/2009 AH: simultaneously estimate velocity when SCLA estimated
%   03/2009 AH: smooth SCLA for unwrapping iteration
%   03/2010 AH: move ps_cal_ifg_std to after merge step
%   12/2012 AH: add gamma option
%   12/2012 DB: add patch_list_file argument as option
%   09/2013 DB: update the stamps version number
%   09/2015 DB: Check if patches do have PS before proceeding with
%               processing.
%   09/2015 DB: Fix when running stamps in a patch folder mode when no PS are left
%   09/2015 AH: allow for non-differentiation of caps by dir
%   01/2016 DB: include stamps_save in step 1-4.
%   08/2016 AH: Fix bug of scn_kriging_flag not being set
%   06/2017 DB: Catching when no PS are left from step 1, allow for re-run
%               when parameters have changed.
%   06/2017 DB: Option to continue from last know processing step
%   08/2107 AH: Removed catch as proceeds also when error in Step 1
%   =================================================================

nfill=40;
fillstr=[repmat('#',1,nfill),'\n'];
skipstr='\n';
msgstr=fillstr;

fprintf(skipstr);
logit(fillstr);
msgstr(round(nfill)/2-12:round(nfill/2)+13)=' StaMPS/MTI Version 4.0b6 ';
logit(msgstr);
msgstr(round(nfill)/2-12:round(nfill/2)+13)='  Beta version, Jun 2018  ';
logit(msgstr);
logit(fillstr);
fprintf(skipstr);

quick_est_gamma_flag=getparm('quick_est_gamma_flag');
reest_gamma_flag=getparm('select_reest_gamma_flag');
unwrap_method=getparm('unwrap_method');
unwrap_prefilter_flag=getparm('unwrap_prefilter_flag');
small_baseline_flag=getparm('small_baseline_flag');
insar_processor=getparm('insar_processor');
scn_kriging_flag=getparm('scn_kriging_flag');

if nargin<1 || isempty(start_step)==1
    start_step=1;
end

if nargin<2 || isempty(end_step)==1
    end_step=8;
end

if nargin<3 || isempty(patches_flag)==1
    if start_step<6
        patches_flag='y';
    else
        patches_flag='n';
    end
end

if nargin<4 || isempty(est_gamma_parm)==1
    est_gamma_parm=0;
end

if nargin<5 || isempty(patch_list_file)     % [DB] allow for own specified patch list file
    patch_list_file = 'patch.list';
    new_patch_file = 0;
else
    % use own file
    new_patch_file = 1;
end

% In support of the multi-core option limit processing to steps 1-5a,
% 5b-upwards, or the old method where all is processed together.
if nargin<6 || isempty(stamps_PART_limitation)
    stamps_PART_limitation=0;
end
stamps_PART1_flag='y';
stamps_PART2_flag='y';
if stamps_PART_limitation==1
    stamps_PART2_flag='n';
end
if stamps_PART_limitation==2
    stamps_PART1_flag='n';
end

if strcmpi(patches_flag,'y')
    if exist(patch_list_file,'file')
        fid=fopen(patch_list_file);
        i=0;
        while 1
            nextline=fgetl(fid);
            if ischar(nextline)
                i=i+1;
                patchdir(i).name=nextline;
            else
                break
            end
        end
	fclose(fid)
    else
        patchdir=dir('PATCH_*');
        patchdir = patchdir(find(~cellfun(@(x) strcmpi(x,'patch_noover.in'),{patchdir(:).name})));
    end
    if isempty(patchdir)
        patches_flag='n';
    else
        ps_parms_default
        patches_flag='y';
    end
end

if ~strcmpi(patches_flag,'y')
    patchdir(1).name='.';
    logit('Will process current directory only')
else
    logit('Will process patch subdirectories')
end

currdir=pwd;

nfill=40;
fillstr=[repmat('#',1,nfill),'\n'];
msgstr=fillstr;

% limit the processing to step 1-5a
start_step_or = start_step;
if strcmpi(stamps_PART1_flag,'y')
  for i=1:length(patchdir)
    if ~isempty(patchdir(i).name)
      cd(patchdir(i).name)
      patchsplit=strsplit(pwd,'/');
      %fprintf(skipstr);
      %logit(sprintf('Processing %s',patchsplit{end}))

      % store if patch dir is empty
      if exist('no_ps_info.mat','file')~=2
         stamps_step_no_ps = zeros([5 1 ]);       % keep for the first 5 steps only
         save('no_ps_info.mat','stamps_step_no_ps')
      end

        % if start_step is 0, then start from the latest stage it was
        if start_step_or==0
            % check the processing stage of stamps for all the patches
            % step 4 find a ps_weed file
            % step 3 find a ps_select file
            % step 2 find a pm file
            % step 1 find a ps file
            % or no PS in case stamps_step_no_ps is found with a 1 in there

            if exist('weed1.mat','file')==2
                   start_step=5;
                   setpsver(2);
            elseif exist('select1.mat','file')==2
                   start_step=4;
            elseif exist('pm1.mat','file')==2
                   start_step=3;
            elseif exist('ps1.mat','file')==2
                   start_step=2;
            else
                start_step=1;
            end

            if start_step>end_step
                fprintf(['\n' patchsplit{end} ': already up to end stage ' num2str(end_step) ' \n'])
            else
                fprintf(['\n' patchsplit{end} ': complete up to stage ' num2str(start_step-1) ' \n'])
            end
        end

      if start_step==1
        msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 1 ';
        fprintf(skipstr);
        logit(fillstr);
        logit(msgstr);
        logit(fillstr)
        logit(['Directory is ',patchsplit{end}])
        fprintf(skipstr);
        if strcmpi(small_baseline_flag,'y')
%             try
                if strcmpi(insar_processor,'gamma') | strcmpi(insar_processor,'snap')
                    sb_load_initial_gamma;
                elseif strcmpi(insar_processor,'gsar')
                    sb_load_initial_gsar;
                elseif  strcmpi(insar_processor,'isce')
                    if exist('data_inc','var')==0
                        % already in patch dir, file contained in the InSAR dir
                        inc_angle = ['..' filesep 'inc_angle.raw'];
                        if exist(inc_angle,'file')~=2
                             inc_angle = ['..' filesep inc_angle];
                        end
                        if exist(inc_angle,'file')==2
                            fprintf('Found inc angle file, will load the data \n')
                            data_inc = (load_isce(inc_angle));
                        else
                            data_inc=[];
                        end
                    end
                    sb_load_initial_isce(data_inc)
                else
                    sb_load_initial;
                end
                load('no_ps_info.mat');
                % reset as we are currently re-processing
                stamps_step_no_ps(1:end)=0;

%             catch
%
%                load('no_ps_info.mat');
%                % reset as we are currently re-processing
%                stamps_step_no_ps(1:end)=0;
%                fprintf('***No PS points left. Updating the stamps log for this****\n')
%                % update the flag indicating no PS left in step 1
%                stamps_step_no_ps(1)=1;
%                psver =1;
%                save('psver.mat','psver')
%
%             end
            save('no_ps_info.mat','stamps_step_no_ps')

        else
%             try
                if strcmpi(insar_processor,'gamma') | strcmpi(insar_processor,'snap')
                    ps_load_initial_gamma;
                elseif strcmpi(insar_processor,'gsar')
                    ps_load_initial_gsar;
                elseif  strcmpi(insar_processor,'isce')
                     if exist('data_inc','var')==0
                        % already in patch dir, file contained in the InSAR dir
                        inc_angle = ['..' filesep 'inc_angle.raw'];
                        if exist(inc_angle,'file')~=2
                             inc_angle = ['..' filesep inc_angle];
                        end
                        if exist(inc_angle,'file')==2
                            fprintf('Found inc angle file, will load the data \n')
                            data_inc = (load_isce(inc_angle));
                        else
                            data_inc=[];
                        end
                    end
                    ps_load_initial_isce(data_inc)

                else
                    ps_load_initial;
                end
                load('no_ps_info.mat');
                % reset as we are currently re-processing
                stamps_step_no_ps(1:end)=0;
%             catch
%                 load('no_ps_info.mat');
%                 % reset as we are currently re-processing
%                 stamps_step_no_ps(1:end)=0;
%                 fprintf('***No PS points left. Updating the stamps log for this****\n')
%                 % update the flag indicating no PS left in step 1
%                 stamps_step_no_ps(1)=1;
%                 save('no_ps_info.mat','stamps_step_no_ps')
%                 psver =1;
%                 save('psver.mat','psver')
%
%             end
            save('no_ps_info.mat','stamps_step_no_ps')
        end
        elseif start_step <=4
            setpsver(1)
      end

        if start_step<=2 & end_step >=2
            msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 2 ';
            fprintf(skipstr);
            logit(fillstr);
            logit(msgstr);
            logit(fillstr)
            logit(['Directory is ',patchsplit{end}])
            fprintf(skipstr);

            % check if step 1 had more than 0 PS points
            load('no_ps_info.mat');
            % reset as we are currently re-processing
            stamps_step_no_ps(2:end)=0;

            % run step 2 when there are PS left in step 1
            if stamps_step_no_ps(1)==0
                if strcmpi(quick_est_gamma_flag,'y')
                    ps_est_gamma_quick(est_gamma_parm);
                else
                    ps_est_gamma(est_gamma_parm);
                end
            else
                stamps_step_no_ps(2)=1;
                fprintf('No PS left in step 1, so will skip step 2 \n')
            end
            save('no_ps_info.mat','stamps_step_no_ps')
        end

        if start_step<=3 & end_step >=3
            msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 3 ';
            fprintf(skipstr);
            logit(fillstr);
            logit(msgstr);
            logit(fillstr)
            logit(['Directory is ',patchsplit{end}])
            fprintf(skipstr);

            % check if step 2 had more than 0 PS points
            load('no_ps_info.mat');
            % reset as we are currently re-processing
            stamps_step_no_ps(3:end)=0;

            % run step 3 when there are PS left in step 2
            if stamps_step_no_ps(2)==0
                if strcmpi(quick_est_gamma_flag,'y') & strcmpi(reest_gamma_flag,'y')
                    ps_select;
                else
                    ps_select(1);
                end
            else
                fprintf('No PS left in step 2, so will skip step 3 \n')
                stamps_step_no_ps(3)=1;
            end
            save('no_ps_info.mat','stamps_step_no_ps')
        end

        if start_step<=4 & end_step >=4
            msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 4 ';
            fprintf(skipstr);
            logit(fillstr);
            logit(msgstr);
            logit(fillstr)
            logit(['Directory is ',patchsplit{end}])
            fprintf(skipstr);

            % check if step 3 had more than 0 PS points
            load('no_ps_info.mat');
            % reset as we are currently re-processing
            stamps_step_no_ps(4:end) =0;       % keep for the first 5 steps only

            % run step 4 when there are PS left in step 3
            if stamps_step_no_ps(3)==0
                if strcmpi(small_baseline_flag,'y')
                    ps_weed(0,1);
                else
                    ps_weed;
                end
            else
                fprintf('No PS left in step 3, so will skip step 4 \n')
                stamps_step_no_ps(4)=1;

            end
            save('no_ps_info.mat','stamps_step_no_ps')
        end

        if start_step<=5 & end_step >=5
            msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 5 ';
            fprintf(skipstr);
            logit(fillstr);
            logit(msgstr);
            logit(fillstr)
            logit(['Directory is ',patchsplit{end}])
            fprintf(skipstr);

            % check if step 4 had more than 0 PS points
            load('no_ps_info.mat');
            % reset as we are currently re-processing
            stamps_step_no_ps(5:end) = 0;       % keep for the first 5 steps only

            % run step 5 when there are PS left in step 3
            if stamps_step_no_ps(4)==0
                ps_correct_phase;
            else
                fprintf('No PS left in step 4, so will skip step 5 \n')
                stamps_step_no_ps(5)=1;
            end
            save('no_ps_info.mat','stamps_step_no_ps')
        end

        cd(currdir)
      end
    end
end

patchsplit=strsplit(pwd,'/');

% check if one can process second part of step 5b and above
if strcmpi(stamps_PART2_flag,'y')
    %%% Loop throught the patches and update the patch.list and keep only those
    %%% that have PS left.
    if patches_flag=='y'
        % go in reverse order such patches can be dropped when needed

        fid = fopen('patch.list_new','w');
        for i=1:length(patchdir)
            % check the file with the PS information
            filename_PS_check = [patchdir(i).name filesep 'no_ps_info.mat'];

            % assume by default to keep patch for backward compatibility
            keep_patch = 1;
            if exist(filename_PS_check,'file')==2
                load(filename_PS_check)
                if sum(stamps_step_no_ps)>=1
                   keep_patch=0;
                end
            end

            % update the patch list.
            if keep_patch==1
                fprintf(fid,[patchdir(i).name '\n']);
            end
            if i==length(patchdir)
               fclose(fid) ;
            end
        end

        % update the files such in futhre the new patch list will be used.
        movefile('patch.list','patch.list_old');
        movefile('patch.list_new','patch.list');
    end

    if start_step<=5 & end_step >=5
        abord_flag=0;
        if patches_flag=='y'
            fprintf(skipstr);
            logit(['Directory is ',patchsplit{end}])
            fprintf(skipstr);
            ps_merge_patches
        else
            % this is processing of an individual patch
            % see if there are any PS left
            if exist('no_ps_info.mat','file')==2
                load('no_ps_info.mat')
                if sum(stamps_step_no_ps)>=1
                   abord_flag=1;
                end
            end
        end

        % see if step 5 can be ran
        if abord_flag==0
            ps_calc_ifg_std;
        else
            fprintf('No PS left in step 4, so will skip step 5 \n')
        end
    end

    if start_step<=6 & end_step >=6
        msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 6 ';
        fprintf(skipstr);
        logit(fillstr);
        logit(msgstr);
        logit(fillstr)
        logit(['Directory is ',patchsplit{end}])
        fprintf(skipstr);

        ps_unwrap
        if strcmpi(small_baseline_flag,'y')
            sb_invert_uw
        end
    end

    if start_step<=7 & end_step >=7
        msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 7 ';
        fprintf(skipstr);
        logit(fillstr);
        logit(msgstr);
        logit(fillstr)
        logit(['Directory is ',patchsplit{end}])
        fprintf(skipstr);

        if strcmpi(small_baseline_flag,'y')
            ps_calc_scla(1,1)   % small baselines
            ps_smooth_scla(1)
            ps_calc_scla(0,1) % single master
        else
            ps_calc_scla(0,1)
            ps_smooth_scla
        end
    end

    if start_step<=8 & end_step >=8
        msgstr(round(nfill)/2-3:round(nfill/2)+4)=' Step 8 ';
        fprintf(skipstr);
        logit(fillstr);
        logit(msgstr);
        logit(fillstr)
        logit(['Directory is ',patchsplit{end}])
        fprintf(skipstr);

        if strcmpi(scn_kriging_flag,'y')
            ps_scn_filt_krig
        else
            ps_scn_filt
        end
    end
end

logit(1);

author: dbekaert

repo: StaMPS

save_path: github-repos/MATLAB/dbekaert-StaMPS


text:
function tfplot(imp);

% tfplot(imp);
%
% Plot the transfer function of the filter whose impulse response is imp.
%
% Multiple transfer functions can be plotted simultaneously where each
% impulse response represents one row of imp.

% Copyright 2012 Evrytania LLC (http://www.evrytania.com)
%
% Written by James Peroulas <james@evrytania.com>
%
% 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/>.

error(nargchk(1,1,nargin));
%error(chk_param(imp,'imp','vector'));

n_pts_min=4096;
if (isvector(imp))
  n_pts=max(n_pts_min,length(imp));
else
  n_pts=max(n_pts_min,size(imp,2));
end

f=linspace(0,2,n_pts+1);
f=f(1:end-1);
f(f>=1)=f(f>=1)-2;
f=fftshift(f);
m=transpose(fftshift(fft(transpose(imp),n_pts)));
plot(f,db20(abs(m)));

xlabel('Frequency/(fs/2)');
ylabel('Power gain (dB)');
title('Transfer function');
zgo;

author: JiaoXianjun

repo: rtl-sdr-LTE


text:
function normalizeUnits(gdata,fontUnitsExcludeTags)
% normalizeUnits(gdata)
% Change 'Units' properties of all handles contained in gdata to
% 'Normalized' (the 'FontUnits' property is unaffected). gdata is a
% structure with fields whose values are handles to GUI components (eg as
% returned by guidata).

if nargin < 2
    fontUnitsExcludeTags = cell(0,1);
end

comps = fieldnames(gdata);
for c = 1:numel(comps)
    tag = comps{c};
    h = gdata.(tag);
    if isprop(h,'Units')
        set(h,'Units','Normalized');
    end
    if isprop(h,'FontUnits') && ~any(strcmp(tag,fontUnitsExcludeTags))
        set(h,'FontUnits','Normalized');
    end
end

end

author: kristinbranson

repo: JAABA


text:
function [vX] = spm_vec(X,varargin)
% Vectorise a numeric, cell or structure array - a compiled routine
% FORMAT [vX] = spm_vec(X)
% X  - numeric, cell or stucture array[s]
% vX - vec(X)
%
% See spm_unvec
%__________________________________________________________________________
%
% e.g.:
% spm_vec({eye(2) 3}) = [1 0 0 1 3]'
%__________________________________________________________________________
% Copyright (C) 2005-2013 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_vec.m 6110 2014-07-21 09:36:13Z karl $

%error('spm_vec.c not compiled - see Makefile')

% initialise X and vX
%--------------------------------------------------------------------------
if nargin > 1
    X = [{X},varargin];
end

% vectorise numerical arrays
%--------------------------------------------------------------------------
if isnumeric(X)
    vX = X(:);

% vectorise logical arrays
%--------------------------------------------------------------------------
elseif islogical(X)
    vX = X(:);

% vectorise structure into cell arrays
%--------------------------------------------------------------------------
elseif isstruct(X)
    vX = [];
    f   = fieldnames(X);
    X    = X(:);
    for i = 1:numel(f)
        vX = cat(1,vX,spm_vec({X.(f{i})}));
    end

% vectorise cells into numerical arrays
%--------------------------------------------------------------------------
elseif iscell(X)
    vX   = [];
    for i = 1:numel(X)
        vX = cat(1,vX,spm_vec(X{i}));
    end
else
    vX = [];
end

author: spm


text:
scale1=[16,32,64,128,256,512,1024];
X1=cumsum(multifractal-mean(multifractal));
for ns=1:length(scale1),
    segments1(ns)=floor(length(X1)/scale1(ns));
    for v=1:segments1(ns),
        Index1=((((v-1)*scale1(ns))+1):(v*scale1(ns)));
        RMS_scale2{ns}(Index1)=RMS{ns}(v).*ones(size(Index1));
    end
end
clear X1 Index1 scale1 segments1
X1=1:length(RMS_scale2{7});
X2=1:8000;
X3=1:length(RMS_scale2{6});
X4=1:length(RMS_scale2{5});
Y1=RMS_scale2{7};
Y2=F(7).*ones(8000,1);
Y3=RMS_scale2{6};
Y4=F(6).*ones(8000,1);
Y5=RMS_scale2{5};
Y6=F(5).*ones(8000,1);
Y7=RMS_scale2{4};
Y8=F(4).*ones(8000,1);
YMatrix1=[RMS_scale2{3};F(3).*ones(1,8000)];
YMatrix2=[RMS_scale2{2};F(2).*ones(1,8000)];
YMatrix3=[RMS_scale2{1};F(1).*ones(1,8000)];
clear RMS_scale2

% Create figure
figure1 = figure('PaperSize',[20.98 29.68],'Color',[1 1 1]);

% Create axes
axes1 = axes('Parent',figure1,'YTick',[0 20],'XTick',zeros(1,0),...
    'Position',[0.09852 0.8515 0.775 0.1274],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes1,[0 45]);
hold(axes1,'all');

% Create plot
plot(X1,Y1,'Parent',axes1,'LineWidth',2,'Color',[0 0 1],...
    'DisplayName','Local fluctuation: RMS in Matlab code 5');

% Create plot
plot(X2,Y2,'Parent',axes1,'LineWidth',2,'Color',[1 0 0],...
    'DisplayName','RMS of local fluctuation: F in Matlab code 5');

% Create axes
axes2 = axes('Parent',figure1,'YTick',[0 20],'XTick',zeros(1,0),...
    'Position',[0.09852 0.7242 0.775 0.1272],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes2,[0 45]);
hold(axes2,'all');

% Create plot
plot(X3,Y3,'Parent',axes2,'LineWidth',2,'Color',[0 0 1]);

% Create plot
plot(X2,Y4,'Parent',axes2,'LineWidth',2,'Color',[1 0 0]);

% Create axes
axes3 = axes('Parent',figure1,'YTick',[0 20],'XTick',zeros(1,0),...
    'Position',[0.09852 0.5968 0.775 0.1271],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes3,[0 44]);
hold(axes3,'all');

% Create plot
plot(X4,Y5,'Parent',axes3,'LineWidth',2,'Color',[0 0 1]);

% Create plot
plot(X2,Y6,'Parent',axes3,'LineWidth',2,'Color',[1 0 0]);

% Create axes
axes4 = axes('Parent',figure1,'YTick',[0 20],'XTick',zeros(1,0),...
    'Position',[0.09852 0.4694 0.775 0.127],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes4,[0 40]);
hold(axes4,'all');

% Create plot
plot(X4,Y7,'Parent',axes4,'LineWidth',2,'Color',[0 0 1]);

% Create plot
plot(X2,Y8,'Parent',axes4,'LineWidth',2,'Color',[1 0 0]);

% Create ylabel
ylabel('RMS in Matlab code 5 (measurement units)','FontSize',14);

% Create axes
axes5 = axes('Parent',figure1,'YTick',[0 10],'XTick',zeros(1,0),...
    'Position',[0.09852 0.3435 0.775 0.1269],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes5,[-0.5 30]);
hold(axes5,'all');

% Create multiple lines using matrix input to plot
pplot1 = plot(X2,YMatrix1,'Parent',axes5,'LineWidth',2);
set(pplot1(1),'Color',[0 0 1]);
set(pplot1(2),'Color',[1 0 0]);

% Create axes
axes6 = axes('Parent',figure1,'YTick',[0 10],'XTick',zeros(1,0),...
    'Position',[0.09852 0.2162 0.775 0.1267],...
    'LineWidth',2,...
    'FontSize',14);
xlim(axes6,[0 8000]);
ylim(axes6,[-0.5 20]);
hold(axes6,'all');

% Create multiple lines using matrix input to plot
pplot2 = plot(X2,YMatrix2,'Parent',axes6,'LineWidth',2);
set(pplot2(1),'Color',[0 0 1]);
set(pplot2(2),'Color',[1 0 0]);

% Create axes
axes7 = axes('Parent',figure1,'YTick',[0 5],...
    'Position',[0.09852 0.09048 0.775 0.1266],...
    'LineWidth',2,...
    'FontSize',14);
ylim(axes7,[0 15]);
hold(axes7,'all');

% Create multiple lines using matrix input to plot
pplot3 = plot(X2,YMatrix3,'Parent',axes7,'LineWidth',2);
set(pplot3(1),'Color',[0 0 1]);
set(pplot3(2),'Color',[1 0 0]);

% Create xlabel
xlabel('time (sample number)','FontSize',14);

% Create legend
legend1 = legend(axes1,'show');
set(legend1,'Position',[0.5252 0.8089 0.3681 0.07268],'LineWidth',2);

% Create textbox
annotation(figure1,'textbox',[0.8714 0.3324 0.0981 0.05573],...
    'String',{'Scale = 64'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8711 0.08541 0.0981 0.05573],...
    'String',{'Scale = 16'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8716 0.2025 0.0981 0.05573],...
    'String',{'Scale = 32'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8711 0.4607 0.102 0.05573],...
    'String',{'Scale = 128'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8716 0.6065 0.1049 0.05573],...
    'String',{'Scale = 256'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8721 0.7507 0.1036 0.05573],...
    'String',{'Scale = 512'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');

% Create textbox
annotation(figure1,'textbox',[0.8719 0.8997 0.1125 0.05573],...
    'String',{'Scale = 1024'},...
    'FontSize',14,...
    'FitBoxToText','off',...
    'LineStyle','none');
clear ans pplot1 pplot2 pplot3 axes1 axes2 axes3 axes4 axes5 axes6 axes7 legend1 figure1 X1 X2 X3 X4 Y1 Y2 Y3 Y4 Y5 Y6 Y7 Y8 YMatrix1 YMatrix2 YMatrix3
v=1;ns=1;


text:
function [specGlobal] = doa_srp(x,method, Param)
%%
if(~any(strcmp(method, {'SRP-PHAT' 'SRP-NON'})))
    error('ERROR[doa_srp]: method参数错误');
end
%% STFT
X = ssl_stft(x.',Param.window, Param.noverlap, Param.nfft, Param.fs);
X = X(2:end,:,:);
%%
if strcmp(method,'SRP-PHAT')
    specGlobal = ssl_srpPhat(X,Param);
else
    specGlobal = ssl_srp_nonlin(X,Param);
end

end

function X=ssl_stft(x,window,noverlap,nfft,fs)

% Inputs:x: nchan x nsampl  window = blackman(wlen);
% Output:X: nbin x nfram x nchan matrix

[nchan,~]=size(x);
[Xtemp,F,T,~] = spectrogram(x(1,:),window,noverlap,nfft,fs); % S nbin x nframe
nbin = length(F);
nframe = length(T);
X = zeros(nbin,nframe,nchan);
X(:,:,1) = Xtemp;
for ichan = 2:nchan
    X(:,:,ichan) = spectrogram(x(ichan,:),window,noverlap,nfft,fs);
end

end

function [specGlobal] = ssl_srpPhat(X,Param)
[~,nFrames,~] = size(X);
specInst = zeros(Param.nGrid, nFrames);

for i = 1:Param.nPairs
    spec = srpPhat_spec(X(Param.freqBins,:,Param.pairId(i,:)), Param.f(Param.freqBins), Param.tauGrid{i}); % NV % [freq x fram x local angle for each pair]
    specSampledgrid = (shiftdim(sum(spec,1)))';
    specCurrentPair = interp1q(Param.alphaSampled{i}', specSampledgrid, Param.alpha(i,:)');
    specInst(:,:) = specInst(:,:) + specCurrentPair;
end

switch Param.pooling
    case 'max'
        specGlobal = shiftdim(max(specInst,[],2));
    case 'sum'
        specGlobal = shiftdim(sum(specInst,2));
end
end

function [specGlobal] = ssl_srp_nonlin(X,Param)

alpha_meth = (10*Param.c)./(Param.d*Param.fs);
[~,nFrames,~] = size(X);
specInst = zeros(Param.nGrid, nFrames);

for i = 1:Param.nPairs
    spec = srpNonlin_spec(X(Param.freqBins,:,Param.pairId(i,:)), Param.f(Param.freqBins), alpha_meth(i), Param.tauGrid{i});
    specSampledgrid = (shiftdim(sum(spec,1)))';
    specCurrentPair = interp1q(Param.alphaSampled{i}', specSampledgrid, Param.alpha(i,:)');
    specInst = specInst + specCurrentPair;
end

switch Param.pooling
    case 'max'
        specGlobal = shiftdim(max(specInst,[],2));
    case 'sum'
        specGlobal = shiftdim(sum(specInst,2));
end
end

author: WenzheLiu-Speech


text:
%% Discrete-time reentry dynamics demonstration with non-linear filters
%
%  Description:
%    In this example various different non-linear filters and smoothers are
%    applied to reentry tracking problem . The filters used in this
%    demonstration are:
%      * Extended Kalman filter
%      * Unscented Kalman filter
%      * Gauss-Hermite Kalman filter (degree 3)
%      * Cubature Kalman filter
%    Additionally, the corresponding smoother results are also presented.
%
%  References:
%    Refer to the Toolbox documentation for details on the model.
%
%  See also:
%    ukf_predict1, ukf_update1, urts_smooth1,
%    ekf_predict1, ekf_update1, erts_smooth1, ghkf_predict, ghkf_update,
%    ghrts_smooth, ckf_predict, ckf_update, crts_smooth
%
%  Author:
%    Copyright (C) 2006 Simo Särkkä
%                  2007 Jouni Hartikainen
%                  2010 Arno Solin
%
%  Licence:
%    This software is distributed under the GNU General Public
%    Licence (version 2 or later); please refer to the file
%    Licence.txt, included with the software, for details.

%% Set up model parameters

  reentry_param;
  make_reentry_data;

  silent = 0;

  % Handles to dynamic and measurement model functions,
  % and to their derivatives.
  func_f  = @reentry_f;
  func_if = @reentry_if;
  func_df = @reentry_df_dx;
  func_h  = @reentry_h;
  func_dh = @reentry_dh_dx;

  % Initial values and space for EKF
  m = m0;
  P = P0;
  Q = L*Qc*L'*dt;

  MM_EKF = zeros(size(m,1),size(Y,2));
  PP_EKF = zeros(size(m,1),size(m,1),size(Y,2));
  VV_EKF = zeros(size(m,1),size(Y,2));
  EE_EKF = zeros(size(m,1),size(Y,2));

  % Check derivatives (should be OK)
  %der_check(func_a, func_da, 1, m0, {dt,b0,H0,Gm0,R0});
  %der_check(func_h, func_dh, 1, m0, {xr,yr});

  clf; clc;
  disp(['This is a demonstration for tracking a reentry vehicle ',...
        'using 1st order EKF and augmented UKF.'])
  disp(' ');

%% Extended Kalman filter

  fprintf('Running EKF...');
  % Filtering with EKF
  for k=1:size(Y,2)
    [m,P] = ekf_predict1(m,P,func_df,Q,func_f,[],{dt,b0,H0,Gm0,R0});
    [m,P] = ekf_update1(m,P,Y(:,k),func_dh,diag([vr va]),func_h,[],{xr,yr});
    MM_EKF(:,k) = m;
    PP_EKF(:,:,k) = P;
    VV_EKF(:,k) = diag(P);
    EE_EKF(:,k) = (X(:,k) - m).^2;
  end

  %
  % Calculate RMSE
  %
  ekf_rmse = sqrt(mean(sum((X(1:2,:)-MM_EKF(1:2,:)).^2)));
  ME_EKF = squeeze(PP_EKF(1,1,:)+PP_EKF(2,2,:));
  fprintf('Done!\n')

  fprintf('Running smoothers...');
  %
  % Smoother 1
  %
  [SM_ERTS,SP_ERTS] = erts_smooth1(MM_EKF,PP_EKF,func_df,Qc*dt,func_f,L,...
                          {dt,b0,H0,Gm0,R0});
  eks_rmse1 = sqrt(mean(sum((X(1:2,:)-SM_ERTS(1:2,:)).^2)));
  ME_ERTS = squeeze(SP_ERTS(1,1,:)+SP_ERTS(2,2,:));

  SV_ERTS = zeros(size(m,1),size(Y,2));
  SE_ERTS = zeros(size(m,1),size(Y,2));
  for k=1:size(Y,2)
    SV_ERTS(:,k) = diag(SP_ERTS(:,:,k));
    SE_ERTS(:,k) = (X(:,k) - SM_ERTS(:,k)).^2;
  end

  %
  % Smoother 2
  %
  [SM_ETF,SP_ETF] = etf_smooth1(MM_EKF,PP_EKF,Y,...
	func_df,Qc*dt,func_if,L,{dt,b0,H0,Gm0,R0},...
	func_dh,diag([vr va]),func_h,[],{xr,yr});

  eks_rmse2 = sqrt(mean(sum((X(1:2,:)-SM_ETF(1:2,:)).^2)));
  ME_ETF = squeeze(SP_ETF(1,1,:)+SP_ETF(2,2,:));
  fprintf('Done!\n');

%% Unscented Kalman filter

  fprintf('Running UKF...');

  % Initial values and space for (augmented) UKF
  m = m0;
  P = P0;
  MM_UKF = zeros(size(m,1),size(Y,2));
  PP_UKF = zeros(size(m,1),size(m,1),size(Y,2));
  VV_UKF = zeros(size(m,1),size(Y,2));
  EE_UKF = zeros(size(m,1),size(Y,2));

  % Filtering with UKF
  for k=1:size(Y,2)
    % Non-augmented UKF
    %[m,P] = ukf_predict1(m,P,func_a,Q,{dt,b0,H0,Gm0,R0});

    % Augmented UKF with separate sigma points
    %[m,P] = ukf_predict2(m,P,func_f,Qc*dt,{dt,b0,H0,Gm0,R0,L});

    %[m,P] = ukf_update2(m,P,Y(:,k),func_h,diag([vr va]),{xr,yr});

    % Augmented UKF with same sigma points for predict and update steps
    [m,P,X_s,w] = ukf_predict3(m,P,func_f,Qc*dt,diag([vr va]),d_param);
    [m,P] = ukf_update3(m,P,Y(:,k),func_h,diag([vr va]),X_s,w,h_param,h_param);

    MM_UKF(:,k) = m;
    PP_UKF(:,:,k) = P;
    VV_UKF(:,k) = diag(P);
    EE_UKF(:,k) = (X(:,k) - m).^2;
  end

  %
  % Calculate RMSE of UKF
  %
  ukf_rmse = sqrt(mean(sum((X(1:2,:)-MM_UKF(1:2,:)).^2)));
  ME_UKF = squeeze(PP_UKF(1,1,:)+PP_UKF(2,2,:));
  fprintf('Done!\n');

  fprintf('Running smoothers...');
  %
  % Smoother 1
  %
  [SM_URTS,SP_URTS] = urts_smooth1(MM_UKF,PP_UKF,func_f,Q,d_param);
  uks_rmse1 = sqrt(mean(sum((X(1:2,:)-SM_URTS(1:2,:)).^2)));
  ME_URTS = squeeze(SP_URTS(1,1,:)+SP_URTS(2,2,:));

  SV_URTS = zeros(size(m,1),size(Y,2));
  SE_URTS = zeros(size(m,1),size(Y,2));
  for k=1:size(Y,2)
    SV_URTS(:,k) = diag(SP_URTS(:,:,k));
    SE_URTS(:,k) = (X(:,k) - SM_URTS(:,k)).^2;
  end

  [SM_URTSb,SP_URTSb] = urts_smooth2(MM_UKF,PP_UKF,func_f,Qc*dt,d_param);
  uks_rmse1b = sqrt(mean(sum((X(1:2,:)-SM_URTSb(1:2,:)).^2)));
  ME_URTSb = squeeze(SP_URTSb(1,1,:)+SP_URTSb(2,2,:));

  %
  % Smoother 2
  %
  [SM_UTF,SP_UTF] = utf_smooth1(MM_UKF,PP_UKF,Y,...
	func_if,Qc*dt,d_param,...
	func_h,diag([vr va]),h_param);

  uks_rmse2 = sqrt(mean(sum((X(1:2,:)-SM_UTF(1:2,:)).^2)));
  ME_UTF = squeeze(SP_UTF(1,1,:)+SP_UTF(2,2,:));

  fprintf('Done!\n');

%% Gauss-Hermite Kalman filter

  fprintf('Running GHKF...');

  % Initial values and space for GHKF
  m = m0;
  P = P0;
  MM_GHKF = zeros(size(m,1),size(Y,2));
  PP_GHKF = zeros(size(m,1),size(m,1),size(Y,2));
  VV_GHKF = zeros(size(m,1),size(Y,2));
  EE_GHKF = zeros(size(m,1),size(Y,2));

  % Filtering with GHKF
  for k=1:size(Y,2)
    [m,P] = ghkf_predict(m,P,func_f,Q,{dt,b0,H0,Gm0,R0},3);
    [m,P] = ghkf_update(m,P,Y(:,k),func_h,diag([vr va]),{xr,yr},3);
    MM_GHKF(:,k) = m;
    PP_GHKF(:,:,k) = P;
    VV_GHKF(:,k) = diag(P);
    EE_GHKF(:,k) = (X(:,k) - m).^2;
  end

  %
  % Calculate RMSE
  %
  ghkf_rmse = sqrt(mean(sum((X(1:2,:)-MM_GHKF(1:2,:)).^2)));
  ME_GHKF = squeeze(PP_GHKF(1,1,:)+PP_GHKF(2,2,:));
  fprintf('Done!\n')

  fprintf('Running smoother...');
  %
  % Smoother
  %
  [SM_GHRTS,SP_GHRTS] = ghrts_smooth(MM_GHKF,PP_GHKF,func_f,Q,...
                          {dt,b0,H0,Gm0,R0},3);
  ghrts_rmse1 = sqrt(mean(sum((X(1:2,:)-SM_GHRTS(1:2,:)).^2)));
  ME_GHRTS = squeeze(SP_GHRTS(1,1,:)+SP_GHRTS(2,2,:));

  SV_GHRTS = zeros(size(m,1),size(Y,2));
  SE_GHRTS = zeros(size(m,1),size(Y,2));
  for k=1:size(Y,2)
    SV_GHRTS(:,k) = diag(SP_GHRTS(:,:,k));
    SE_GHRTS(:,k) = (X(:,k) - SM_GHRTS(:,k)).^2;
  end

  fprintf('Done!\n');

%% Cubature Kalman filter

  fprintf('Running CKF...');

  % Initial values and space for CKF
  m = m0;
  P = P0;
  MM_CKF = zeros(size(m,1),size(Y,2));
  PP_CKF = zeros(size(m,1),size(m,1),size(Y,2));
  VV_CKF = zeros(size(m,1),size(Y,2));
  EE_CKF = zeros(size(m,1),size(Y,2));

  % Filtering with CKF
  for k=1:size(Y,2)
    [m,P] = ckf_predict(m,P,func_f,Q,{dt,b0,H0,Gm0,R0});
    [m,P] = ckf_update(m,P,Y(:,k),func_h,diag([vr va]),{xr,yr});
    MM_CKF(:,k) = m;
    PP_CKF(:,:,k) = P;
    VV_CKF(:,k) = diag(P);
    EE_CKF(:,k) = (X(:,k) - m).^2;
  end

  %
  % Calculate RMSE
  %
  ckf_rmse = sqrt(mean(sum((X(1:2,:)-MM_CKF(1:2,:)).^2)));
  ME_CKF = squeeze(PP_CKF(1,1,:)+PP_CKF(2,2,:));
  fprintf('Done!\n')

  fprintf('Running smoother...');
  %
  % Smoother
  %
  [SM_CRTS,SP_CRTS] = crts_smooth(MM_CKF,PP_CKF,func_f,Q,...
                          {dt,b0,H0,Gm0,R0});
  crts_rmse1 = sqrt(mean(sum((X(1:2,:)-SM_CRTS(1:2,:)).^2)));
  ME_CRTS = squeeze(SP_CRTS(1,1,:)+SP_CRTS(2,2,:));

  SV_CRTS = zeros(size(m,1),size(Y,2));
  SE_CRTS = zeros(size(m,1),size(Y,2));
  for k=1:size(Y,2)
    SV_CRTS(:,k) = diag(SP_CRTS(:,:,k));
    SE_CRTS(:,k) = (X(:,k) - SM_CRTS(:,k)).^2;
  end

  fprintf('Done!\n');

%% Visualize methods

  if ~silent
    aa = 0.02*(-1:0.1:4);
    cx = R0 * cos(aa);
    cy = R0 * sin(aa);
    % Plot EKF estimate
    plot(xr,yr,'ko',cx,cy,'r-',X(1,:),X(2,:),'g-',...
         MM_EKF(1,:),MM_EKF(2,:),'k--');
    legend('Radar','Earth','True','Estimate');
    title('Filtering result with EKF');
    disp(' ');
    disp('Filtering result with EKF is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_1>');
    pause;

    % Error for x_1 with EKF
    semilogy(T,EE_EKF(1,:),'g-',T,VV_EKF(1,:),'b--',...
	     T,SE_ERTS(1,:),'r-',T,SV_ERTS(1,:),'k--');
    legend('EKF-RMSE','EKF-STDE',...
	   'ERTS-RMSE1','ERTS-STDE1');
    title('RMSE of estimating x_1 with EKF and ERTS')
    clc;
    disp('RMSE of estimating x_1 with EKF and ERTS is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_5>');
    pause

    % Error for x_5 with EKF
    semilogy(T,EE_EKF(5,:),'g-',T,VV_EKF(5,:),'b--',...
	     T,SE_ERTS(5,:),'r-',T,SV_ERTS(5,:),'k--');
    legend('EKF-RMSE','EKF-STDE',...
	   'ERTS-RMSE1','ERTS-STDE1');
    title('RMSE of estimating x_5 with EKF and ERTS')
    clc;
    disp('RMSE of estimating x_5 with EKF and ERTS is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_5>');
    pause

    % Plot UKF estimate
    plot(xr,yr,'ko',cx,cy,'r-',X(1,:),X(2,:),'g-',...
         MM_UKF(1,:),MM_UKF(2,:),'k--');
    legend('Radar','Earth','True','Estimate');
    clc;
    disp('Filtering result with UKF is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_1>');
    pause;

    % Error for x_1 with UKF
    semilogy(T,EE_UKF(1,:),'g-',T,VV_UKF(1,:),'b--',...
             T,SE_URTS(1,:),'r-',T,SV_URTS(1,:),'k--');
    legend('UKF-RMSE','UKF-STDE','URTS-RMSE','URTS-STDE');
    title('RMSE of estimating x_1 with UKF and URTS')
    clc;
    disp('RMSE of estimating x_1 with EKF and ERTS is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_5>');
    pause;

    % Error for x_5 with UKF
    semilogy(T,EE_UKF(5,:),'g-',T,VV_UKF(5,:),'b--',...
             T,SE_URTS(5,:),'r-',T,SV_URTS(5,:),'k--');
    legend('UKF-RMSE','UKF-STDE','URTS-RMSE','URTS-STDE');
    title('RMSE of estimating x_5 with UKF and URTS')
    clc;
    disp('RMSE of estimating x_5 with EKF and ERTS is now displayed.');

    % Plot GHKF estimate
    plot(xr,yr,'ko',cx,cy,'r-',X(1,:),X(2,:),'g-',...
         MM_GHKF(1,:),MM_GHKF(2,:),'k--');
    legend('Radar','Earth','True','Estimate');
    clc;
    disp('Filtering result with GHKF is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_1>');
    pause;

    % Error for x_1 with GHKF
    semilogy(T,EE_UKF(1,:),'g-',T,VV_UKF(1,:),'b--',...
             T,SE_URTS(1,:),'r-',T,SV_URTS(1,:),'k--');
    legend('GHKF-RMSE','GHKF-STDE','GHRTS-RMSE','GHRTS-STDE');
    title('RMSE of estimating x_1 with GHKF and GHRTS')
    clc;
    disp('RMSE of estimating x_1 with GHKF and GHRTS is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_5>');
    pause;

    % Error for x_5 with GHKF
    semilogy(T,EE_GHKF(5,:),'g-',T,VV_GHKF(5,:),'b--',...
             T,SE_GHRTS(5,:),'r-',T,SV_GHRTS(5,:),'k--');
    legend('GHKF-RMSE','GHKF-STDE','GHRTS-RMSE','GHRTS-STDE');
    title('RMSE of estimating x_5 with GHKF and GHRTS')
    clc;
    disp('RMSE of estimating x_5 with EKF and GHRTS is now displayed.');

    % Plot CKF estimate
    plot(xr,yr,'ko',cx,cy,'r-',X(1,:),X(2,:),'g-',...
         MM_CKF(1,:),MM_CKF(2,:),'k--');
    legend('Radar','Earth','True','Estimate');
    clc;
    disp('Filtering result with CKF is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_1>');
    pause;

    % Error for x_1 with CKF
    semilogy(T,EE_UKF(1,:),'g-',T,VV_UKF(1,:),'b--',...
             T,SE_URTS(1,:),'r-',T,SV_URTS(1,:),'k--');
    legend('CKF-RMSE','CKF-STDE','CRTS-RMSE','CRTS-STDE');
    title('RMSE of estimating x_1 with CKF and CRTS')
    clc;
    disp('RMSE of estimating x_1 with CKF and CRTS is now displayed.');
    disp(' ');
    disp('<press any key to see the estimation error of x_5>');
    pause;

    % Error for x_5 with CKF
    semilogy(T,EE_CKF(5,:),'g-',T,VV_CKF(5,:),'b--',...
             T,SE_CRTS(5,:),'r-',T,SV_CRTS(5,:),'k--');
    legend('CKF-RMSE','CKF-STDE','CRTS-RMSE','CRTS-STDE');
    title('RMSE of estimating x_5 with CKF and CRTS')
    clc;
    disp('RMSE of estimating x_5 with EKF and CRTS is now displayed.');

  end

%% Show RMSE errors for each method

  disp(' ');
  disp('RMS errors:');
  fprintf('EKF-RMSE   = %.6f [%.6f]\n',ekf_rmse,sqrt(mean(ME_EKF)));
  fprintf('ERTS-RMSE  = %.6f [%.6f]\n',eks_rmse1,sqrt(mean(ME_ERTS)));
  fprintf('ETF-RMSE   = %.6f [%.6f]\n',eks_rmse2,sqrt(mean(ME_ETF)));
  fprintf('UKF-RMSE   = %.6f [%.6f]\n',ukf_rmse,sqrt(mean(ME_UKF)));
  fprintf('URTS1-RMSE = %.6f [%.6f]\n',uks_rmse1,sqrt(mean(ME_URTS)));
  fprintf('URTS2-RMSE = %.6f [%.6f]\n',uks_rmse1b,sqrt(mean(ME_URTSb)));
  fprintf('UTF-RMSE   = %.6f [%.6f]\n',uks_rmse2,sqrt(mean(ME_UTF)));

  % Cubature methods
  fprintf('GHKF-RMSE  = %.6f [%.6f]\n',ghkf_rmse,sqrt(mean(ME_GHKF)));
  fprintf('GHRTS-RMSE = %.6f [%.6f]\n',ghrts_rmse1,sqrt(mean(ME_GHRTS)));
  fprintf('CKF-RMSE   = %.6f [%.6f]\n',ckf_rmse,sqrt(mean(ME_CKF)));
  fprintf('CRTS-RMSE  = %.6f [%.6f]\n',crts_rmse1,sqrt(mean(ME_CRTS)));

author: EEA-sensors

repo: ekfukf

save_path: github-repos/MATLAB/EEA-sensors-ekfukf


text:
function [net, info] = TrainModel(varargin)

opts.batchNormalization = false ;
opts.network = [] ;
opts.networkType = 'dagnn' ;
opts.PatchSize = 64;
opts.PatchNum = 50;
opts.Model = 'AlexNet';
opts.database = 'ChallengeDB_release';
opts.trainingpropertion = 0.8;
opts.seed = 1;
opts.Net = @NetIntial_S_CNN;
opts.quantizationMethod = 'uniform';
opts.bins = 5;
opts.beta = 10;
opts.PatchSize = 224;
opts.PatchStep = 16;
opts.epoch = 10;
[opts, varargin] = vl_argparse(opts, varargin) ;

sfx = opts.networkType ;
if opts.batchNormalization, sfx = [sfx '-bnorm'] ; end
[opts, varargin] = vl_argparse(opts, varargin) ;

if opts.bins == 1
    opts.lossType = 'MSE';
else
    opts.lossType = 'CE';
end

opts.expDir = fullfile('data',[opts.Model '_' opts.database '_TrainPropertion'...
              num2str(opts.trainingpropertion) '_PatchSize' num2str(opts.PatchSize)...
              '_PatchNum' num2str(opts.PatchNum) '_Seed' num2str(opts.seed),...
              '_bins' num2str(opts.bins), '_beta', num2str(opts.beta)]) ;
opts.dataDir = fullfile('databases', opts.database) ;
opts.imdbPath = fullfile('data', ['imdb_' opts.database '_PatchSize' num2str(opts.PatchSize)...
              '_PatchNum' num2str(opts.PatchNum) '.mat']);
opts.train = struct();
opts = vl_argparse(opts, varargin) ;
if ~isfield(opts.train, 'gpus'), opts.train.gpus = [1]; end;

% --------------------------------------------------------------------
%                                 model intialization
% --------------------------------------------------------------------

if isempty(opts.network)
    if strcmp(opts.Model,'AlexNet')
      net = load(fullfile('pretrained_models','imagenet-caffe-alex.mat'));
      net = vl_simplenn_tidy(net);
      net.layers = net.layers(1:end-2);
      net.layers{end+1} = struct('name','dropout1','type', 'dropout','rate', 0.5) ;
      net.layers{end+1} = struct('name','pred','type', 'conv', 'weights',...
                          {{randn(1, 1, 4096, opts.bins, 'single')*0.01,zeros(1,opts.bins,'single')}},...
                          'stride', 1, 'pad', 0) ;
      net = dagnn.DagNN.fromSimpleNN(net, 'canonicalNames', true) ;
      net.renameVar('input','data');
      net.addLayer('myLoss', dagnn.myLoss('lossType',opts.lossType), {'x21','label'}, 'objective') ;
      move(net, 'gpu')
      net.meta.inputSize = [opts.PatchSize opts.PatchSize 3] ;
      net.meta.trainOpts.learningRate = logspace(-3, -4, opts.epoch);
      net.meta.trainOpts.numEpochs = numel(net.meta.trainOpts.learningRate) ;
      net.meta.trainOpts.batchSize = 256 ;
      net.meta.trainOpts.numSubBatches = 1 ;
        net.meta.trainOpts.weightDecay = 0.0005;
    elseif strcmp(opts.Model,'ResNet')
        netStruct = load(fullfile('pretrained_models','imagenet-resnet-50-dag.mat')) ;
        net = dagnn.DagNN.loadobj(netStruct) ;
        net.removeLayer('fc1000');
        net.removeLayer('prob');
        dropoutBlock = dagnn.DropOut('rate',0.5);
        net.addLayer('dropout',dropoutBlock,{'pool5'},{'pool5d'},{});
        fc_size = [1 1 2048 opts.bins];
        fc5Block = dagnn.Conv('size',fc_size,'hasBias',true,'stride',[1,1],'pad',[0,0,0,0]);
        net.addLayer('pred',fc5Block,{'pool5d'},{'pred'},{'fc5_filter','fc5_bias'});
        f_ind = net.layers(end).paramIndexes(1);
        b_ind = net.layers(end).paramIndexes(2);
        he_gain = sqrt(2/(fc_size(1)*fc_size(2)*fc_size(3)));
        net.params(f_ind).value = (1/2)*he_gain*randn(fc_size, 'single');
        net.params(b_ind).value = zeros(fc_size(4),1, 'single');
        net.addLayer('myLoss', dagnn.myLoss('lossType',opts.lossType), {'pred','label'}, 'objective') ;
        move(net, 'gpu')
        net.meta.inputSize = {'input', [224,224,3,64]} ;
        net.meta.trainOpts.learningRate = logspace(-3, -4, opts.epoch);
        net.meta.trainOpts.numEpochs = opts.epoch;
        net.meta.trainOpts.batchSize = 64 ;
        net.meta.trainOpts.numSubBatches = 1 ;
        net.meta.trainOpts.weightDecay = 0.0001 ;
    elseif strcmp(opts.Model,'S_CNN')
        net = opts.Net('batchNormalization', opts.batchNormalization,...
                       'networkType', opts.networkType, 'lossType', opts.lossType,...
                       'class',opts.bins,'epoch',opts.epoch) ;
    end
end

%% for pre-processing, no use currently
filter_size = 7;
opts.filter_ = gpuArray(fspecial('gaussian',filter_size,1.0));

% --------------------------------------------------------------------
%                                 dataset  preparison
% --------------------------------------------------------------------

if exist(opts.imdbPath, 'file')
  imdb = load(opts.imdbPath) ;
else
    switch opts.database
        case 'LIVE'
            imdb = getImdbLIVE(opts) ;
        case 'TID2013'
            imdb = getImdbTID2013(opts) ;
        case 'CSIQ'
            imdb = getImdbCSIQ(opts) ;
        case 'ChallengeDB_release'
            imdb = getImdbChalleng(opts) ;
        otherwise
            error('Unknown database!');
    end

    if opts.bins > 1
        if strcmp(opts.quantizationMethod,'LloydMax')
            [partition,codebook] = lloyds(imdb.images.score,opts.bins);
            imdb.weights = codebook;
            D = bsxfun(@minus,imdb.images.score,imdb.weights');
            D = exp(-opts.beta*(D.^2));
            D = bsxfun(@rdivide,D,sum(D));
            imdb.images.labels = D;
        else
            partition = 1.0 / opts.bins;
            imdb.weights = partition/2:partition:1.001;
            D = bsxfun(@minus,imdb.images.score,imdb.weights');
            D = exp(-opts.beta*(D.^2));
            D = bsxfun(@rdivide,D,sum(D));
            imdb.images.labels = D;
        end
        scores_train = imdb.images.score(1:opts.PatchNum:end)';
        D = D(:,1:opts.PatchNum:end)';
        svm_model = svmtrain(scores_train, D, '-s 4 -t 0');
        [scores_pred, ~,~] = svmpredict(scores_train, D, svm_model,'-q');
        svm_err = (sum(abs(scores_pred - scores_train))) / length(scores_train);
    else
        imdb.images.labels = imdb.images.score;
    end
end
dataMean = mean(imdb.images.data, 4);
imdb.images.data_mean = dataMean;

try
    net.meta.classes.weights = imdb.weights;
    net.meta.svm_model = svm_model;
    net.meta.svm_err = svm_err;
end
net.meta.normalization.averageImage = dataMean;

% --------------------------------------------------------------------
%                          model training
% --------------------------------------------------------------------

switch opts.networkType
  case 'simplenn', trainfn = @cnn_train_modified ;
  case 'dagnn', trainfn = @cnn_train_dag_modified ;
end

[net, info] = trainfn(net, imdb, getBatch(opts), ...
  'expDir', opts.expDir, ...
  net.meta.trainOpts, ...
  opts.train, ...
  'val', find(imdb.images.set == 3)) ;

% --------------------------------------------------------------------
function fn = getBatch(opts)
% --------------------------------------------------------------------
switch lower(opts.networkType)
  case 'simplenn'
    fn = @(x,y) getSimpleNNBatch(x,y) ;
  case 'dagnn'
    bopts = struct('numGpus', numel(opts.train.gpus)) ;
    fn = @(x,y) getDagNNBatch(bopts,x,y) ;
end

% --------------------------------------------------------------------
function inputs = getDagNNBatch(opts, imdb, batch)
% --------------------------------------------------------------------
images = imdb.images.data(:,:,:,batch) ;
images = bsxfun(@minus, single(images), imdb.images.data_mean) ;
labels = single(imdb.images.labels(:,batch));
if opts.numGpus > 0
  images = gpuArray(images) ;
end
inputs = {'data', images, 'label', labels} ;

author: HuiZeng


text:
function f = floor(f)
%FLOOR   Pointwise floor function of a CLASSICFUN.
%   G = FLOOR(F) returns the CLASSICFUN G such that G(X) = FLOOR(F(x)) for each x in
%   F.domain.
%
%   If F is complex, then the G = FLOOR(REAL(F)) + 1i*FLOOR(IMAG(F)).
%
%   Note that FLOOR() assumes the output G(X) is a constant. If it is not, then
%   garbage is returned with no warning.
%
% See also CEIL, ROUND, FIX.

% Copyright 2017 by The University of Oxford and The Chebfun Developers.
% See http://www.chebfun.org/ for Chebfun information.

% FLOOR() the ONEFUN:
f.onefun = floor(f.onefun);

end

author: chebfun

repo: chebfun

save_path: github-repos/MATLAB/chebfun-chebfun


text:
function cropDose(doseNum,structureNum,margin)
%function cropDose(doseNum,structureNum,margin)
%
%This function crops the dose doseNum to lie within the structure
%structureNum plus margin.
%
%Usage:
% doseNum      = 1; %final
% structureNum = 8; %skin
% margin       = 2; %cm
% cropDose(doseNum,structureNum,margin)
%
%APA, 07/03/2010
%
% Copyright 2010, Joseph O. Deasy, on behalf of the CERR development team.
%
% This file is part of The Computational Environment for Radiotherapy Research (CERR).
%
% CERR development has been led by:  Aditya Apte, Divya Khullar, James Alaly, and Joseph O. Deasy.
%
% CERR has been financially supported by the US National Institutes of Health under multiple grants.
%
% CERR is distributed under the terms of the Lesser GNU Public License.
%
%     This version of CERR 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.
%
% CERR 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 CERR.  If not, see <http://www.gnu.org/licenses/>.

global stateS planC
indexS = planC{end};

% for command line help document
if ~exist('doseNum') || ~exist('structureNum') || ~exist('margin')
    prompt = {'Enter the dose number';'Enter the structure number'; 'Enter the margin (cm)'};
    dlg_title = 'Crop scan';
    num_lines = 1;
    def = {'';'';''};
    outPutQst = inputdlg(prompt,dlg_title,num_lines,def);
    if isempty(outPutQst{1}) || isempty(outPutQst{2}) || isempty(outPutQst{3})
        warning('Need to enter all the inputs');
        return
    else
        doseNum         = str2num(outPutQst{1});
        structureNum    = str2num(outPutQst{2});
        margin          = str2num(outPutQst{3});
    end
end

%Get associated scan number
scanNum = getStructureAssociatedScan(structureNum);
assocScanUID = planC{indexS.dose}(doseNum).assocScanUID;
scanNumDose = getAssociatedScan(assocScanUID);
tmDose = getTransM('dose',doseNum,planC);
tmScan = getTransM('scan',scanNum,planC);
if isempty(scanNumDose) && ~isequal(tmDose,tmScan)
    error('This function currently supports dose and structure with same transformation matrix')
end

%Get dose grid
[xDoseVals, yDoseVals, zDoseVals] = getDoseXYZVals(planC{indexS.dose}(doseNum));

%Get structure boundary
rasterSegments = getRasterSegments(structureNum,planC);
zMin = min(rasterSegments(:,1)) - margin;
zMax = max(rasterSegments(:,1)) + margin;
xMin = min(rasterSegments(:,3)) - margin;
xMax = max(rasterSegments(:,4)) + margin;
yMin = min(rasterSegments(:,2)) - margin;
yMax = max(rasterSegments(:,2)) + margin;

%Get min, max indices of DoseArray
[JLow, jnk] = findnearest(xDoseVals, xMin);
[jnk, JHigh] = findnearest(xDoseVals, xMax);
[jnk, ILow] = findnearest(yDoseVals, yMin);
[IHigh, jnk] = findnearest(yDoseVals, yMax);
[KLow, jnk] = findnearest(zDoseVals, zMin);
[jnk, KHigh] = findnearest(zDoseVals, zMax);

%Crop DoseArray
planC{indexS.dose}(doseNum).doseArray = planC{indexS.dose}(doseNum).doseArray(IHigh:ILow,JLow:JHigh,KLow:KHigh);

%Reassign zvalues
planC{indexS.dose}(doseNum).zValues = planC{indexS.dose}(doseNum).zValues(KLow:KHigh);

%Reassign coordinates of left corner
planC{indexS.dose}(doseNum).coord1OFFirstPoint = xDoseVals(JLow);
planC{indexS.dose}(doseNum).coord2OFFirstPoint = yDoseVals(IHigh);

%Reassign Dimension
planC{indexS.dose}(doseNum).sizeOfDimension2 = length(IHigh:ILow);
planC{indexS.dose}(doseNum).sizeOfDimension1 = length(JLow:JHigh);
planC{indexS.dose}(doseNum).sizeOfDimension3 = length(KLow:KHigh);

%Create new UID since this dose has changed
planC{indexS.dose}(doseNum).doseUID = createUID('dose');

stateS.doseSet = doseNum;
stateS.doseSetChanged = 1;

sliceCallBack('refresh');

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [ images_array ] = readImages( image_dir, max )
%readImages Summary
% Reads images in the give directory and returns a cell array.

% Read images
file_names = fullfile(image_dir, '*.jpg');
image_files = dir(file_names);

% Number of frames for testing
n = min([length(image_files) max]);

% Preallocation
images_array = cell(n, 1);

for k = 1:n
    base_name = image_files(k).name;
    full_name = fullfile(image_dir, base_name);
    fprintf(1, 'Now reading %s\n', full_name);
    im = imread(full_name);
    im = imresize(im, [360 640]);
    images_array{k} = im;
end

end

author: ishit


text:
function [vol, ext, endian] = er_ldtvolume(stem,ext)
% [vol ext endian] = er_ldtvolume(stem,<ext>)
%
% Loads a volume in tSeries format and returns a 4D structure
% of dimension Nslices X Nrows X Ncols X Ndepth (where Ndepth = # of
% temporal frames).
%
% stem should be the path where a bunch of tSeries are stored. E.g.:
% stem = 'X:\mri\mySession\Inplane\Original\tSeries\Scan1'. Within this
% directory the code looks for files named 'tSeries[#].mat', the default
% format for mrLoadRet.
%
% 'ext' is currently a dummy param, but since some older tSeries are saved
% as .dat files, maybe I'll make it work.
%
% '$Id: er_ldtvolume.m,v 1.1 2004/03/11 01:29:55 sayres Exp $'
%
% 06/18/03 ras: created from fmri_ldbvolume in an attempt to integrate
% fsfast into mrLoadRet.

vol = [];
ext = 'tSeries';
endian = 0;

if(nargin ~= 1 & nargin ~= 2)
  fprintf(2,'USAGE: er_ldtvolume(stem,<ext>)\n');
  qoe; error; return;
end

% First, check whether tSeries are in stem path %
test = filterdir('tSeries',stem);
if isempty(test)
    fprintf(2,'Couldn''t find any tSeries in specified path!!!');
    qoe; error; return;
end

% count slices
nslices = size(test);

% fprintf('er_ldtvolume: found %d slices\n',nslices);

% concatenate into volume
for slice = 1:nslices
  % fprintf('Loading Slice %3d\n',slice);
  fname = fullfile(stem,['tSeries' num2str(slice) '.mat']);

  if ~exist(fname,'file')
      warning(['Could not find file ', fname]);
  end

  z = er_ldtfile(fname);
  if(size(z,1) == 1 & size(z,3) == 1)
    vol(slice,:,:)  = z;
  else
    vol(slice,:,:,:) = z;
  end
end

return;

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
% This demo efficiently computes the emd_hat between
% histograms where the ground distance is non
% symmetric and the histograms are of non equal size.
% emd_hat is described in the paper:
% A Linear Time Histogram Metric for Improved SIFT Matching
% Ofir Pele, Michael Werman
%  ECCV 2008
% The efficient algorithm is described in the paper:
%  Fast and Robust Earth Mover's Distances
%  Ofir Pele, Michael Werman
%  ICCV 2009
clc; close all; clear all;

D= [0 3 5;
    10 0 30;
    90 80 0;
   80 80 80];

P= [1; 0; 0; 2];
Q= [0; 1; 1];

[dist F]= emd_hat_mex_nes(P,Q,D,0,3)
[dist F]= emd_hat_mex_nes(Q,P,D',0,3)

% Copyright (c) 2009-2012, Ofir Pele
% All rights reserved.

% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are
% met:
%    * Redistributions of source code must retain the above copyright
%    notice, this list of conditions and the following disclaimer.
%    * Redistributions in binary form must reproduce the above copyright
%    notice, this list of conditions and the following disclaimer in the
%    documentation and/or other materials provided with the distribution.
%    * Neither the name of the The Hebrew University of Jerusalem nor the
%    names of its contributors may be used to endorse or promote products
%    derived from this software without specific prior written permission.

% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
% IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
% THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
% PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR
% CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
% EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
% PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
% PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
% LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
% NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
% SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

author: cvzoya

repo: saliency

save_path: github-repos/MATLAB/cvzoya-saliency