text:
function SmoothY=fastsmooth(Y,w,type,ends)
% fastsmooth(Y,w,type,ends) smooths vector Y with smooth
%  of width w. Version 3.0, October 2016.
% The argument "type" determines the smooth type:
%   If type=1, rectangular (sliding-average or boxcar)
%   If type=2, triangular (2 passes of sliding-average)
%   If type=3, pseudo-Gaussian (3 passes of sliding-average)
%   If type=4, pseudo-Gaussian (4 passes of same sliding-average)
%   If type=5, multiple-width (4 passes of different sliding-average)
% The argument "ends" controls how the "ends" of the signal
% (the first w/2 points and the last w/2 points) are handled.
%   If ends=0, the ends are zero.  (In this mode the elapsed
%     time is independent of the smooth width). The fastest.
%   If ends=1, the ends are smoothed with progressively
%     smaller smooths the closer to the end. (In this mode the
%     elapsed time increases with increasing smooth widths).
% fastsmooth(Y,w,type) smooths with ends=0.
% fastsmooth(Y,w) smooths with type=1 and ends=0.
% Examples:
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3)= [0 1 4 7 10 7 4 1 1 0]
%
% fastsmooth([1 1 1 10 10 10 1 1 1 1],3,1,1)= [1 1 4 7 10 7 4 1 1 1]
%
% x=1:100;
% y=randn(size(x));
% plot(x,y,x,fastsmooth(y,5,3,1),'r')
% xlabel('Blue: white noise.    Red: smoothed white noise.')
%
% Copyright (c) 2012, Thomas C. O'Haver
%
% 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.
if nargin==2, ends=0; type=1; end
if nargin==3, ends=0; end
  switch type
    case 1
       SmoothY=sa(Y,w,ends);
    case 2
       SmoothY=sa(sa(Y,w,ends),w,ends);
    case 3
       SmoothY=sa(sa(sa(Y,w,ends),w,ends),w,ends);
    case 4
       SmoothY=sa(sa(sa(sa(Y,w,ends),w,ends),w,ends),w,ends);
    case 5
       SmoothY=sa(sa(sa(sa(Y,round(1.6*w),ends),round(1.4*w),ends),round(1.2*w),ends),w,ends);
  end

function SmoothY=sa(Y,smoothwidth,ends)
w=round(smoothwidth);
SumPoints=sum(Y(1:w));
s=zeros(size(Y));
halfw=round(w/2);
L=length(Y);
for k=1:L-w,
   s(k+halfw-1)=SumPoints;
   SumPoints=SumPoints-Y(k);
   SumPoints=SumPoints+Y(k+w);
end
s(k+halfw)=sum(Y(L-w+1:L));
SmoothY=s./w;
% Taper the ends of the signal if ends=1.
  if ends==1,
  startpoint=(smoothwidth + 1)/2;
  SmoothY(1)=(Y(1)+Y(2))./2;
  for k=2:startpoint,
     SmoothY(k)=mean(Y(1:(2*k-1)));
     SmoothY(L-k+1)=mean(Y(L-2*k+2:L));
  end
  SmoothY(L)=(Y(L)+Y(L-1))./2;
  end

author: CHLNDDEV


text:
%
% Calculates the time necessary for the calculation of the correlation dimension.
%
%
%
%
org = [2];
startfd;
radm=[];
rasm=[];
tim=[];
nvt=[];

for m= 250:250:3750

    E=ran(1:m,:);
    tic;
    pdc3nofig;
    tim=toc;
    nvt =[nvt; tim];
    E = ran;

end

figure;
ax = gca;
plot(m,nvt,'ko');
xlabel('Number of Points', 'fontsize',12);
ylabel('Seconds', 'fontsize',12);
Title('D-value Computation Time', 'fontsize',14);
set(ax,'fontsize',11);

author: CelsoReyes


text:
function obj = t3_hand_eye_func_new(in1,in2,in3)
coefs_tq1_1 = in2(1);
coefs_tq1_2 = in2(4);
coefs_tq1_3 = in2(7);
coefs_tq1_4 = in2(10);
coefs_tq1_5 = in2(13);
coefs_tq1_6 = in2(16);
coefs_tq1_7 = in2(19);
coefs_tq1_8 = in2(22);
coefs_tq1_9 = in2(25);
coefs_tq2_1 = in2(2);
coefs_tq2_2 = in2(5);
coefs_tq2_3 = in2(8);
coefs_tq2_4 = in2(11);
coefs_tq2_5 = in2(14);
coefs_tq2_6 = in2(17);
coefs_tq2_7 = in2(20);
coefs_tq2_8 = in2(23);
coefs_tq2_9 = in2(26);
coefs_tq3_1 = in2(3);
coefs_tq3_2 = in2(6);
coefs_tq3_3 = in2(9);
coefs_tq3_4 = in2(12);
coefs_tq3_5 = in2(15);
coefs_tq3_6 = in2(18);
coefs_tq3_7 = in2(21);
coefs_tq3_8 = in2(24);
coefs_tq3_9 = in2(27);
coefs_tq1_10 = in2(28);
coefs_tq1_11 = in2(31);
coefs_tq2_10 = in2(29);
coefs_tq2_11 = in2(32);
coefs_tq3_10 = in2(30);
coefs_tq3_11 = in2(33);
pinvG3_1 = in1(3);
pinvG3_2 = in1(6);
pinvG3_3 = in1(9);
q0 = in3(1,:);
q1 = in3(2,:);
q2 = in3(3,:);
q3 = in3(4,:);
obj = q0.^2.*(coefs_tq1_1.*pinvG3_1+coefs_tq2_1.*pinvG3_2+coefs_tq3_1.*pinvG3_3)+q1.^2.*(coefs_tq1_5.*pinvG3_1+coefs_tq2_5.*pinvG3_2+coefs_tq3_5.*pinvG3_3)+q2.^2.*(coefs_tq1_8.*pinvG3_1+coefs_tq2_8.*pinvG3_2+coefs_tq3_8.*pinvG3_3)+q3.^2.*(coefs_tq1_10.*pinvG3_1+coefs_tq2_10.*pinvG3_2+coefs_tq3_10.*pinvG3_3)+coefs_tq1_11.*pinvG3_1+coefs_tq2_11.*pinvG3_2+coefs_tq3_11.*pinvG3_3+q0.*q1.*(coefs_tq1_2.*pinvG3_1+coefs_tq2_2.*pinvG3_2+coefs_tq3_2.*pinvG3_3)+q0.*q2.*(coefs_tq1_3.*pinvG3_1+coefs_tq2_3.*pinvG3_2+coefs_tq3_3.*pinvG3_3)+q0.*q3.*(coefs_tq1_4.*pinvG3_1+coefs_tq2_4.*pinvG3_2+coefs_tq3_4.*pinvG3_3)+q1.*q2.*(coefs_tq1_6.*pinvG3_1+coefs_tq2_6.*pinvG3_2+coefs_tq3_6.*pinvG3_3)+q1.*q3.*(coefs_tq1_7.*pinvG3_1+coefs_tq2_7.*pinvG3_2+coefs_tq3_7.*pinvG3_3)+q2.*q3.*(coefs_tq1_9.*pinvG3_1+coefs_tq2_9.*pinvG3_2+coefs_tq3_9.*pinvG3_3);

author: zarathustr

repo: LibQPEP

save_path: github-repos/MATLAB/zarathustr-LibQPEP


text:
function title = p06_title ( )

%*****************************************************************************80
%
%% P06_TITLE returns a title for problem 6.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    15 February 2012
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Harald Niederreiter, Kevin McCurley,
%    Optimization of functions by quasi-random search methods,
%    Computing,
%    Volume 22, Number 2, 1979, pages 119-123.
%
%  Parameters:
%
%    Output, string TITLE, a title for the problem.
%
  title = 'f(x) = - sin(1/x(1)+1/x(2)+1/x(3)+1/x(4))';

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% center distance to closest larva
function [data,units] = compute_dcentral(trx,n)

larvae = trx.exp2flies{n};
nlarvae = numel(larvae);
data = cell(1,nlarvae);

for i1 = 1:nlarvae,
  larva1 = larvae(i1);
  % access closestlarva to ensure that dcentral is computed
  trx(larva1).closestlarva_central;
  data{i1} = trx(larva1).dcentral;
end
units = parseunits('mm');

author: kristinbranson

repo: JAABA


text:
function PQ = paddedsize(AB,CD,PARAM)
%PADDEDSIZE Computes padded sizes useful for FFT-based filtering.
%   PQ = PADDEDSIZE(AB), where AB is a two-element vector, computes the
%   two-element vector PQ = 2*AB.
%
%   PQ = PADDEDSIZE(AB,'PWR2') computes the vector PQ such that PQ(1) =
%   PQ(2) = 2^nextpow2(2*m), where m is MAX(AB).
%
%   PQ = PADDEDSIZE(AB,CD), where AB and CD are two-element size vectors,
%   computes the two-element size vector PQ. The elements of PQ are the
%   smallest even integers greater than or equal to AB + CD - 1.
%
%   PQ = PADDEDSIZE(AB,CD,'PWR2') computes the vector PQ such that PQ(1)
%   = PQ(2) = 2^nextpow2(2*m), where m is MAX([AB,CD]).
%
%   Copyright 2002-2020 Gatesmark
%
%   This function, and other functions in the DIPUM Toolbox, are based
%   on the theoretical and practical foundations established in the
%   book Digital Image Processing Using MATLAB, 3rd ed., Gatesmark
%   Press, 2020.
%
%   Book website: http://www.imageprocessingplace.com
%   License: https://github.com/dipum/dipum-toolbox/blob/master/LICENSE.txt

if nargin == 1
   PQ  = 2*AB;
elseif nargin == 2 && ~ischar(CD)
   PQ = AB + CD - 1;
   PQ = 2*ceil(PQ/2);
elseif nargin == 2
   m = max(AB); % Maximum dimension.

   % Find power-of-2 at least twice m.
   P = 2^nextpow2(2*m);
   PQ = [P,P];
elseif (nargin == 3) && strcmpi(PARAM,'pwr2')
   m = max([AB,CD]); % Maximum dimension.
   P = 2^nextpow2(2*m);
   PQ = [P,P];
else
   error('Wrong number of inputs.')
end

author: dipum

repo: dipum-toolbox

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


text:
%%
% Purpose:
%
% This algorithm decrypts and encrypts images based on keys
%
% INPUTS:
%
%
%
% imMat: If the nargin is one, imMat is the filename of the image which
% you would like to encrypt. If the nargin is 4, imMat is the encrypted
% matrix
%
% key1, key2, key3: These are number values used to encrypt and decrypt
% images.
% key3 and key2 and key1 can be left blank if the nargin is one and  you would like
% the computer to make the keys for you. if the nargin is 4 key1 is
% a vector at baseto that is the encrypted vector of the original random numbers.
% key2 is the base of the original vector (baseKey) **Please
% read the help section of baseVecExpand.
%
% OUTPUTS:
%
% cryptPic: If the nargin = 1, cryptPic is a matrix of encrypted
% numbers representing an encrypted image. If the nargin = 4, imMat
% is a uint8 array which is a decrypted image. This function will also
% display the decrypted image.
%
% keyVec: If the nargin = 1, the computer will use the random number generator to
% make three keys. These three keys are then changed to a random base, baseto
% using baseVecExpand. The starting base is baseKey.
% keyVec is a vector containing the  key vector at baseto. if nargin = 4,
% keyVec is not returned
%
% baseKey: If the nargin = 1 baseKey is the base of the
% original vector. please read the help section of baseVecExpand. baseKey
% is not returned if the nargin = 4.
%
% COMMENTS:
%
% Please change the formulas used to encrypt and decrypt the images into
% your own unique formulas so that nobody can see your images.
function [cryptPic keyVec baseto baseKey] = picCrypt(imMat,key,key2,key3)
if nargin == 1
   key = rand;
   key2 = rand;
   key3 = rand;
   baseto = rand*1000;
   baseKey = max(floor(1000000*[key key2 key3])) + rand*100;
   keyVec = baseVecExpand(floor(1000000*[key key2 key3]),baseKey,baseto);
   pic = 10000*im2double(imread(imMat));
   cryptPic = ((floor(pic)/key2).^key)+key3;
   imagesc(im2uint8(cryptPic));
   axis equal;
   axis tight;
elseif nargin == 4
   key10 = baseVecExpand(key,key3,key2)/1000000;
   cryptPic = im2uint8(((imMat - key10(3)).^(1/key10(1))*key10(2))/10000); % CHANGE THIS FORMULA (this is the opposite of the formula above
   imagesc(cryptPic);
   axis equal;
   axis tight;
else
    error('You must enter the correct number of parameters');
end

author: Sable

repo: mcbench-benchmarks


text:
%% shiftV
% Below is a demonstration of the features of the |shiftV| function

%%
clear; close all; clc;

%% Syntax
% |[Vv]=shiftV(V,v);|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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 v1 = uminus(v1)
% overloads unitary minus

v1.x = -v1.x;
v1.y = -v1.y;
v1.z = -v1.z;

author: mtex-toolbox

repo: mtex

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


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_clutter';

RadarCoords.lat = 50.346069;
RadarCoords.lon = -4.113670;

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

%% Focus on region of high clutter

% Convert LLA limits to NED
[y_lim, x_lim] = geodetic2ned(lat_lim,...
                              lon_lim, ...
                              0,...
                              RadarCoords.lat,...
                              RadarCoords.lon,...
                              0,...
                              referenceEllipsoid('wgs84'));
V_bounds = [x_lim(1);... % X-min |
            x_lim(2);... % X-max | Surveillance region bounding
            y_lim(1);... % Y-min | box coordinates (m)
            y_lim(2)]';  % Y-max |

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

% Recording settings
Record = 0;                 % 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,'_mon.avi'); % Set to the desired path and name of produced recording

% Model parameter shortcuts
lambdaV = 3; % Expected number of clutter measurements over entire surveillance region
V = (abs(V_bounds(2)-V_bounds(1))*abs(V_bounds(4)-V_bounds(3))); % Total area of surveillance region
P_D = 0.5;    % Probability of detection
timestep_duration = duration(0,0,2);

%% Instantiation of necessary components

% Instantiate a Dynamic model (CV model with q = 1 m/s^2)
transition_model = OrnsteinUhlenbeckModelX('NumDims',2,...
                             'VelocityErrVariance',0.1,...
                             'DampingCoefficient',0.01,...
                             'TimestepDuration',timestep_duration);
% Instantiate an Observation model (Variance of 50m^2 on each coordinate)
measurement_model = LinearGaussianX('NumMeasDims', 2,...
                                    'NumStateDims', 4,...
                                    'MeasurementErrVariance', 10^2,...
                                    'Mapping', [1 3]);
clutter_model = PoissonRateUniformPosition2X('ClutterRate',lambdaV,...
                                            'Limits',[V_bounds(1:2);...
                                                      V_bounds(3:4)]);
detection_model = ConstantDetectionProbabilityX('DetectionProbability',P_D);

birth_model = DistributionBasedBirthModelX('Distribution', UniformDistributionX([V_bounds(1:2); ...
                                                                                [-7 7 ];...
                                                                                V_bounds(3:4);...
                                                                                [-7 7 ]]),...
                                           'BirthIntensity', 1);

% Compile the State-Space model
model = StateSpaceModelX(transition_model, measurement_model,...
                        'Clutter',clutter_model,...
                        'Detection', detection_model,...
                        'Birth', birth_model);

%% Base Filter
obs_covar= measurement_model.covar();
PriorState = GaussianStateX(zeros(4,1), transition_model.covar() + blkdiag(obs_covar(1,1), 5, obs_covar(2,2), 5));
base_filter = ExtendedKalmanFilterX('Model', model, 'StatePrior', PriorState);

%% Data Associator
config.ClutterModel = clutter_model;
config.Clusterer = NaiveClustererX();
config.Gater = EllipsoidalGaterX(2,'GateLevel',10)';
config.DetectionModel = detection_model;
jpdaf = GlobalNearestNeighbourDataAssocX(config);
% jpdaf = JointProbabilisticDataAssocX(config);

%% Track Initiator

% Initiate Data Associator
config_pdaf.ClutterModel = clutter_model;
% config.Clusterer = NaiveClustererX();
config_pdaf.Gater = EllipsoidalGaterX(2,'GateLevel',5);
config_pdaf.DetectionProbability = P_D;
% pdaf = ProbabilisticDataAssocX(config_pdaf);
pdaf = NearestNeighbourDataAssocX(config_pdaf);

% Initiate Tag Generator
tag_gen = RandSampleTagGeneratorX(1:1000000);

% Prepare initiator parameters
config_ti.TagGenerator = tag_gen;
config_ti.InitFilter = base_filter;
config_ti.DataAssociator = pdaf;
config_ti.ConfirmThreshold = [4,10];
config_ti.DeleteThreshold = [8,10];
% CovarThreshold = 4*PriorState.Covar;
% config_ti.CustomDeleteConditionFcn = ...
%     @(x,t) t.Filter.StatePosterior.Covar(1,1)>CovarThreshold(1,1) ...
%            || t.Filter.StatePosterior.Covar(3,3)>CovarThreshold(3,3);

% Create the track initiator
myti = MofN_TrackInitiatorX(config_ti);

%% Track Deleter
config_td = struct('Fieldname', 'ExistenceProbability',...
                   'ReferenceValue', 0.1,...
                   'ReferenceOperand', 'lt');
mytd = FieldBasedDeleterX(config_td);

TrackList = {};
DeletedTracks = {};

%% Create plot windows
if(ShowPlots)

    % Map plot
    figure('units','normalized','outerposition',[0 0 .5 1])
    ax(1) = gca;
    if PlotLatLon
        plot_google_map('Axis',ax(1),'APIKey','AIzaSyBXKujdtXRZiqya1soVS9pxBzYR4g7aGvM','Resize',3,'Scale',2,'MapType','satellite');
        axis(ax(1),[lon_lim(1) lon_lim(2) lat_lim(1) lat_lim(2)])
    end
%     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=2:N
    fprintf('Iteration = %d/%d\n================>\n',k,N);

    %% Extract DataList at time k
    MeasurementList = MeasurementScans(k);
    timestamp_km1 = MeasurementScans(k-1).Timestamp;
    timestamp_k = MeasurementList.Timestamp;
    dt = timestamp_k - timestamp_km1;
    transition_model.TimestepDuration = dt;
    fprintf('Timestamp = %s\n================>\n',timestamp_k);

    %% Process JPDAF
    jpdaf.MeasurementList = MeasurementList;
    jpdaf.TrackList = TrackList;
    jpdaf.predictTracks();
    jpdaf.associate(TrackList, MeasurementList);
    jpdaf.updateTracks();

%     if(MeasurementList.NumMeasurements>0)
%         % Convert measurements to LLA and plot them
%         [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
%                                    MeasurementList.Vectors(1,:),...
%                                    0,...
%                                    RadarCoords.lat,...
%                                    RadarCoords.lon,...
%                                    0,...
%                                    referenceEllipsoid('wgs84'));
%         plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
%
%         for i=1:numel(lat)
%             plots(end+1) = text(ax(1), lon(:,i)+0.0001,lat(:,i)+0.00027,num2str(i),'FontSize',8,'Color','w');
%         end
% %             plot(ax(1), RadarCoords.lon,RadarCoords.lat,...
% %                                '-s','MarkerSize',20,...
% %                                'MarkerEdgeColor','red',...
% %                                'MarkerFaceColor',[1 .6 .6]);
%     end
%     pause(0.01);
     %% Perform Track initiation
    [TrackList, ~, del_tracks] = myti.initiateTracks(jpdaf.TrackList, MeasurementList, jpdaf.AssocWeightsMatrix);
    DeletedTracks = [DeletedTracks, del_tracks];

    del_tracks = [];
    for t=1:numel(TrackList)
        track = TrackList{t};
        xk = track.State.Mean;
        [ylim, xlim] = geodetic2ned(lat_lim,...
                                lon_lim, ...
                                0,...
                                RadarCoords.lat,...
                                RadarCoords.lon,...
                                0,...
                                referenceEllipsoid('wgs84'));
        if xk(1,:)>x_lim(2) || xk(1,:)<x_lim(1) || xk(3,:)>y_lim(2) || xk(3,:)<y_lim(1)
            del_tracks{end+1} = track;
            TrackList{t} = [];
        end
    end
    TrackList = TrackList(~cellfun('isempty',TrackList));
    DeletedTracks = [DeletedTracks, del_tracks];
    %% Perform Track deletion
%     TrackList = mytd.deleteTracks(TrackList);

     % 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(MeasurementList.NumMeasurements>0)
            % Convert measurements to LLA and plot them
            [lat,lon,~] = ned2geodetic(MeasurementList.Vectors(2,:),...
                                       MeasurementList.Vectors(1,:),...
                                       0,...
                                       RadarCoords.lat,...
                                       RadarCoords.lon,...
                                       0,...
                                       referenceEllipsoid('wgs84'));
            plots(end+1) = plot(ax(1), lon,lat,'y*','MarkerSize', 10);
        end

        % Plot region of low PD
        x = [-3422.85261310741,-3191.58441229017,-2993.55608122595,-3228.08558459284, -3422.85261310741];
        y = [1.472966334316646e+03,1.123217741935625e+03,1.243952901078573e+03,1.579667558371468e+03, 1.472966334316646e+03];
        [y,x,~] = ned2geodetic(y,...
                               x,...
                               0,...
                               RadarCoords.lat,...
                               RadarCoords.lon,...
                               0,...
                               referenceEllipsoid('wgs84'));
        plots(end+1) = plot(ax(1), x, y, 'm-');

        % Plot region of high clutter
        x = [-469.710602736631,-2572.92295458513,-2516.31234127506,-1320.74163107899,-835.307363807233, -469.710602736631];
        y = [-4915.41679013513,-4065.50122858976,-2858.60918927224,-1380.19032924249,-1261.08679763585, -4915.41679013513];
        [y,x,~] = ned2geodetic(y,...
                               x,...
                               0,...
                               RadarCoords.lat,...
                               RadarCoords.lon,...
                               0,...
                               referenceEllipsoid('wgs84'));
        plots(end+1) = plot(ax(1), x, y, 'r-');

        % Plot all existing tracks
        for j=1:numel(TrackList)
            track = TrackList{j};
            % Convert track trajectory to LLA and plot it
            means = [track.Trajectory.Mean];
            if PlotLatLon
                [lat,lon,~] = ned2geodetic(means(3,:),...
                                           means(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(track.State.Mean(4,end),...
                                                   track.State.Mean(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);
                if ShowTrackInfo
                    speed_kmph = sqrt(track.State.Mean(4,end)^2+track.State.Mean(2,end)^2)*3.6;
                    speed_knot = speed_kmph/1.852;
                    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(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','w');
                end
            else
                traj_length = size(lon,2);
                if(traj_length>NumPersistFrames)
                    start = traj_length-NumPersistFrames;
                else
                    start = 1;
                end
                plots(end+1) = plot(ax(1), means(1,start:end),means(3,start:end),'-.k','LineWidth',2);
                plots(end+1) = plot(ax(1), means(1,end),means(3,end),'ks','MarkerSize',15);

                x_vel = track.State.Mean(3,end)*cos(track.State.Mean(4,end));
                y_vel = track.State.Mean(3,end)*sin(track.State.Mean(4,end));
                plots(end+1) = quiver(ax(1), means(1,end),means(2,end),20*x_vel,20*y_vel,'r','LineWidth',1.5);

                if ShowTrackInfo
                    speed_kmph = track.State.Mean(3,end)*3.6;
                    speed_knot = speed_kmph/1.852;
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)+50,strcat("Sog:",num2str(speed_knot,2)," kt"),'FontSize',8,'Color','k');
                    plots(end+1) = text(ax(1), means(1,end)+60,means(2,end)-50,strcat("PoE:",num2str(track.ExistenceProbability*100,3)," %"),'FontSize',8,'Color','k');
                end
            end
        end

        % Add axis labels
%         xlabel('Longitude')
%         ylabel('Latitude')
        title(datestr(timestamp_k));
%         set(gca,'visible','off')
        pause(0.0001)
        % Store video frame
        if(Record)
            F(k) = getframe(ax(1));
        end
    end

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 value = r8_sin ( x )

%*****************************************************************************80
%
%% R8_SIN evaluates the sine of an R8 argument.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    27 September 2011
%
%  Author:
%
%    Original FORTRAN77 version by Wayne Fullerton.
%    MATLAB version by John Burkardt.
%
%  Reference:
%
%    Wayne Fullerton,
%    Portable Special Function Routines,
%    in Portability of Numerical Software,
%    edited by Wayne Cowell,
%    Lecture Notes in Computer Science, Volume 57,
%    Springer 1977,
%    ISBN: 978-3-540-08446-4,
%    LC: QA297.W65.
%
%  Parameters:
%
%    Input, real X, the argument.
%
%    Output, real VALUE, the sine of X.
%
  persistent ntsn
  persistent pi2rec
  persistent pihi
  persistent pilo
  persistent pirec
  persistent sincs
  persistent xmax
  persistent xsml
  persistent xwarn

  pi2rec = 0.63661977236758134307553505349006;
  pihi = 3.140625;
  pilo = 9.6765358979323846264338327950288E-04;
  pirec = 0.31830988618379067153776752674503;

  if ( isempty ( ntsn ) )

    sincs = [ ...
     -0.374991154955873175839919279977323464, ...
     -0.181603155237250201863830316158004754, ...
     +0.005804709274598633559427341722857921, ...
     -0.000086954311779340757113212316353178, ...
     +0.000000754370148088851481006839927030, ...
     -0.000000004267129665055961107126829906, ...
     +0.000000000016980422945488168181824792, ...
     -0.000000000000050120578889961870929524, ...
     +0.000000000000000114101026680010675628, ...
     -0.000000000000000000206437504424783134, ...
     +0.000000000000000000000303969595918706, ...
     -0.000000000000000000000000371357734157, ...
     +0.000000000000000000000000000382486123, ...
     -0.000000000000000000000000000000336623, ...
     +0.000000000000000000000000000000000256 ]';

    ntsn = r8_inits ( sincs, 15, 0.1 * r8_mach ( 3 ) );
    xsml = sqrt ( 2.0 * r8_mach ( 3 ) );
    xmax = 1.0 / r8_mach( 4 );
    xwarn = sqrt ( xmax );

  end

  y = abs ( x );

  if ( xmax < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_SIN - Warning!\n' );
    fprintf ( 1, '  No precision because |X| is big.\n' );
    value = 0.0;
    return
  end

  if ( xwarn < y )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'R8_SIN - Warning!\n' );
    fprintf ( 1, '  Answer < half precision because |X| is big.\n' );
  end

  value = x;
  if ( y < xsml )
    return
  end

  xn = r8_aint ( y * pirec + 0.5 );
  n2 = r8_aint ( mod ( xn, 2.0 ) + 0.5 );

  sgn = x;
  if ( n2 ~= 0 )
    sgn = - sgn;
  end

  f = ( y - xn * pihi ) - xn * pilo;

  xn = 2.0 * ( f * pi2rec ) * ( f * pi2rec ) - 1.0;

  value = f + f * r8_csevl ( xn, sincs, ntsn );

  if ( sgn < 0.0 )
    value = - value;
  end

  if ( value < - 1.0 )
    value = - 1.0;
  elseif ( 1.0 < value )
    value = + 1.0;
  end

  return
end


text:
clear;
clc
close all;

addpath('../../library/nav_lib');

%% attidude

true_attitude = deg2rad([4 -6 30])';
roll_meas_err = deg2rad(0.1);
pitch_meas_err = deg2rad(-0.15);

%% Earth's Magnetic Field
flux = 45;
declination= deg2rad(10);
inclination = deg2rad(40);

delta_alpha_nE= deg2rad(0.3);

%% Other Souces of Magnetism
MAn = [2 -1.8 3]'; % local magnetic anomalies(N,E,D)
hard = [-1 -0.5 1]'; % hard-iron magnetism(body-frame axes)
soft = [0.012 -0.014 0.006; -0.008 0.007 0.017; 0.003 -0.019 -0.011]; % soft-iron magnetism, body-rame-axes

%%
Cn2b = ch_eul2m(true_attitude);

%% flux desity of Earth's magnetic filed
m_En = [cos(declination)*cos(inclination)  sin(declination)*cos(inclination) sin(inclination)]'*flux;

%% Total magnetic flux density
m_mb = hard + (eye(3) + soft)*Cn2b*(m_En + MAn);

%% Calculate Magnetic Heading Measurement
roll_meas = true_attitude(1) + roll_meas_err;
pitch_meas =  true_attitude(2) +pitch_meas_err;

yaw = atan2(-m_mb(2)*cos(roll_meas) + m_mb(3)*sin(roll_meas), m_mb(1)*cos(pitch_meas) + m_mb(2)*sin(pitch_meas)*sin(roll_meas) + m_mb(3)*cos(roll_meas)*sin(pitch_meas));

%% Calculate True Heading
data_base_magnetic_declination = declination + delta_alpha_nE

%% True heading
true_heading = yaw+ data_base_magnetic_declination;
%rad2deg(true_heading)

%% Heading error
error = true_heading - true_attitude(3);
rad2deg(error)

author: yandld


text:
function [colChar]=xlsColNum2Str(colNum)
%XLSCOLNUM2STR takes in an array of numbers and returns a cellular array
%of the same size with cell of corresponding Excel column names.
%
%For example:
%n=[1  10;
%   53 256]
%c=xlsColNum2Str(n);
%c={'A' , 'J';
%   'BA', 'IV'}
%Note: up to Excel 2003 the number of columns was limited to 256, as of
%Excel 2007 the number of columns has increased to 16,384 or 'XFD'
%This function is designed to take accept any integer so proper handling
%of the number of columns should be taken care of outside this function

    colChar=cell(size(colNum)); %blank cell array

    % find max number of characters (AA n=2)
    numOfChars=ceil(max(colNum)/26)-1;
    n=1;

    while numOfChars>=1
        numOfChars=ceil(numOfChars/26)-1;
        n=n+1;
    end

    remainder=num2cell(colNum);

    for s=n:-1:1
        if s>1
            %find limits
            % if n=2 then the columns go from AA to ZZ or 27 to 702
            L=sum(26.^(1:s-1))+1; % lower limit
            U=sum(26.^(1:s));   %upper limit
            %place current character to right of previous
            colChar(colNum>=L & colNum<=U) = ...
                cellfun(@(x,y) ([x char(ceil((y-(L-1))/26^(s-1))+64)]),...
                colChar(colNum>=L & colNum<=U),...      % x
                remainder(colNum>=L & colNum<=U),...    % y
                'UniformOutput',false);
            %calculate the remaining string
            %for example if last string was 'ABA' the 'A' was placed to the
            %right of the previous string and now 'BA' is remaining
            remainder(colNum>=L & colNum<=U)=...
                cellfun(@(x,y) (y-26^(s-1)*(double(x(end))-64)),...
                colChar(colNum>=L & colNum<=U),...
                remainder(colNum>=L & colNum<=U),'UniformOutput',false);
            colNum=cell2mat(remainder);
        else
             colChar=cellfun(@(x,y) ([x char(y+64)]),...
                 colChar,remainder,'UniformOutput',false);
        end
    end
end

author: cerr

repo: CERR

save_path: github-repos/MATLAB/cerr-CERR


text:
function [planar] = constructplanargrad(cfg, grad)

% CONSTRUCTPLANARGRAD constructs a planar gradiometer array from an axial gradiometer
% definition. This can be used to compute the planar field gradient for a known
% (estimated) source configuration.
%
% Use as
%   [grad_planar] = constructplanargrad(cfg, grad_axial)
%
% Where cfg contains the following configuration details
%   cfg.baseline_axial   = number (default is 5)
%   cfg.baseline_planar  = number (default is 0.5)
%   cfg.planaraxial      = 'no' or 'yes' (default)
%
% The option planaraxial='yes' specifies that the planar gradiometers
% should consist of axial gradiometers, to make them comparable with
% Ole Jensens planar gradient computation. If planaraxial='no', the
% planar gradiometers will be more or less similar to the Neuromag
% system.
%
% The input grad can be a CTF type axial gradiometer definition, but
% just as well be a magnetometer definition. This function only assumes
% that
%   grad.coilpos
%   grad.coilori
%   grad.label
% exist and that the first Nlabel channels in pnt and ori should be
% used to compute the position of the coils in the planar gradiometer
% channels.

% Copyright (C) 2004, Robert Oostenveld
%
% This file is part of FieldTrip, see http://www.fieldtriptoolbox.org
% for the documentation and details.
%
%    FieldTrip 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.
%
%    FieldTrip 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 FieldTrip. If not, see <http://www.gnu.org/licenses/>.
%
% $Id$

if ~isfield(cfg, 'planaraxial'),     cfg.planaraxial = 'yes';   end
if ~isfield(cfg, 'baseline_axial'),  cfg.baseline_axial  = 5;   end
if ~isfield(cfg, 'baseline_planar'), cfg.baseline_planar = 0.5; end

Nchan = length(grad.label);

% these will hold all the coil positions
lo_posx = zeros(Nchan,3);
lo_negx = zeros(Nchan,3);
lo_posy = zeros(Nchan,3);
lo_negy = zeros(Nchan,3);
hi_posx = zeros(Nchan,3);
hi_negx = zeros(Nchan,3);
hi_posy = zeros(Nchan,3);
hi_negy = zeros(Nchan,3);

for chan=1:Nchan
  % Attach a local coordinate system to this gradiometer:
  %   the origin at the location of its bottom coil
  %   the z-axis pointing outwards from the head
  %   the x-axis pointing horizontal w.r.t. the head
  %   the y-axis pointing vertical, i.e. approximately towards the vertex
  this_o = grad.chanpos(chan,:);
  this_z = grad.chanori(chan,:);
  this_z = this_z / norm(this_z);
  this_x = cross([0 0 1], this_z);
  if all(this_x==0)
    this_x = [1 0 0];
  else
    this_x = this_x / norm(this_x);
  end
  this_y = cross(this_z, this_x);

  % compute the position of all the 8 coils per channel
  lo_posx(chan,:) = this_o + (cfg.baseline_planar/2) * this_x;
  lo_negx(chan,:) = this_o - (cfg.baseline_planar/2) * this_x;
  lo_posy(chan,:) = this_o + (cfg.baseline_planar/2) * this_y;
  lo_negy(chan,:) = this_o - (cfg.baseline_planar/2) * this_y;
  hi_posx(chan,:) = lo_posx(chan,:) + cfg.baseline_axial * this_z;
  hi_negx(chan,:) = lo_negx(chan,:) + cfg.baseline_axial * this_z;
  hi_posy(chan,:) = lo_posy(chan,:) + cfg.baseline_axial * this_z;
  hi_negy(chan,:) = lo_negy(chan,:) + cfg.baseline_axial * this_z;
end

% start with an empty planar gradiometer definition
planar = [];

if strcmp(cfg.planaraxial, 'yes')
  % combine all the 8 coils into a single sensor
  planar.coilpos = [
    lo_posx
    lo_negx
    lo_posy
    lo_negy
    hi_posx
    hi_negx
    hi_posy
    hi_negy
  ];

  % the orientation of all the coils of a single sensor should be the same
  planar.coilori = [
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
  ];

  e = eye(Nchan);
  z = zeros(Nchan);

  % the linear combination matrix should be 2*Nchan x 8*Nchan
  planar.tra = [
    e -e  z  z -e  e  z  z    % this is for the horizontal gradients
    z  z  e -e  z  z -e  e    % this is for the vertical gradients
  ];

else
  % combine only the 4 lower coils into a single sensor
  planar.coilpos = [
    posx
    negx
    posy
    negy
  ];

  % the orientation of all the coils of a single gradiometer should be the same
  planar.coilori = [
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
    grad.coilori(1:Nchan,:)
  ];

  e = eye(Nchan);
  z = zeros(Nchan);

  % the linear combination matrix should be 2*Nchan x 4*Nchan
  planar.tra = [
    e -e  z  z    % this is for the horizontal gradients
    z  z  e -e    % this is for the vertical gradients
  ];

end

for chan=1:Nchan
  planar.label{chan      } = [grad.label{chan} '_dH'];
  planar.label{chan+Nchan} = [grad.label{chan} '_dV'];
end

planar.label = planar.label(:);
planar.tra   = planar.tra / cfg.baseline_planar;
planar.chanpos = [grad.chanpos; grad.chanpos];
planar.chanori = [grad.chanori; grad.chanori];

try
  planar.unit  = grad.unit;
end

% add information about the version of this function to the configuration
cfg.version.name = mfilename('fullpath');
cfg.version.id   = '$Id$';

% rememember the exact configuration details in the output
planar.cfg = cfg;

author: fieldtrip

repo: fieldtrip

save_path: github-repos/MATLAB/fieldtrip-fieldtrip


text:
function input2conv(in)
    global config mem;
    for m = 1:config.batch_size
        for n = 1:config.chs
            mem.layer_inputs{1}((n-1)*config.kernel_size(1, 1)*config.kernel_size(1, 2)+1:n*config.kernel_size(1, 1)*config.kernel_size(1, 2), (m-1)*(size(mem.layer_inputs{1}, 2)/config.batch_size)+1:m*size(mem.layer_inputs{1}, 2)/config.batch_size) = ...
                    config.SCALE_INPUT(config.IM2COL(in(:,:,n,m), [config.kernel_size(1, 1), config.kernel_size(1, 2)]));
        end
    end
end

author: jimmy-ren


text:
function sel = filtertable( table, filter )

fnames = fieldnames(filter);
sel = true(size(table, 1), 1);

for fi = 1:numel(fnames)
  fname = fnames{fi};
  fvalue = filter.(fname);
  sel = sel & ismember(table.(fname), fvalue);
end

end

author: hpatches

repo: hpatches-benchmark


text:
function y = tapas_softmax_mu3_sim(r, infStates, p)
% Simulates observations from a Boltzmann distribution with volatility as temperature
%
% --------------------------------------------------------------------------------------------------
% Copyright (C) 2017-2019 Christoph Mathys, TNU, UZH & ETHZ
%
% This file is part of the HGF toolbox, which is released under the terms of the GNU General Public
% Licence (GPL), version 3. You can redistribute it and/or modify it under the terms of the GPL
% (either version 3 or, at your option, any later version). For further details, see the file
% COPYING or <http://www.gnu.org/licenses/>.

% Predictions or posteriors?
pop = 1; % Default: predictions
if r.c_obs.predorpost == 2
    pop = 3; % Alternative: posteriors
end

% Assumed structure of infStates:
% dim 1: time (ie, input sequence number)
% dim 2: HGF level
% dim 3: choice number
% dim 4: 1: muhat, 2: sahat, 3: mu, 4: sa

% Number of choices
nc = size(infStates,3);

% Belief trajectories at 1st level
states = squeeze(infStates(:,1,:,pop));

% Log-volatility trajectory
mu3 = squeeze(infStates(:,3,1,3));

% Inverse decision temperature
be = exp(-mu3);
be = repmat(be,1,nc);

% Partition functions
Z = sum(exp(be.*states),2);
Z = repmat(Z,1,nc);

% Softmax probabilities
prob = exp(be.*states)./Z;

% Initialize random number generator
if isnan(r.c_sim.seed)
    rng('shuffle');
else
    rng(r.c_sim.seed);
end

% Draw responses
n = size(infStates,1);
y = NaN(n,1);

for j=1:n
    y(j) = find(mnrnd(1, prob(j,:)));
end

end

author: translationalneuromodeling

repo: tapas


text:
function cvx_optval = norm( x, p )

%   Disciplined convex programming information:
%       NORM is convex, except when P<1, so an error will result if
%       these non-convex "norms" are used within CVX expressions. NORM
%       is nonmonotonic, so its input must be affine.

%
% Argument map
%

persistent remap1 remap2
if isempty( remap2 ),
    remap1 = cvx_remap( 'log-convex' );
    remap2 = cvx_remap( 'affine', 'log-convex' );
end

%
% Check arguments
%

error( nargchk( 1, 2, nargin ) );
if nargin < 2,
    p = 2;
elseif ~isequal( p, 'fro' ) && ( ~isnumeric( p ) || ~isreal( p ) || p < 1 ),
    error( 'Second argument must be a real number between 1 and Inf, or ''fro''.' );
end
if ndims( x ) > 2,
    error( 'norm is not defined for N-D arrays.' );
end

[m,n] = size(x);
if m == 1 || n == 1 || isequal( p, 'fro' ),

    %
    % Vector norms
    %

    if isempty( x ),
        cvx_optval = cvx( 0 );
        return
    end
    if isequal( p, 'fro' ),
        p = 2;
    end
    x = svec( x, p );
    if length( x ) == 1,
        p = 1;
    end
    xc = cvx_classify( x );
    if ~all( remap2( xc ) ),
        error( 'Disciplined convex programming error:\n    Cannot perform the operation norm( {%s}, %g )', cvx_class( x ), p );
    end
    switch p,
        case 1,
            cvx_optval = sum( abs( x ) );
        case Inf,
            cvx_optval = max( abs( x ) );
        otherwise,
            tt = remap1( xc );
            if all( tt ),
                cvx_optval = ( sum( x .^ p ) ) .^ (1/p);
            else
                if nnz( tt ) > 1,
                    tt = tt ~= 0;
                    xx = cvx_subsref( x, tt );
                    xx = ( sum( xx .^ p ) ) .^ (1/p);
                    x  = [ cvx_subsref( x, ~tt ) ; cvx_accept_convex( xx ) ];
                end
                n = length( x );
                if p == 2,
                    cvx_begin
                        epigraph variable z
                        { x, z } == lorentz( n, [], ~isreal( x ) );
                    cvx_end
                else
                    if isreal( x ),
                        cmode = 'abs';
                    else
                        cmode = 'cabs';
                    end
                    cvx_begin
                        epigraph variable z
                        variable y( n )
                        { [ y, z*ones(n,1) ], x } == geo_mean_cone( [n,2], 2, [1/p,1-1/p], cmode );
                        sum( y ) == z;
                    cvx_end
                end
            end
    end

else

    %
    % Matrix norms
    %

    if ~cvx_isaffine( x ),
        error( 'Disciplined convex programming error:\n    Cannot perform the operation norm( {%s}, %g )\n   when the first argument is a matrix.', cvx_class( xt ), p );
    end
    switch p,
        case 1,
            cvx_optval = max( sum( abs( x ), 1 ), [], 2 );
        case Inf,
            cvx_optval = max( sum( abs( x ), 2 ), [], 1 );
        case 2,
            cvx_optval = sigma_max( x );
        otherwise,
            error( 'The only matrix norms available are 1, 2, Inf, and ''fro''.' );
    end

end

% Copyright 2010 Michael C. Grant and Stephen P. Boyd.
% See the file COPYING.txt for full copyright information.
% The command 'cvx_where' will show where this file is located.

repo: REST


text:
function sudoku_test06 ( )

%*****************************************************************************80
%
%% SUDOKU_TEST06 demonstrates SUDOKU_ROTATE.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    18 February 2006
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'SUDOKU_TEST06:\n' );
  fprintf ( 1, '  SUDOKU_ROTATE rotates a Sudoku by 90 degrees;\n' );

  s = [ ...
       3,4,1,9,2,7,5,6,8; ...
       6,9,2,1,8,5,7,3,4; ...
       8,5,7,4,6,3,1,9,2; ...
       1,3,4,2,9,6,8,7,5; ...
       2,7,8,5,3,4,6,1,9; ...
       5,6,9,7,1,8,4,2,3; ...
       4,2,5,3,7,1,9,8,6; ...
       9,1,6,8,4,2,3,5,7; ...
       7,8,3,6,5,9,2,4,1];

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the initial Sudoku puzzle:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  s = sudoku_rotate ( s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the rotated Sudoku:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  s = sudoku_rotate ( s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the rotated Sudoku:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  s = sudoku_rotate ( s );

  fprintf ( 1, '\n' );
  fprintf ( 1, '  Here is the rotated Sudoku:\n' );
  fprintf ( 1, '\n' );

  sudoku_print ( s );

  return
end


text:
%% runGmsh
% Below is a demonstration of the features of the |runGmsh| function

%%
clear; close all; clc;

%% Syntax
% |;|

%% Description
% UNDOCUMENTED
%% Examples
%
%%
%
% <<gibbVerySmall.gif>>
%
% _*GIBBON*_
% <www.gibboncode.org>
%
% _Kevin Mattheus Moerman_, <gibbon.toolbox@gmail.com>
%%
% _*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/>.


text:
function r8_walsh_1d_test ( )

%*****************************************************************************80
%
%% R8_WALSH_1D_TEST tests R8_WALSH_1D;
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    19 April 2009
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'R8_WALSH_1D_TEST\n' );
  fprintf ( 1, '  R8_WALSH_1D evaluates 1D Walsh functions:\n' );
  fprintf ( 1, '\n' );
  fprintf ( 1, '  X  W(+2) W(+1) W(0) W(-1) W(-2) W(-3)\n' );
  fprintf ( 1, '\n' );

  for i = 0 : 32

    x = i / 4.0;

    wp2 = r8_walsh_1d ( x,  2 );
    wp1 = r8_walsh_1d ( x,  1 );
    w0  = r8_walsh_1d ( x,  0 );
    wm1 = r8_walsh_1d ( x, -1 );
    wm2 = r8_walsh_1d ( x, -2 );
    wm3 = r8_walsh_1d ( x, -3 );

    fprintf ( 1, '  %10f  %4f  %4f  %4f  %4f  %4f  %4f\n', ...
      x, wp2, wp1, w0, wm1, wm2, wm3 );

  end

  return
end


text:
% counting.m
%
%
%
% author: Young-Seok Kweon
% created: 2020.11.12
%% init
clc; clear; close all;
%% load

list=dir('Dataset\*.mat');

%% counting

name=[];
for i=1:length(list)
    temp=list(i).name;
    temp=split(temp,'_');
    name(i)=str2num(temp{1});
end
%%
count(max(name))=0;
for i=1:length(list)
    count(name(i)) = count(name(i))+1;
end

for i=unique(count)
    if i==0
        continue;
    end

    n(i)=sum(count==i);
end

%%
dir_='F:\wsc\datasets\wsc-dataset-0.1.0.csv';

[num, txt, raw]=xlsread(dir_);

year=count;
c=count;
year(:)=0;
c(:)=0;
for i=1:length(list)
    if c(num(i,1))~=0
        year(num(i,1))=year(num(i,1))+num(i,3)-c(num(i,1));
        c(num(i,1))=num(i,3);
    else
        c(num(i,1))=num(i,3);
    end
end
year=year./count;
for i=unique(count)
    if i==0
        continue;
    end

    n(i)=mean(year(count==i));
end

author: PatternRecognition

repo: OpenBMI


text:
function [costunit, selectedf] = choosefactory(factorycap,factholdcost,transcost,setupcost,factory,trkcap,salespoint,selectedsp);

    costPerUnitValues = inf(1,length(factorycap));

    for i = 1:length(factorycap)
        if (factory{i}(1,2)>=salespoint{selectedsp}(1,2))
            costPerUnitValues(1,i) = transcost(i,selectedsp) / trkcap + ( ( factorycap(i) - salespoint{selectedsp}(1,2) ) / factorycap(i) ) * factholdcost(i);
        elseif(factory{i}(1,1)==0 && (factorycap(i)+factory{i}(1,2))>=salespoint{selectedsp}(1,2))
            costPerUnitValues(1,i) = (setupcost(i) * (1 - factory{i}(1,1))) / salespoint{selectedsp}(1,2) + transcost(i,selectedsp) / trkcap + ( ( factorycap(i) - salespoint{selectedsp}(1,2) ) / factorycap(i) ) * factholdcost(i);
        end

    end;

    [costunit,selectedf] = min(costPerUnitValues);


text:
function a = rutis2_inverse ( )

%*****************************************************************************80
%
%% RUTIS2_INVERSE returns the inverse of the RUTIS2 matrix.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    20 October 2007
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Joan Westlake,
%    A Handbook of Numerical Matrix Inversion and Solution of
%    Linear Equations,
%    John Wiley, 1968,
%    ISBN13: 978-0471936756,
%    LC: QA263.W47.
%
%  Parameters:
%
%    Output, real A(4,4), the matrix.
%

%
%  Note that the matrix entries are listed by row.
%
  a(1:4,1:4) = [ ...
     28.0,  -22.0, -1.0, -1.0; ...
    -22.0,   28.0, -1.0, -1.0; ...
     -1.0,   -1.0, 17.0, -8.0; ...
     -1.0,   -1.0, -8.0, 17.0 ];

  a(1:4,1:4) = a(1:4,1:4) / 50.0;

  return
end


text:
function Mz = ir_equations(params, seqFlag, approxFlag)
%IR_EQUATIONS Analytical equations for the longitudinal magnetization of
%steady-state inversion recovery experiments with either a gradient echo
%(GRE-IR) or spin-echo (SE-IR) readouts.
%   Reference: Barral, J. K., Gudmundson, E. , Stikov, N. , Etezadi?Amoli,
%   M. , Stoica, P. and Nishimura, D. G. (2010), A robust methodology for
%   in vivo T1 mapping. Magn. Reson. Med., 64: 1057-1067.
%   doi:10.1002/mrm.22497
%
%   params: struct with the required parameters for the sequence and
%   approximation. See below for list.
%
%   seqFlag: String. Either 'GRE-IR' or 'SE-IR'
%   approxFlag: Integer between 1 and 4.
%       1: General equation (no approximation).
%       2: Ideal 180 degree pulse approximation of case 1.
%       3: Ideal 90 degree pulse approximation of case 2, and readout term
%          absorbed into constant.
%       4: Long TR (TR >> T1) approximation of case 3.
%
%   **PARAMS PROPERTIES**
%   All times in seconds, all angles in degrees.
%  'GRE-IR'
%       case 1: T1, TR, TI, EXC_FA, INV_FA, constant (optional)
%       case 2: T1, TR, TI, EXC_FA, constant (optional)
%       case 3: T1, TR, TI, constant (optional)
%       case 4: T1, TI, constant (optional)
%
%  'SE-IR'
%       case 1: Same as 'GRE-IR' case + SE_FA, TE
%       case 2: Same as 'GRE-IR' case + TE
%       case 3: Same as 'GRE-IR' case + TE
%       case 4: Same as 'GRE-IR' case
%

switch seqFlag
    case 'GRE-IR'
        switch approxFlag
            case 1 % General
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    EXC_FA = params.EXC_FA; % Excitation pulse in deg
                    INV_FA = params.INV_FA; % Inversion pulse in deg

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* ( (1-cosd(INV_FA).*exp(-TR/T1) - (1-cosd(INV_FA)).*exp(-TI./T1)) ./ (1-cosd(INV_FA).*cosd(EXC_FA).*exp(-TR./T1)) );
                catch
                    error('ir_equations.GRE-IR.case1: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
            case 2 % Ideal 180 pulse
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    EXC_FA = params.EXC_FA; % in deg

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* (1 - 2*exp(-TI./T1) + exp(-TR./T1));
                catch
                    error('ir_equations.GRE-IR.case2: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
            case 3 % Ideal 90 pulse.
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* (1 - 2*exp(-TI./T1) + exp(-TR./T1));
                catch
                    error('ir_equations.GRE-IR.case3: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end

            case 4 % Long TR (TR >> T1)
                try
                    T1 = params.T1;
                    TI = params.TI;

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* (1 - 2*exp(-TI./T1));
                catch
                    error('ir_equations.GRE-IR.case4: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end

            otherwise
                error('ir_equations: Unknown flag. Run `help ir_equations` for more info.')
        end
    case 'SE-IR'
        switch approxFlag
            case 1 % General equation
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    TE = params.TE;

                    EXC_FA = params.EXC_FA; % Excitation pulse in deg
                    INV_FA = params.INV_FA; % Inversion pulse in deg
                    SE_FA = params.SE_FA; % Inversion pulse in deg

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* ( (1-cosd(INV_FA).*cosd(SE_FA).*exp(-TR/T1) - cosd(INV_FA).*(1-cosd(SE_FA)).*exp(-(TR-(TE/2))./T1) - (1-cosd(INV_FA)).*exp(-TI./T1)) ./ (1-cosd(INV_FA).*cosd(EXC_FA).*cosd(SE_FA).*exp(-TR./T1)) );
                catch
                    error('ir_equations.SE-IR.case1: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
            case 2 % Ideal 180 pulses
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    TE = params.TE;

                    EXC_FA = params.EXC_FA; % Excitation pulse in deg

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* ( (1-exp(-TR/T1) + 2.*exp(-(TR-(TE/2))./T1) - 2.*exp(-TI./T1)) ./ (1-cosd(EXC_FA).*exp(-TR./T1)) );
                catch
                    error('ir_equations.SE-IR.case2: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
            case 3 % Ideal 90 pulse
                try
                    T1 = params.T1;
                    TR = params.TR;
                    TI = params.TI;

                    TE = params.TE;

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* (1-exp(-TR/T1) + 2.*exp(-(TR-(TE/2))./T1) - 2.*exp(-TI./T1));
                catch
                    error('ir_equations.SE-IR.case3: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
            case 4 % Long TR
                try
                    T1 = params.T1;
                    TI = params.TI;

                    try
                        constant = params.constant;
                    catch
                        constant = 1;
                    end

                    Mz= constant .* (1 - 2.*exp(-TI./T1));
                catch
                    error('ir_equations.SE-IR.case4: Incorrect parameters for given flag.  Run `help ir_equations` for more info.')
                end
        end
    otherwise
        error('ir_equations.seqFlag: Incorrect seqFlag arguement. Must be either GRE-IR or SE-IR.')
end

author: qMRLab

repo: qMRLab

save_path: github-repos/MATLAB/qMRLab-qMRLab


text:
function output=user5sym()

%**************************************************************************
 %*************************************************************************
                %//input signal form transmitter side for symbol based\\
 clear all;
close all;
clc;
% us=input('enter the number of users');
% nb1=input('number of bits for user1');
% m1=input('enter the message sequence1');
% nb2=input('number of bits for user2');
% m2=input('enter the message sequence2');
% nb3=input('number of bits for user3');
% m3=input('enter the message sequence3');
% nb4=input('number of bits for user4');
% m4=input('enter the message sequence4');
% nb5=input('number of bits for user5');
% m5=input('enter the message sequence5');
us=3;
% nb1=5;
% nb2=5;
% nb3=5;
% nb4=5;
% nb5=5;
% m1=[1 0 0 0 0];
% m2=[1 0 0 1 0];
% m3=[1 0 1 1 0];
% m4=[1 0 1 0 1 ];
% m5=[1 0 1 0 0];

nb1=10;
nb2=10;
nb3=10;
nb4=10;
nb5=10;
m1=[1 0 0 0 0 1 0 1 0 0];
m2=[1 0 0 1 0 0 1 0 0 0];
m3=[1 0 1 0 0 1 0 1 0 1];
m4=[1 0 1 0 1 0 0 0 1 0];
m5=[1 0 1 0 0 0 1 0 0 0];

                    %//message1 into polar form\\

for i=1:nb1
    if m1(i)==1
        mb1(i)=m1(i);
    else
        mb1(i)=-1;
    end
end
  % disp(mb1);

                %//message2 into polar form\\

for i=1:nb2
    if m2(i)==1
        mb2(i)=m2(i);
    else
        mb2(i)=-1;
    end
end
  %disp(mb2);

                %//message3 into polar form\\

for i=1:nb3
    if m3(i)==1
        mb3(i)=m3(i);
    else
        mb3(i)=-1;
    end
end
  %disp(mb3);

       %//message4 into polar form\\

for i=1:nb4
    if m4(i)==1
        mb4(i)=m4(i);
    else
        mb4(i)=-1;
    end
end
  %disp(mb4);

       %//message5 into polar form\\

for i=1:nb5
    if m5(i)==1
        mb5(i)=m5(i);
    else
        mb5(i)=-1;
    end
end
  %disp(mb5);

                  %//generation of maximal length sequence1\\

  f1=1;f2=1;f3=0;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f2,f5);
    f5=f4;
    op1(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op1);

                    %//msequence1 into polar form\\

    for i=1:31

         if op1(i)==1
             pb1(i)=1;
         else
         pb1(i)=-1;
        end
    end
    %disp(op1);
    %disp(pb1);

                    %//generation of maximal length sequence2\\

  f1=1;f2=0;f3=0;f4=1;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f3,f5);
    f5=f4;
    op2(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op2);

                    %//msequence2 into polar form\\

    for i=1:31
        if op2(i)==1
             pb2(i)=1;
         else
         pb2(i)=-1;
        end
    end
    %disp(op2);
    %disp(pb2);

                     %//generation of maximal length sequence3\\

  f1=1;f2=1;f3=1;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op3(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op3);

                    %//msequence3 into polar form\\

    for i=1:31

         if op3(i)==1
             pb3(i)=1;
         else
         pb3(i)=-1;
        end
    end
    %disp(op3);
    %disp(pb3);

                     %//generation of maximal length sequence4\\

  f1=0;f2=1;f3=1;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op4(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op4);

                    %//msequence4 into polar form\\

    for i=1:31

         if op4(i)==1
             pb4(i)=1;
         else
         pb4(i)=-1;
        end
    end
    %disp(op4);
    %disp(pb4);

                     %//generation of maximal length sequence5\\

  f1=1;f2=1;f3=0;f4=0;f5=1;
m=(2^5)-1;
for j=1:m
    p=xor(f4,f5);
    f5=f4;
    op5(j)=f5;
    f4=f3;
    f3=f2;
    f2=f1;
    f1=p;

end
  %disp(op5);

                    %//msequence5 into polar form\\

    for i=1:31

         if op5(i)==1
             pb5(i)=1;
         else
         pb5(i)=-1;
        end
    end
    %disp(op5);
    %disp(pb5);

                    %//spreading signals from transmitter side\\

                         %//first transmit bit\\

k=1;
for i=1:nb1
    for j=1:31
        tb1(k)=mb1(i)*pb1(j);
        k=k+1;
    end
end
%disp(tb1);

                     %//second transmit bit\\

k=1;
for i=1:nb2
    for j=1:31
        tb2(k)=mb2(i)*pb2(j);
        k=k+1;
    end
end
%disp(tb2);

                    %//third transmit bit\\

k=1;
for i=1:nb3
    for j=1:31
        tb3(k)=mb3(i)*pb3(j);
        k=k+1;
    end
end
%disp(tb3);

                %//fourth transmit bit\\

k=1;
for i=1:nb4
    for j=1:31
        tb4(k)=mb4(i)*pb4(j);
        k=k+1;
    end
end
%disp(tb4);

             %//fifth transmit bit\\

k=1;
for i=1:nb5
    for j=1:31
        tb5(k)=mb5(i)*pb5(j);
        k=k+1;
    end
end
%disp(tb5);

                 %//Addition of five  signals transmitted in the channel\\

n=1;
for i=1:k-1
    tb(n)=tb1(n)+tb2(n)+tb3(n)+tb4(n)+tb5(n);
    n=n+1;
end
%disp(tb);
chn=awgn(tb,10);

                         %//%receiver side\\
                    %//reconstruction of user2 signal\\

s2=0;
t2=1;
for i=0:31:k-32
    for j=1:31
        ob2(t2)=chn(i+j)*pb2(j)+s2;
        s2=ob2(t2);
        if ob2(t2)>0
            ob2(t2)=1;
        else ob2(t2)=0;
        end
    end
    t2=t2+1;
    s2=0;
end
%disp('second message is')
%disp(ob2);
n1=awgn(ob2,10);
 figure(1);
grid on;
subplot(2,2,1);
plot(ob2);
title('Received signal with out noise for SB');
xlabel('Time');
ylabel('Amplitude');
%g=awgn(ob2,10);
%disp(g);
subplot(2,2,2);
plot(chn);
title('Received signal with noise for SB');
xlabel('Time');
ylabel('Amplitude');

                    %//reconstruction of received signal from noise\\

for i=1:nb2
     if chn(i)>0
            g2(i)=1;
        else g2(i)=0;
        end
    end

% disp(g2);

                 % //calculating optimum weight using minimum variance\\

 ob=randint(1,31);
q=length(10);
%q1=input('enter the sequence');
s=randsrc(1,q);
c1=complex(ob,s);
%disp(c);
d=conj(c1);
d1=d*d';
%disp(d1);
m=mean(d1);
lamda=5;
teta=45;
k=((2*pi)/lamda);
K=4;
for i=0:K-1
    v=(exp(j*k*i*d*sin(teta)))';
end
%disp(v);
v1=conj(v)';
teta=30;
k=((2*pi)/lamda);
K=5;
for i=0:K-1
    eta=exp(sqrt(-1)*k*i*d*sin(teta))';
end
%u1=[2+2j 3-1j 4+3j 5+2j]'
%m=randsrc(1,124);
%a=randsrc(1,124);
l=length(10);
%l1=input('enter the sequence');
u1=complex(ob2,l);
u=eta*u1;
u2=conj(u)';
u3=u*u2;
Ru=mean(u3);

%Ru1=imresize(Ru,[124 124]);
%a1=inv(Ru1);
%a1=Ru1';
a1=Ru';
%disp(a1);
%calculating beta value
%g=1;
dem1=v*v1*a1;
beta1=dem1.^-1;
 %optimum weight
%Wopt=beta*a1*v';
temp1=beta1*v';
Wopt1=temp1*a1;
%disp(Wopt1);
 subplot(2,2,3);
plot(Wopt1);
title('Received signal with optimum weight for SB');
xlabel('Time');
ylabel('weight');

                    %//beamformer output\\

o1=conj(Wopt1)*ob2;
%disp('Beamformer output for SB')
%disp(o1);
 subplot(2,2,4);
plot(o1);
title('Beamformer output  for SB');
xlabel('Time');
ylabel('Beam former output');

                    %//SINR CALCULATION FOR SB CONFIGURATION\\

%s=length(m2);
%q1=input('enter the sequence');
s=randsrc(1,q);
c=complex(ob2,s);

hop1=conj(pb2)';
hob1=conj(c)';
nr3=hop1*hob1';
sum2=nr3'*Wopt1;
%disp('sum2');
%disp(sum2);
cm2=sum2'*sum2;
signal1=mean(cm2);
%disp('cm2');
%disp(cm2);
hop2=conj(pb2)';
hx2=conj(n1)';
nr5=hop2*hx2';
 sum3=nr5'*Wopt1;
 %disp('sum');
 %disp(sum);
  cm3=sum3'*sum3;
  noise1=mean(cm3);
 %disp('cm3');
 %disp(cm3);
 SSINR=signal1/noise1;
 disp('SSINR');
 disp(SSINR);

                %//GRAPH FOR DOA VERSUS SINR\\

    DOA=-80:10:80;
    figure(4);
    plot(DOA,SSINR,'B-*');
    title('simulation for DOA versus SSINR');
    xlabel('DOA in degree');
    ylabel('SINR in db');

                %// BER CALCULATION FOR SB CONFIGURATION\\

       k1=biterr(m2,g2);
       output=k1;
       disp('biterror rate for Sb configuration');
       disp(k1);

                 %//GRAPH FOR NUMBER USERS VERSUS BER\\

     figure(5);
     plot(us,k1,'R-*');

     title('simulation for NUMBER USERS VERSUS BER SYMBOL BASED) ');
     xlabel('NUMBER OF USERS');
     ylabel('BER');


text:
function varargout = spm_mesh_utils(action,varargin)
% A gateway function for surface mesh-related compiled algorithms
%
% FORMAT [N, D] = spm_mesh_utils('neighbours',A)
% Return an array of first-order neighbours given an adjacency matrix
%
% FORMAT Fi = spm_mesh_utils('neighbouringfaces',F,i)
% Return the indices of the neighbouring triangles of a given triangle
%
% FORMAT D = spm_mesh_utils('dijkstra',N,D,i,dmax)
% Compute geodesic distance on a triangular mesh using Dijkstra algorith
%
% FORMAT V = spm_mesh_utils('volume',M)
% Compute the volume of a closed surface mesh
%__________________________________________________________________________
% Copyright (C) 2010-2017 Wellcome Trust Centre for Neuroimaging

% Guillaume Flandin
% $Id: spm_mesh_utils.m 7240 2017-12-19 12:06:59Z guillaume $

%-This is merely the help file for the compiled routine
error('spm_mesh_utils.c not compiled - see Makefile')

author: spm


text:
function gridlines_test03 ( )

%*****************************************************************************80
%
%% GRIDLINES_TEST03 tests GRID_TRIANGULAR.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 July 2010
%
%  Author:
%
%    John Burkardt
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'GRIDLINES_TEST03\n' );
  fprintf ( 1, '  GRID_TRIANGULAR draws a triangular grid.\n' );

  axis ( [ 0.0, 1.0, 0.0, 1.0 ] );

  xval = [ 0.1, 0.95, 0.25 ];
  yval = [ 0.1, 0.2, 0.90 ];
  ngrid = 10;
  grid_triangular ( xval, yval, ngrid );

  xval = [ 0.8, 0.6, 0.7 ];
  yval = [ 0.9, 0.9, 0.7 ];
  ngrid = 3;
  grid_triangular ( xval, yval, ngrid );

  xlabel ( 'X' )
  ylabel ( 'Y' )
  title ( 'Triangular grids' )
  axis equal

  fprintf ( 1, '\n' );
  fprintf ( 1, 'Press return to continue.\n' )
  hold off
  pause
  clf

  return
end


text:
classdef  VideoMaker_TopOpt_levelSet3D_DoubleMirroring < VideoMaker_TopOpt_levelSet
  properties

  end

  methods (Access = public)

        function Make_video_design_variable(obj,output_video_name)
            post = Postprocess_TopOpt_levelSet();
            field2print = post.levelSet_name;
            componentfield = post.levelSet_name_component;
            obj.Make_video_characteristic_function(field2print,componentfield,output_video_name)
        end

  end

  methods (Access = private)
         function Make_video_characteristic_function(obj,field2print,componentfield,output_video_name_in)
            file_tcl_name = 'tcl_gid.tcl';
            file_list = obj.create_file_list(obj.iterations_to_print,obj.file_name,obj.files_folder);

            [output_video_name] = obj.replace_special_character(output_video_name_in);
            [output_photo]=obj.replace_special_character(strcat(obj.files_folder,'\',obj.file_name,'.png'));
            [file_list] = obj.replace_special_character(file_list);

            min_value = -1e-32;

            file_tcl_name_with_path = fullfile(obj.files_folder,file_tcl_name);
            file_path_in = fullfile(pwd,'FEM','PostProcess','Make_Video_characteristic3D_DoubleMirroring.tcl');
            filepath = obj.replace_special_character(file_path_in);

            fid = fopen(file_tcl_name_with_path,'w+');

            fprintf(fid,'GiD_Process PostProcess \n');
            fprintf(fid,'%s\n',['set arg1 "',file_list,'"']);
            fprintf(fid,'%s\n',['set arg2 "',output_video_name,'"']);
            fprintf(fid,'%s\n',['set arg3 "',field2print,'"']);
            fprintf(fid,'%s\n',['set arg4 "',componentfield,'"']);
            fprintf(fid,'%s\n',['set arg5 "',num2str(min_value),'"']);
            fprintf(fid,'%s\n',['set arg6 "',output_photo,'"']);

            fprintf(fid,'%s\n',['source "',filepath,'"']);
            fprintf(fid,'%s\n',['Make_Video_characteristic $arg1 $arg2 $arg3 $arg4 $arg5 $arg6']);
            fprintf(fid,'%s\n',['GiD_Process Mescape Quit']);
            fclose(fid);
            obj.execute_tcl_files(obj.gidPath,file_tcl_name_with_path)
        end

  end
end

author: SwanLab

repo: Swan

save_path: github-repos/MATLAB/SwanLab-Swan


text:
function printx_ppp(x,rtk)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Copyright (C) 2020-2025, by Kai Chen, All rights reserved.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global glc
opt=rtk.opt; nf=rtk.NF;

fprintf('pos   = ');fprintf('%14.5f %14.5f %14.5f',x(1),x(2),x(3));fprintf('\n');

if rtk.NP==9
    fprintf('vel   = ');fprintf('%9.5f%9.5f%9.5f',x(4),x(5),x(6));fprintf('\n');
    fprintf('acc   = ');fprintf('%9.5f%9.5f%9.5f',x(7),x(8),x(9));fprintf('\n');
end

fprintf('clk   = ');
for i=1:glc.NSYS
    fprintf('%9.5f',x(rtk.ic+i));
end
fprintf('\n');

if rtk.NT>0
    fprintf('trop  = ');
    if opt.tropopt==glc.TROPOPT_EST
        fprintf('%9.5f',x(rtk.it+1));
    else
        fprintf('%9.5f',x(rtk.it+2));
        fprintf('%9.5f',x(rtk.it+3));
    end
    fprintf('\n');
end

if rtk.NI>0
    fprintf('iono  = ');
    for i=1:glc.MAXSAT
        if abs(x(rtk.ii+i))>=1e-6,fprintf('%9.5f',x(rtk.ii+i));end
    end
    fprintf('\n');
end

if rtk.ND>0
    fprintf('DCB   = ');
    fprintf('%9.5f',x(rtk.id+1));
    fprintf('\n');
end

if rtk.NB>0
    fprintf('bias1 = ');
    for i=1:glc.MAXSAT
        if abs(x(rtk.ib+i))>=1e-6,fprintf('%10.5f',x(rtk.ib+i));end
    end
    fprintf('\n');
    if nf>=2
        fprintf('bias2 = ');
        for i=1:glc.MAXSAT
            if abs(x(rtk.ib+glc.MAXSAT+i))>=1e-6
                fprintf('%10.5f',x(rtk.ib+glc.MAXSAT+i));
            end
        end
        fprintf('\n');
    end
    if nf>=3
        fprintf('bias3 = ');
        for i=1:glc.MAXSAT
            if abs(x(rtk.ib+glc.MAXSAT*2+i))>=1e-6
                fprintf('%10.5f',x(rtk.ib+glc.MAXSAT*2+i));
            end
        end
        fprintf('\n');
    end
end

return

repo: GINav


text:
function data = imgSTD(img)
    data = struct('intensity',[],'red',[],'green',[],'blue',[],'mean_range',[],'std_range',[],'mean_std',[],'std_std',[],'mean_entropy',[]);
    dim = size(img,1) / 4;
    for j = 1:4
        xend = j*dim;
        xstart = xend - 3;
        for m = 1:4

            intensity = [];
            red = [];
            green = [];
            blue = [];
            mean_range = [];
            std_range = [];
            mean_std = [];
            std_std = [];
            mean_entropy = [];

            yend = m*dim;
            ystart = yend - dim + 1;
            for k = xstart:xend
                for l = ystart:yend
                    intensity = [intensity img{k,l}.intensity];
                    red = [red img{k,l}.red];
                    green = [green img{k,l}.green];
                    blue = [blue img{k,l}.blue];
                    mean_range = [mean_range img{k,l}.mean_range];
                    std_range = [std_range img{k,l}.std_range];
                    mean_std = [mean_std img{k,l}.mean_std];
                    std_std = [std_std img{k,l}.std_std];
                    mean_entropy = [mean_entropy img{k,l}.mean_entropy];
                end
            end
            data.intensity = [data.intensity std(intensity)];
            data.red = [data.red std(red)];
            data.green = [data.green std(green)];
            data.blue = [data.blue std(blue)];
            data.mean_range = [data.mean_range std(mean_range)];
            data.std_range = [data.std_range std(std_range)];
            data.mean_std = [data.mean_std std(mean_std)];
            data.std_std = [data.std_std std(std_std)];
            data.mean_entropy = [data.mean_entropy std(mean_entropy)];
        end
    end
end

author: zhangqianqianQQ

repo: MachineVisionAlgorithm


text:
function [V,t,Err] = evoked(data,Fs,win,width,plt,err)
% Function to calculate the evoked response given continuous data in the
% form time x channels
% Usage [V,t,Err] = evoked(data,Fs,win,width,plt,err)
%
% Inputs
%   Note that all times can be in arbitrary units. But the units have to be
%   consistent. So, if win is in secs, width is in secs and Fs has to be Hz.
%   If win is in samples, so is width and Fs=1.
%
%    data(times, channels/trials or a single vector)      (required)
%    Fs  sampling frequency            (required)
%    win   subsection of data to be used. Default all available data
%    width (s) of smoothing kernel. Default 50 samples
%    plt plot 'n' for no plot, otherwise plot color. Default blue colored lines.
%    err = 0/1. Default 1=calculate bootstrap errorbars.
%
% Outputs
%    V = evoked potential
%    t = times of evaluation
%    Err = bootstrap statdard deviation

if nargin < 2;error('Data, sampling frequency required');end
data=change_row_to_column(data);
N=size(data,1);
data=data';
if nargin <3; win = [0 (N-1)/Fs];end
if nargin <4; width = 50/Fs;end
if nargin <5; plt = 'b';end
if nargin <6;err = 1;end
T=win;
if isempty(T); T = [0 (N-1)/Fs];end
if isempty(width); width = 50/Fs;end
if isempty(plt); plt = 'b';end
if isempty(err);err = 1;end

t = min(T):1/Fs:max(T);
if nargin >= 5
  indx = find(t>T(1) & t<T(2));
  t = t(indx);
  data = data(:,indx);
end

if width > (t(length(t))-t(1))/2
  disp('Width is too large for data segment: should be in seconds')
  disp('Turn off smoothing')
  width = 0;
end

s = t(2)-t(1);
N = fix(width/s);
NT = length(data(:,1));

if NT > 1
    mdata = mean(data);
else
    mdata = data;
end
if N > 4
  smdata = locsmooth(mdata,N,fix(N/2));
else
  smdata = mdata;
end

% if errorbars requested then do a bootstrap over trials...

Err = 0;
if NT < 4;
  disp('Too few trials: no errorbars calculated')
  err = 0;
end

if err ~= 0 && NT > 1
  Nboot = 10;
  bevk = 0;
  sevk = 0;
  for b=1:Nboot
    indx = floor(NT*rand(1,NT)) + 1;
    evktmp = mean(data(indx,:));
    if N > 4
      evktmp = locsmooth(evktmp,N,fix(N/2));
    end
    bevk = bevk + evktmp;
    sevk = sevk + evktmp.^2;
  end
  stdevk = sqrt((sevk/Nboot - bevk.^2/Nboot^2));
  Err = stdevk;
end

V = smdata;
if plt ~= 'n'
  plot(t,smdata,plt)
  hold on
  mn = mean(smdata);
  ax = get(gca,'xlim');
  line(ax,mn*[1 1],'color','k')
  if err
    line(ax,(mn+2*mean(stdevk))*[1 1],'color','r')
    line(ax,(mn-2*mean(stdevk))*[1 1],'color','r')
    hold off
  end
end

author: buzsakilab

repo: buzcode

save_path: github-repos/MATLAB/buzsakilab-buzcode


text:
%PARSEEXPR Parses a string expression.
%
%   [ S_EXPR, CSUBEXPR ] = PARSEEXPR( S_EXPR, S_SUBEXPR, C_ALLOWEXPR )
%   Parses the string expression in S_EXPR and returns a cell array
%   CSUBEXPR with variable expressions to process. Also returns a modified
%   expression S_EXPR which has the sub variable names substituted with
%   'S_SUBEXPR(:,I)' where I is the index in CSUBEXPR to the corresponding
%   variable. This enables a call to EVAL(S_EXPR) after the expressions in
%   CSUBEXPR has been computed and copied to the S_SUBEXPR array. The
%   optional array C_ALLOWEXPR allows strings as subexpressions (which
%   function names would otherwise be removed).
%
%   See also SYMVARBE

% Copyright 2013-2022 Precise Simulation, Ltd.

author: precise-simulation

repo: featool-multiphysics


text:
%#eml
function [min,med,max] = get_median(inbuf)

numpixels = length(inbuf);

tbuf = inbuf;

for ii=eml.unroll(1:numpixels)
    if bitand(ii,uint32(1)) == 1
        tbuf = compare_stage1(tbuf);
    else
        tbuf = compare_stage2(tbuf);
    end
end

min = tbuf(1);
med = tbuf(ceil(numpixels/2));
max = tbuf(numpixels);

end

function outbuf = compare_stage1(inbuf)
numpixels = length(inbuf);
tbuf = compare_stage(inbuf(1:numpixels-1));
outbuf = [tbuf(:)' inbuf(numpixels)];
end

function outbuf = compare_stage2(inbuf)
numpixels = length(inbuf);
tbuf = compare_stage(inbuf(2:numpixels));
outbuf = [inbuf(1) tbuf(:)'];
end

function [outbuf] = compare_stage(inbuf)

step = 2;
numpixels = length(inbuf);

outbuf = inbuf;

for ii=eml.unroll(1:step:numpixels)
    t = compare_pixels([inbuf(ii), inbuf(ii+1)]);
    outbuf(ii) = t(1);
    outbuf(ii+1) = t(2);
end

end

function outbuf = compare_pixels(inbuf)
if (inbuf(1) > inbuf(2))
    outbuf = [inbuf(1), inbuf(2)];
else
    outbuf = [inbuf(2), inbuf(1)];
end
end

author: Sable


text:
function m = getmem
% return total system memory in kb

try
  m = memory;
  m = m.MemAvailableAllArrays / 1024;
catch %#ok<CTCH>
  [r,s] = system('free');
  m = sscanf(s(strfind(s,'Mem:')+5:end),'%d',3);
  if isempty(m)
    m = sscanf(s(strfind(s,'Speicher:')+9:end),'%d',3);
  end
  if isempty(m)
    m = 300 * 1024;
  else
    m = m(3);
  end
  % reset error
  lasterr('');
end
end


text:
%RPY2JAC Jacobian from RPY angle rates to angular velocity
%
% J = RPY2JAC(EUL) is a Jacobian matrix (3x3) that maps roll-pitch-yaw angle
% rates to angular velocity at the operating point RPY=[R,P,Y].
%
% J = RPY2JAC(R, P, Y) as above but the roll-pitch-yaw angles are passed
% as separate arguments.
%
% Notes::
% - Used in the creation of an analytical Jacobian.
%
% See also EUL2JAC, SerialLink.JACOBN.

% Copyright (C) 1993-2015, by Peter I. Corke
%
% This file is part of The Robotics Toolbox for MATLAB (RTB).
%
% RTB is free software: you can redistribute it and/or modify
% it under the terms of the GNU Lesser General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% RTB 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 Lesser General Public License for more details.
%
% You should have received a copy of the GNU Leser General Public License
% along with RTB.  If not, see <http://www.gnu.org/licenses/>.
%
% http://www.petercorke.com

function J = rpy2jac(r, p, y)

    if length(r) == 3
        % rpy2jac([r,p,y])
        p = r(2);
        y = r(3);
        r = r(1);
    elseif nargin ~= 3
        error('RTB:rpy2jac:badarg', 'bad arguments');
    end
	J = [
        1  0       sin(p)
        0  cos(r)  -cos(p)*sin(r)
        0  sin(r)  cos(p)*cos(r)
        ];

author: Allopart

repo: rbpf-gmapping


text:
function pdf = dirichlet_pdf ( x, n, a )

%*****************************************************************************80
%
%% DIRICHLET_PDF evaluates the Dirichlet PDF.
%
%  Discussion:
%
%    PDF(X)(N,A) = Product ( 1 <= I <= N ) X(I)**( A(I) - 1 )
%      * Gamma ( A_SUM ) / A_PROD
%
%    where
%
%      0 < A(I) for all I;
%      0 <= X(I) for all I;
%      Sum ( 1 <= I <= N ) X(I) = 1;
%      A_SUM = Sum ( 1 <= I <= N ) A(I).
%      A_PROD = Product ( 1 <= I <= N ) Gamma ( A(I) )
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license.
%
%  Modified:
%
%    30 October 2004
%
%  Author:
%
%    John Burkardt
%
%  Parameters:
%
%    Input, real X(N), the argument of the PDF.  Each X(I) should
%    be greater than 0.0, and the X(I)'s must add up to 1.0.
%
%    Input, integer N, the number of components.
%
%    Input, real A(N), the probabilities for each component.
%    Each A(I) should be positive.
%
%    Output, real PDF, the value of the PDF.
%
  tol = 0.0001;

  for i = 1 : n
    if ( x(i) <= 0.0 )
      fprintf ( 1, '\n' );
      fprintf ( 1,  'DIRICHLET_PDF - Fatal error!\n' );
      fprintf ( 1,  '  X(I) <= 0.\n' );
    end
  end

  x_sum = sum ( x(1:n) );

  if ( tol < abs ( x_sum - 1.0 ) )
    fprintf ( 1, '\n' );
    fprintf ( 1,  'DIRICHLET_PDF - Fatal error!\n' );
    fprintf ( 1,  '  SUM X(I) =/= 1.\n' );
  end

  a_sum = sum ( a(1:n) );

  a_prod = 1.0;
  for i = 1 : n
    a_prod = a_prod * gamma ( a(i) );
  end

  pdf = gamma ( a_sum ) / a_prod;
  for i = 1 : n
    pdf = pdf * x(i)^( a(i) - 1.0 );
  end

  return
end

author: johannesgerer

repo: jburkardt-m


text:
% Copyright (C) 1994, 1995, 1996, 1997, 1998, 2000, 2002, 2004, 2005,
%               2006, 2007, 2008, 2009 John W. Eaton
%
% This file is part of Octave.
%
% Octave 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.
%
% Octave 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 Octave; see the file COPYING.  If not, see
% <http://www.gnu.org/licenses/>.

% -*- texinfo -*-
% @deftypefn {Function File} {} postpad (@var{x}, @var{l}, @var{c})
% @deftypefnx {Function File} {} postpad (@var{x}, @var{l}, @var{c}, @var{dim})
% @seealso{prepad, resize}
% @end deftypefn

% Author: Tony Richardson <arichard@stark.cc.oh.us>
% Created: June 1994

function y = postpad (x, l, c, dim)

if nargin < 2 || nargin > 4
  %print_usage ();
  error('wrong number of input arguments, should be between 2 and 4');
end

if nargin < 3 || isempty(c)
  c = 0;
else
  if ~isscalar(c)
    error ('postpad: third argument must be empty or a scalar');
  end
end

nd = ndims(x);
sz = size(x);
if nargin < 4
  % Find the first non-singleton dimension
  dim  = 1;
  while dim < nd+1 && sz(dim)==1
    dim = dim + 1;
  end
  if dim > nd
    dim = 1;
  elseif ~(isscalar(dim) && dim == round(dim)) && dim > 0 && dim< nd+1
    error('postpad: dim must be an integer and valid dimension');
  end
end

if ~isscalar(l) || l<0
  error ('second argument must be a positive scalar');
end

if dim > nd
  sz(nd+1:dim) = 1;
end

d = sz(dim);

if d >= l
  idx = cell(1,nd);
  for i = 1:nd
    idx{i} = 1:sz(i);
  end
  idx{dim} = 1:l;
  y = x(idx{:});
else
  sz(dim) = l-d;
  y = cat(dim, x, c * ones(sz));
end


text:
function [vw, slices] = selectSliceCoords(vw)
%
%    [vw, selectedPositions] = selectSliceCoords(vw)
%
% Allow the user to use the mouse to set the three sliceNum editable text
% fields in the Volume view.
%
% The current selected position is stored in vw.ui.sliceNumFields.
% The values are also returned as the 2nd argument.
%

if ~strcmp(vw.viewType,'Volume') & ~strcmp(vw.viewType,'Gray')
  myErrorDlg('selectSliceCoords only for Volume view.');
end

% Get mouse input
[x,y] = ginput(1);
x=round(x);
y=round(y);

% Interpret mouse click according to current slice orientation
curSlice = viewGet(vw, 'Current Slice');
sliceOri=getCurSliceOri(vw);
switch sliceOri
  case 1				% axi (y=cor pos, x=sag pos)
    slices=[curSlice y x];
  case 2 				% cor (y=axi pos, x=sag pos)
    slices=[y curSlice x];
  case 3 				% sag (y=axi pos, x=cor pos)
    slices=[y x curSlice];
end

% Set the sliceNum editable text fields
for n=1:3
  set(vw.ui.sliceNumFields(n),'String',num2str(slices(n)));
end

return;

% Debug/Test

[VOLUME{1},selectedPosition] = selectSliceCoords(VOLUME{1});

author: vistalab

repo: vistasoft

save_path: github-repos/MATLAB/vistalab-vistasoft


text:
function [s] = gsp_filter_synthesis(G, filters, c, param)
%GSP_FILTER_SYNTHESIS Synthesis operator of a gsp filterbank
%   Usage:  s = gsp_filter_synthesis(G, filters, c);
%           s = gsp_filter_synthesis(G, filters, c, param);
%
%   Input parameters:
%         G         : Graph structure.
%         filters   : Set of spectral graph filters.
%         c         : Transform coefficients
%         param     : Optional parameter
%   Output parameters:
%         signal    : sythesis signal
%
%   'gsp_filter_synthesis(G,filters,c)' computes the synthesis
%   operator for coefficient $c$, where the atoms of the transform
%   dictionary are generalized translations of each graph spectral filter
%   to each vertex on the graph.
%
%   .. f = D * c
%
%   .. math:: f =  D c
%
%   where the columns of $D$ are $g_{i,m}=T_i g_m$, and $T_i$ is a
%   generalized translation operator applied to each filter
%   $\hat{g}_m(\cdot)$.
%
%   Each column of *c* is the response of the signal to one filter.
%
%   Example:::
%
%         Nf = 4;
%         G = gsp_sensor(30);
%         G = gsp_estimate_lmax(G);
%         G = gsp_estimate_lmax(G);
%         g = gsp_design_mexican_hat(G, Nf);
%         f = zeros(G.N,1);
%         f(1) = 1;
%         f = G.L^2*f;
%         ff = gsp_filter_analysis(G,g,f);
%         f2 = gsp_filter_synthesis(G,g,ff);
%         paramplot.show_edges = 1;
%         figure()
%         subplot(211)
%         gsp_plot_filter(G,g)
%         subplot(223)
%         gsp_plot_signal(G,f,paramplot);
%         subplot(224)
%         gsp_plot_signal(G,f2,paramplot);
%
%   Additional parameters
%   ---------------------
%
%   * *param.method*  : Select the method to be used for the computation.
%     * 'exact'     : Exact method using the graph Fourier matrix
%     * 'cheby'     : Chebyshev polynomial approximation
%     * 'lanczos'   : Lanczos approximation
%     Default: if the Fourier matrix is present: 'exact' otherwise 'cheby'
%   * *param.order* : Degree of the Chebyshev approximation
%     (default=30).
%   * *param.verbose* : Verbosity level (0 no log - 1 display warnings)
%     (default 1).
%
%   See also: gsp_filter_analysis gsp_filter_inverse
%
%   References: hammond2011wavelets
%

% Author: Nathanael Perraudin
% Testing: test_filter
% Date: 19 March 2014

% TODO: Perfect

% Read input parameters
if nargin < 4
    param = struct;
end

if iscell(G)
    NG = numel(G);
    s = cell(NG,1);
    for ii = 1:NG
        warning('Check what happen here')
       s{ii} = gsp_filter_synthesis(G{ii}, filters{ii}, c{ii}, param);
%         if iscell(s)
%             c{ii} = gsp_filter_analysis(G{ii}, fi{ii}, s{ii}, param);
%         else
%             c{ii} = gsp_filter_analysis(G{ii}, fi{ii}, s, param);
%         end
    end
    return
end

if isnumeric(filters)
    Nf = size(filters,2);
else
    Nf = numel(filters);
end

if isfield(param, 'exact')
    warning('param.exact is not used anymore. Please use param.method instead');
    if param.exact
        param.method = 'exact';
    else
        param.method = 'cheby';
    end
end

if ~isfield(param,'method')
    if gsp_check_fourier(G)
        param.method = 'exact';
    else
        param.method = 'cheby';
    end
end

if ~isfield(param,'order'); param.order = 30; end
if ~isfield(param,'verbose'); param.verbose = 1; end

if isfield(param, 'cheb_order')
    param.order = param.cheb_order;
    warning('param.cheb_order is not used anymore. Please use param.order instead');
end

switch param.method
    case 'exact'

        if ~gsp_check_fourier(G)
            if param.verbose
                warning(['GSP_FILTER_SYNTHESIS: The Fourier matrix is not ',...
                    'available. The function will compute it for you. ',...
                    'However, if you apply many time this function, you ',...
                    'should precompute it using the function: ',...
                    'gsp_compute_fourier_basis']);
            end
            G=gsp_compute_fourier_basis(G);
        end
        if isnumeric(filters)
            fie = filters;
        else
            fie = gsp_filter_evaluate(filters,G.e);
        end
%         Nv = size(c,2);
%         s =zeros(G.N,size(c,2));
%         for ii=1:Nf
%             s = s + G.U * ...
%                 (repmat(fie(:,ii),1,Nv) ...
%                 .* (G.U' * c((1:G.N)+G.N * (ii-1),:)));
%         end

        chat = gsp_gft(G,gsp_vec2mat(c,numel(filters)));
        shat = squeeze(sum(bsxfun(@times, fie, chat), 2));
        s = gsp_igft(G, shat);

    case 'cheby'
        if ~isfield(G,'lmax');
            G = gsp_estimate_lmax(G);
            if param.verbose
                warning(['GSP_FILTER_ANALYSIS: The variable lmax is not ',...
                    'available. The function will compute it for you. ',...
                    'However, if you apply many time this function, you ',...
                    'should precompute it using the function: ',...
                    'gsp_estimate_lmax']);
            end
        end

        cheb_coeffs = gsp_cheby_coeff(G, filters,...
                param.order, param.order +1);

        s=zeros(G.N,size(c,2));

        for ii=1:Nf
            s = s + gsp_cheby_op(G,cheb_coeffs(:,ii),c((1:G.N)+G.N * (ii-1),:));
        end

    case 'lanczos'
        s=zeros(G.N,size(c,2));
        if ~iscell(filters)
            filters = {filters};
        end
        for ii=1:Nf
            s = s + gsp_lanczos_op(G, filters{ii}, c((1:G.N)+G.N * (ii-1),:), param);
        end

    otherwise
        error('Unknown method: please select exact, cheby or lanczos');
end

end

repo: gspbox


text:
%% Defines the class interface for the ECGtask derived classes
%
% Description:
% Abstract class for defining ECGtask interface
%
% Arguments: %
%
% Output:
%
% Examples:
%
%
% See also ECGtask_QRS_detection, ECGtask_ECG_delineation
%
% Author: Mariano Llamedo Soria llamedom@electron.frba.utn.edu.ar
% Version: 0.1 beta
% Last update: 14/5/2014
% Birthdate  : 20/3/2012
% Copyright 2008-2015
%
classdef ECGtask < handle

    properties(Abstract = true, GetAccess = public, Constant)
        name
        target_units
        doPayload
        % if user = memory;
        % memory_constant is the fraction respect to user.MaxPossibleArrayBytes
        % which determines the maximum input data size.

    end

    properties(Abstract = true, GetAccess = public, SetAccess = private)
        memory_constant
        started
    end

    properties(Abstract = true)
        progress_handle
        tmp_path
    end

    methods (Abstract)

% ECG_header is the header for this signal
% ECG_annotations are the QRS locations available for this signal
        function Start(obj, ECG_header, ECG_annotations)
        endfunction
% ECG is the ECG signal
% ECG_start_offset is the location of ECG(1,:) within the whole signal
% ECG_sample_start_end_idx are the start and end samples within ECG to
%    generate valid results.
% ECG_header is the header for this signal
% ECG_annotations are the QRS locations available for this signal
% ECG_annotations_start_end_idx are the start and end indexes corresponding
%    to the first and last element of ECG_annotations in the current
%    iteration.
        function payload = Process(ECG, ECG_start_offset, ECG_sample_start_end_idx, ECG_header, ECG_annotations, ECG_annotations_start_end_idx )
        endfunction
% payload is the result of the task, after the final processing done by this method
% ECG_header is the header for this signal
        function payload = Finish(obj, payload, ECG_header)
        endfunction

% plA, plB are two payloads produced by the Process method.
        function payload = Concatenate(plA, plB)
        endfunction

    end

end

author: marianux

repo: ecg-kit

save_path: github-repos/MATLAB/marianux-ecg-kit


text:
function [pathname] = uigetdir2(start_path, dialog_title, selections)
% Pick multiple directories and/or files

import javax.swing.JFileChooser;

if nargin == 0 || isempty(start_path) || isnumeric(start_path) % Allow a null argument.
    start_path = pwd;
end

jchooser = javaObjectEDT('javax.swing.JFileChooser', start_path);

jchooser.setFileSelectionMode(JFileChooser.FILES_AND_DIRECTORIES);
if nargin > 1
    jchooser.setDialogTitle(dialog_title);
end

jchooser.setMultiSelectionEnabled(true);

% set selected file
if nargin > 2,
  if ~iscell(selections),
    selections = {selections};
  end
  selections = setdiff(selections,'');
  selections(cellfun(@isempty,selections)) = [];
  if ~isempty(selections),
    clear jselections;
    for i = 1:numel(selections),
      jselections(i) = javaObjectEDT('java.io.File',selections{i}); %#ok<AGROW>
    end
    jchooser.setSelectedFiles(jselections);
  end
end

status = jchooser.showOpenDialog([]);

if status == JFileChooser.APPROVE_OPTION
  jFile = jchooser.getSelectedFiles();
	pathname{size(jFile, 1)}=[];
  for i=1:size(jFile, 1)
		pathname{i} = char(jFile(i).getAbsolutePath);
	end

elseif status == JFileChooser.CANCEL_OPTION
    pathname = [];
else
    error('Error occured while picking file.');
end

author: kristinbranson

repo: JAABA


text:
function w = simplex_projection(v, b)
% Simplex_Projection Projects point onto simplex of specified radius.
%
% w = simplex_projection(v, b) returns the vector w which is the solution
%   to the following constrained minimization problem:
%
%    min   ||w - v||_2
%    s.t.  sum(w) <= b, w >= 0.
%
%   That is, performs Euclidean projection of v to the positive simplex of
%   radius b.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This file is part of OLPS: http://OLPS.stevenhoi.org/
% Original authors: John Duchi (jduchi@cs.berkeley.edu)
% Contributors: Bin LI, Steven C.H. Hoi
% Change log:
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if (b < 0)
  error('Radius of simplex is negative: %2.3f\n', b);
end
v = (v > 0) .* v;
u = sort(v,'descend');
sv = cumsum(u);
rho = find(u > (sv - b) ./ (1:length(u))', 1, 'last');
theta = max(0, (sv(rho) - b) / rho);
w = max(v - theta, 0);

author: OLPS

repo: OLPS

save_path: github-repos/MATLAB/OLPS-OLPS


text:
function DoMultiEchoChk

global VCtl

if VCtl.TEPerTR ==1
    return;
end

if ~isfield(VCtl,'ME_TEs')
    error('Please load MultiEcho tab. Multiple TE values must be provided before proceeding.');
end

if VCtl.TEPerTR ~= length(VCtl.ME_TEs)
    error('The number of TE values must match TEPerTR.');
end

end

author: leoliuf

repo: MRiLab

save_path: github-repos/MATLAB/leoliuf-MRiLab


text:
function TestDereverbNet_Regression_Rpi
addpath('local');
dnn1 = load('nnet/Dereverb.noCMN.DeltaByEqn.MbSize20.U52013.771-LSTM-2048-771.L2_3E-4.LR_5E-3/nnet.itr10.LR4.39E-4.CV1630.073.mat');
dnn2 = load('nnet/DereverbMask.noCMN.DeltaByEqn.MbSize20.U52013.771-LSTM-1500-771.L2_3E-4.LR_5E-3/nnet.itr11.LR1.46E-4.CV1801.505.mat');

[layer1, para1] = PrepareProcessing(dnn1);
[layer2, para2] = PrepareProcessing(dnn2);

wavroot = 'D:\Data\RpiRecording';
wavfiles = findFiles(wavroot, 'wav');

for i=2:length(wavfiles)
    [wav, fs] = audioread(wavfiles{i});
    wav = StoreWavInt16(wav(:,2))';
    DataTmp(1).data{1} = wav(fs*10+1: min(fs*20,length(wav)) );

    output{1} = FeatureTree2(DataTmp, para1, layer1);
    output{2} = FeatureTree2(DataTmp, para2, layer2);

    for j = 1:2
        noisyLogSpec{j} = gather(output{j}{1}{1});
        enhancedLogSpec{j} = gather(output{j}{1}{2});
        noisySTFT{j} = gather(output{j}{1}{3});
        phase = angle(noisySTFT{j});
        noisy_wav = abs2wav(exp(noisyLogSpec{j}(1:257,:)/2)', phase', 400, 240);
        enhanced_wav{j} = abs2wav(exp(enhancedLogSpec{j}(1:257,:)/2)', phase', 400, 240);
    end
    imagesc(CMN([noisyLogSpec{1}; enhancedLogSpec{1}(1:257,:); enhancedLogSpec{2}(1:257,:)]')');
    title(regexprep(wavfiles{i}, '_', '\\_'));  pause(0.1);
    soundsc(noisy_wav, 16000); pause
    soundsc(enhanced_wav{1}, 16000);                pause
    soundsc(enhanced_wav{2}, 16000);                pause
end

end

function [layer, para] = PrepareProcessing(dnn)

layer = dnn.layer;
para = dnn.para;
para.topology.useFileName = 1;
para.local.seglen = 100;
para.local.segshift = 100;
para.useGPU = 0;

input_idx = ReturnLayerIdxByName(layer, 'input');
layer = layer(1:input_idx(end)-1);
para.IO = RemoveIOStream(para.IO, 2);
if isfield(para.topology, 'useCMN') && para.topology.useCMN==0
    log_idx = ReturnLayerIdxByName(layer, 'Log');
    para.out_layer_idx = [log_idx(1) length(layer)];
else
    cmn_idx = ReturnLayerIdxByName(layer, 'CMN');
    para.out_layer_idx = [cmn_idx(1) length(layer)];
end
stft_idx = ReturnLayerIdxByName(layer, 'stft');
para.out_layer_idx = [para.out_layer_idx stft_idx(1) 1];

reverb_root = ChoosePath4OS({'D:/Data/REVERB_Challenge', '/media/xiaoxiong/OS/data1/G/REVERB_Challenge'});   % you can set two paths, first for windows OS and second for Linux OS.
para.local.wavroot = [reverb_root];
para.local.wsjcam0root = ChoosePath4OS({'D:/Data/wsjcam0', '/media/xiaoxiong/OS/data1/G/wsjcam0_wav'});
para.local.wsjcam0ext = ChoosePath4OS({'wav', 'wv1'});
para.local.useFileName = 1;      % if set to 0, load all training data to memory. otherwise, only load file names.
para.local.loadLabel = 0;

end

author: singaxiong

repo: SignalGraph


text:
function mean = pareto_mean ( a, b )

%*****************************************************************************80
%
%% PARETO_MEAN returns the mean of the Pareto 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 MEAN, the mean of the PDF.
%
  if ( b <= 1.0 )
    fprintf ( 1, '\n' );
    fprintf ( 1, 'PARETO_MEAN - Fatal error!\n' );
    fprintf ( 1, '  For B <= 1, the mean does not exist.\n' );
    mean = 0.0;
    return
  end

  mean = b * a / ( b - 1.0 );

  return
end