Commit d92d1364 authored by Ling Sun's avatar Ling Sun
Browse files

initialization

parents
function [b,h] = binary_fstat(f,f0,P,a0)
% Phase - Ph = f0*(t+a0*sin(2*pi*t/P));
% Frequency - F = f0+2*pi*a0*f0/P*cos(2*pi*t/P);
% Fstat width = [-c c]/P where c --- frequency magnitude (Carson's rule)
% f0 is fixed to setup matched filter for each 1-Hz sub-band
c=(2*pi*a0*f0);
% Make it slightly wider than c
N=ceil(c)+100;
n=-N:N;
b0=besselj(n,c);
M=length(f);
% Identify the indices to place the matched filter in the 1-Hz band
% k indicates the interpolated indices at query points n/P using nearest interpolation, 1 is the extrapval for out of range values
k=interp1(f-mean(f),1:M,n/P,'nearest',1);
% Initialise the 1-Hz band matched filter
b=zeros(M,1);
% Place the filter
b(k)=abs(b0).^2;
% h is comb matched filter, not used in the search
h=zeros(M,1);
h(k)=1;
function j=argmax(L);
[~,j]=max(L);
\ No newline at end of file
function [mb,nb] = bestblk(siz,k)
%BESTBLK Block size with minimum padding.
% SIZ = BESTBLK([M N],K) returns a block size which minimizes the padding
% required along the outer partial blocks, while block processing an M by
% N image. K is a scalar specifying the maximum row and column dimensions
% for the block; if the argument is omitted, it defaults to 100. SIZ is a
% 1-by-2 vector containing row and column dimensions for the block.
%
% [MB,NB] = BESTBLK([M N],K) returns the row and column dimensions in MB
% and NB, respectively.
%
% Example
% -------
% siz = bestblk([640 800], 72)
%
% See also BLOCKPROC.
% Copyright 1993-2011 The MathWorks, Inc.
if nargin==1, k = 100; end % Default block size
%
% Find possible factors of siz that make good blocks
%
% Define acceptable block sizes
m = floor(k):-1:floor(min(ceil(siz(1)/10),ceil(k/2)));
n = floor(k):-1:floor(min(ceil(siz(2)/10),ceil(k/2)));
% Choose that largest acceptable block that has the minimum padding.
[~,ndx] = min(ceil(siz(1)./m).*m-siz(1)); blk(1) = m(ndx);
[~,ndx] = min(ceil(siz(2)./n).*n-siz(2)); blk(2) = n(ndx);
if nargout == 2,
mb = blk(1);
nb = blk(2);
else
mb = blk;
end
function a = col2im(b,block,mat,kind)
%COL2IM Rearrange matrix columns into blocks.
% A = COL2IM(B,[M N],[MM NN],'distinct') rearranges each column of B into a
% distinct M-by-N block to create the matrix A of size MM-by-NN. If B =
% [A11(:) A21(:) A12(:) A22(:)], where each column has length M*N, then A =
% [A11 A12; A21 A22] where each Aij is M-by-N.
%
% A = COL2IM(B,[M N],[MM NN],'sliding') rearranges the row vector B into a
% matrix of size (MM-M+1)-by-(NN-N+1). B must be a vector of size
% 1-by-(MM-M+1)*(NN-N+1). B is usually the result of processing the output
% of IM2COL(...,'sliding') using a column compression function (such as
% SUM).
%
% COL2IM(B,[M N],[MM NN]) is the same as COL2IM(B,[M N],[MM NN],'sliding').
%
% Class Support
% -------------
% B can be logical or numeric. A has the same class as B.
%
% Example
% -------
% B = reshape(uint8(1:25),[5 5])'
% C = im2col(B,[1 5])
% A = col2im(C,[1 5],[5 5],'distinct')
%
% See also BLOCKPROC, COLFILT, IM2COL, NLFILTER.
% Copyright 1993-2015 The MathWorks, Inc.
% I/O Spec
% ========
% IN
% B - any numeric class or logical
% M,N,MM,NN - double, integer
% OUT
% A - same class as B
validateattributes(b,{'numeric' 'logical'},{'nonsparse'},mfilename,'B',1);
validateattributes(block,{'double'},{'integer' 'real' 'positive'},...
mfilename,'[M N]',2);
validateattributes(mat,{'double'},{'integer' 'real' 'positive'},...
mfilename,'[MM NN]',3);
if nargin < 4, % Try to determine which block type is assumed.
kind = 'sliding';
end
if ~ischar(kind),
error(message('images:col2im:wrongBlockType'));
end
kind = [lower(kind) ' ']; % Protect against short string
if kind(1)=='d', % Distinct
% Check argument sizes
[m,n] = size(b);
if prod(block)~=m, error(message('images:col2im:wrongSize')); end
% Find size of padded A.
mpad = rem(mat(1),block(1)); if mpad>0, mpad = block(1)-mpad; end
npad = rem(mat(2),block(2)); if npad>0, npad = block(2)-npad; end
mpad = mat(1)+mpad; npad = mat(2)+npad;
if mpad*npad/prod(block)~=n,
error(message('images:col2im:inconsistentSize'));
end
mblocks = mpad/block(1);
nblocks = npad/block(2);
aa = mkconstarray(class(b), 0, [mpad npad]);
x = mkconstarray(class(b), 0, block);
rows = 1:block(1); cols = 1:block(2);
for i=0:mblocks-1,
for j=0:nblocks-1,
x(:) = b(:,i+j*mblocks+1);
aa(i*block(1)+rows,j*block(2)+cols) = x;
end
end
a = aa(1:mat(1),1:mat(2));
elseif kind(1)=='s', % sliding
a = reshape(b,mat(1)-block(1)+1,mat(2)-block(2)+1);
else
error(message('images:col2im:unknownBlockType', deblank(kind)))
end
function b = colfilt(varargin)
%COLFILT Columnwise neighborhood operations.
% COLFILT processes distinct or sliding blocks as columns. COLFILT can
% perform similar operations to BLOCKPROC and NLFILTER, but often
% executes much faster.
%
% B = COLFILT(A,[M N],BLOCK_TYPE,FUN) processes the image A by rearranging
% each M-by-N block of A into a column of a temporary matrix, and then
% applying the function FUN to this matrix. FUN must be a
% FUNCTION_HANDLE. COLFILT zero pads A, if necessary.
%
% Before calling FUN, COLFILT calls IM2COL to create the temporary
% matrix. After calling FUN, COLFILT rearranges the columns of the matrix
% back into M-by-N blocks using COL2IM.
%
% BLOCK_TYPE is a string with one these values:
%
% 'distinct' for M-by-N distinct blocks
% 'sliding' for M-by-N sliding neighborhoods
%
% B = COLFILT(A,[M N],'distinct',FUN) rearranges each M-by-N distinct
% block of A in a temporary matrix, and then applies the function FUN to
% this matrix. FUN must return a matrix of the same size as the temporary
% matrix. COLFILT then rearranges the columns of the matrix returned by
% FUN into M-by-N distinct blocks.
%
% B = COLFILT(A,[M N],'sliding',FUN) rearranges each M-by-N sliding
% neighborhood of A into a column in a temporary matrix, and then applies
% the function FUN to this matrix. FUN must return a row vector containing
% a single value for each column in the temporary matrix. (Column
% compression functions such as SUM return the appropriate type of
% output.) COLFILT then rearranges the vector returned by FUN into a
% matrix of the same size as A.
%
% B = COLFILT(A,[M N],[MBLOCK NBLOCK],BLOCK_TYPE,FUN) processes the
% matrix A as above, but in blocks of size MBLOCKS-by-NBLOCKS to save
% memory. Note that using the [MBLOCK NBLOCK] argument does not change the
% result of the operation.
%
% B = COLFILT(A,'indexed',...) processes A as an indexed image, padding
% with zeros if the class of A is uint8 or uint16, or ones if the class of
% A is double or single.
%
% Note
% ----
% To save memory, COLFILT may divide A into subimages and process one
% subimage at a time. This implies that FUN may be called multiple
% times, and that the first argument to FUN may have a different number
% of columns each time.
%
% Class Support
% -------------
% The input image A can be of any class supported by FUN. The class of B
% depends on the class of the output from FUN.
%
% Example
% -------
% I = imread('tire.tif');
% figure, imshow(I)
% I2 = uint8(colfilt(I,[5 5],'sliding',@mean));
% figure, imshow(I2)
%
% See also BLOCKPROC, COL2IM, FUNCTION_HANDLE, IM2COL, NLFILTER.
% Copyright 1993-2011 The MathWorks, Inc.
% Obsolete syntax:
% B = COLFILT(A,[M N],BLOCK_TYPE,FUN,P1,P2,...) passes the additional
% parameters P1,P2,..., to FUN. COLFILT calls FUN using:
%
% Y = FUN(X,P1,P2,...)
% I/O Spec
% ========
% IN
% A - any class supported by FUN
% M,N - double, integer scalars
% FUN - function handle
% BLOCK_TYPE - string; either 'sliding' or 'distinct'
% MBLOCK,NBLOCK - doube, integer scalars
% OUT
% B - class of output depends on the class produced
% by FUN
narginchk(4,Inf)
a = varargin{1};
if ischar(varargin{2}) && strcmpi(varargin{2},'indexed')
indexed = 1; % It's an indexed image
varargin(2) = [];
else
indexed = 0;
end
nhood = varargin{2};
if nargin==4,
% COLFILT(A, [m n], kind, fun)
kind = varargin{3};
fun = varargin{4};
block = bestblk(size(a));
params = cell(0,0);
else
if ~ischar(varargin{3}),
% COLFILT(A, [m n], block, kind, fun, P1, ...)
block = varargin{3};
kind = varargin{4};
fun = varargin{5};
params = varargin(6:end);
else
% COLFILT(A, [m n], kind, fun, P1, ...)
kind = varargin{3};
fun = varargin{4};
block = bestblk(size(a));
params = varargin(5:end);
end
end
fun = fcnchk(fun, length(params));
if ~ischar(kind),
error(message('images:colfilt:invalidBlockType'))
end
kind = [lower(kind) ' ']; % Protect against short string
if kind(1)=='s', % Sliding
if all(block>=size(a)), % Process the whole matrix at once.
% Expand A
[ma,na] = size(a);
if indexed && isa(a, 'double'),
aa = ones(size(a)+nhood-1);
else
aa = mkconstarray(class(a), 0, (size(a)+nhood-1));
end
aa(floor((nhood(1)-1)/2)+(1:ma),floor((nhood(2)-1)/2)+(1:na)) = a;
% Convert neighborhoods of matrix A to columns
x = im2col(aa,nhood,'sliding');
% Apply fun to column version of a
b = reshape(feval(fun,x,params{:}), size(a));
else % Process the matrix in blocks of size BLOCK.
% Expand A: Add border, pad if size(a) is not divisible by block.
[ma,na] = size(a);
mpad = rem(ma,block(1)); if mpad>0, mpad = block(1)-mpad; end
npad = rem(na,block(2)); if npad>0, npad = block(2)-npad; end
if indexed && isa(a, 'double'),
aa = ones(size(a) + [mpad npad] + (nhood-1));
else
aa = mkconstarray(class(a), 0, ...
size(a) + [mpad npad] + (nhood-1));
end
aa(floor((nhood(1)-1)/2)+(1:ma),floor((nhood(2)-1)/2)+(1:na)) = a;
%
% Process each block in turn.
%
m = block(1) + nhood(1)-1;
n = block(2) + nhood(2)-1;
mblocks = (ma+mpad)/block(1);
nblocks = (na+npad)/block(2);
% Figure out return type
chunk = a(1:nhood(1), 1:nhood(2));
temp = feval(fun, chunk(:), params{:});
b = mkconstarray(class(temp), 0, [ma+mpad,na+npad]);
arows = (1:m); acols = (1:n);
brows = (1:block(1)); bcols = (1:block(2));
mb = block(1); nb = block(2);
for i=0:mblocks-1,
for j=0:nblocks-1,
x = im2col(aa(i*mb+arows,j*nb+acols),nhood);
b(i*mb+brows,j*nb+bcols) = ...
reshape(feval(fun,x,params{:}),block(1),block(2));
end
end
b = b(1:ma,1:na);
end
elseif kind(1)=='d', % Distinct
if all(block>=size(a)), % Process the whole matrix at once.
% Convert neighborhoods of matrix A to columns
x = im2col(a,nhood,'distinct');
% Apply fun to column version of A and reshape
b = col2im(feval(fun,x,params{:}),nhood,size(a),'distinct');
else % Process the matrix in blocks of size BLOCK.
% Expand BLOCK so that it is divisible by NHOOD.
mpad = rem(block(1),nhood(1)); if mpad>0, mpad = nhood(1)-mpad; end
npad = rem(block(2),nhood(2)); if npad>0, npad = nhood(2)-npad; end
block = block + [mpad npad];
% Expand A: Add border, pad if size(A) is not divisible by BLOCK.
[ma,na] = size(a);
mpad = rem(ma,block(1)); if mpad>0, mpad = block(1)-mpad; end
npad = rem(na,block(2)); if npad>0, npad = block(2)-npad; end
if indexed && isa(a, 'double'),
aa = ones(size(a) + [mpad npad]);
else
aa = mkconstarray(class(a), 0, size(a) + [mpad npad]);
end
aa((1:ma),(1:na)) = a;
%
% Process each block in turn.
%
m = block(1); n = block(2);
mblocks = (ma+mpad)/block(1);
nblocks = (na+npad)/block(2);
% Figure out return type
chunk = a(1:nhood(1), 1:nhood(2));
temp = feval(fun, chunk(:), params{:});
b = mkconstarray(class(temp), 0, [ma+mpad,na+npad]);
rows = 1:block(1); cols = 1:block(2);
for i=0:mblocks-1,
ii = i*m+rows;
for j=0:nblocks-1,
jj = j*n+cols;
x = im2col(aa(ii,jj),nhood,'distinct');
b(ii,jj) = col2im(feval(fun,x,params{:}),nhood,block,'distinct');
end
end
b = b(1:ma,1:na);
end
else
error(message('images:colfilt:unknownBlockType', deblank(kind)))
end
function [f,msg] = fcnchk(fun,varargin)
%FCNCHK Check FUNFUN function argument.
%
% FCNCHK will not accept string expressions for FUN in a future
% release. Use anonymous functions for FUN instead.
%
% FCNCHK(FUN,...) returns an inline object based on FUN if FUN
% is a string containing parentheses, variables, and math
% operators. FCNCHK simply returns FUN if FUN is a function handle,
% or a MATLAB object with an feval method (such as an inline object).
% If FUN is a string name of a function (e.g. 'sin'), FCNCHK returns a
% function handle to that function.
%
% FCNCHK is a helper function for FMINBND, FMINSEARCH, FZERO, etc. so they
% can compute with string expressions in addition to functions.
%
% FCNCHK(FUN,...,'vectorized') processes the string (e.g., replacing
% '*' with '.*') to produce a vectorized function.
%
% When FUN contains an expression then FCNCHK(FUN,...) is the same as
% INLINE(FUN,...) except that the optional trailing argument 'vectorized'
% can be used to produce a vectorized function.
%
% [F,ERR] = FCNCHK(...) returns a struct array ERR. This struct is empty
% if F was constructed successfully. ERR can be used with ERROR to throw
% an appropriate error message if F was not constructed successfully.
%
% See also ERROR, INLINE, @, FUNCTION_HANDLE.
% Copyright 1984-2012 The MathWorks, Inc.
msgident = '';
nin = nargin;
if nin > 1 && strcmp(varargin{end},'vectorized')
vectorizing = true;
nin = nin - 1;
else
vectorizing = false;
end
if ischar(fun)
fun = strtrim_local_function(fun);
% Check for non-alphanumeric characters that must be part of an
% expression.
if isempty(fun),
f = inline('[]');
elseif ~vectorizing && isidentifier_local_function(fun)
f = str2func(fun); % Must be a function name only
% Note that we avoid collision of f = str2func(fun) with any local
% function named fun, by uglifying the local function's name
if isequal('x',fun)
warning(message('MATLAB:fcnchk:AmbiguousX'));
end
else
if vectorizing
f = inline(vectorize(fun),varargin{1:nin-1});
var = argnames(f);
f = inline([formula(f) '.*ones(size(' var{1} '))'],var{1:end});
else
f = inline(fun,varargin{1:nin-1});
end
end
elseif isa(fun,'function_handle')
f = fun;
% is it a MATLAB object with a feval method?
elseif isobject(fun)
% delay the methods call unless we know it is an object to avoid
% runtime error for compiler
[meths,cellInfo] = methods(class(fun),'-full');
if ~isempty(cellInfo) % if fun is a MATLAB object
meths = cellInfo(:,3); % get methods names from cell array
end
if any(strmatch('feval',meths))
if vectorizing && any(strmatch('vectorize',meths))
f = vectorize(fun);
else
f = fun;
end
else % no feval method
f = '';
msgident = 'MATLAB:fcnchk:objectMissingFevalMethod';
end
else
f = '';
msgident = 'MATLAB:fcnchk:invalidFunctionSpecifier';
end
% If no errors and nothing to report then we are done.
if nargout < 2 && isempty(msgident)
return
end
% compute MSG
if isempty(msgident)
msg.message = '';
msg.identifier = '';
msg = msg(zeros(0,1)); % make sure msg is the right dimension
else
msg.identifier = msgident;
msg.message = getString(message(msg.identifier));
end
if nargout < 2
if ~isempty(msg)
error(message(msg.identifier));
end
end
%------------------------------------------
function s1 = strtrim_local_function(s)
%STRTRIM_LOCAL_FUNCTION Trim spaces from string.
% Note that we avoid collision with line 45: f = str2func('strtrim')
% by uglifying the local function's name
if isempty(s)
s1 = s;
else
% remove leading and trailing blanks (including nulls)
c = find(s ~= ' ' & s ~= 0);
s1 = s(min(c):max(c));
end
%-------------------------------------------
function tf = isidentifier_local_function(str)
% Note that we avoid collision with line 45: f = str2func('isidentifier')
% by uglifying the local function's name
tf = false;
if ~isempty(str)
first = str(1);
if (isletter(first))
letters = isletter(str);
numerals = (48 <= str) & (str <= 57);
underscore = (95 == str);
tf = all(letters | numerals | underscore);
end
end
function b=im2col(varargin)
%IM2COL Rearrange image blocks into columns.
% B = IM2COL(A,[M N],'distinct') rearranges each distinct
% M-by-N block in the image A into a column of B. IM2COL pads A
% with zeros, if necessary, so its size is an integer multiple
% of M-by-N. If A = [A11 A12; A21 A22], where each Aij is
% M-by-N, then B = [A11(:) A21(:) A12(:) A22(:)].
%
% B = IM2COL(A,[M N],'sliding') converts each sliding M-by-N
% block of A into a column of B, with no zero padding. B has
% M*N rows and will contain as many columns as there are M-by-N
% neighborhoods in A. If the size of A is [MM NN], then the
% size of B is (M*N)-by-((MM-M+1)*(NN-N+1). Each column of B
% contains the neighborhoods of A reshaped as NHOOD(:), where
% NHOOD is a matrix containing an M-by-N neighborhood of
% A. IM2COL orders the columns of B so that they can be
% reshaped to form a matrix in the normal way. For example,
% suppose you use a function, such as SUM(B), that returns a
% scalar for each column of B. You can directly store the
% result in a matrix of size (MM-M+1)-by-(NN-N+1) using these
% calls:
%
% B = im2col(A,[M N],'sliding');
% C = reshape(sum(B),MM-M+1,NN-N+1);
%
% B = IM2COL(A,[M N]) uses the default block type of
% 'sliding'.
%
% B = IM2COL(A,'indexed',...) processes A as an indexed image,
% padding with zeros if the class of A is uint8 or uint16, or
% ones if the class of A is double.
%
% Class Support
% -------------
% The input image A can be numeric or logical. The output matrix
% B is of the same class as the input image.
%
% Example
% -------
% Calculate the local mean using a [2 2] neighborhood with zero padding.
%
% A = reshape(linspace(0,1,16),[4 4])'
% B = im2col(A,[2 2])
% M = mean(B)
% newA = col2im(M,[1 1],[3 3])
%
% See also BLOCKPROC, COL2IM, COLFILT, NLFILTER.
% Copyright 1993-2012 The MathWorks, Inc.
[a, block, kind, padval] = parse_inputs(varargin{:});
if strcmp(kind, 'distinct')
% Pad A if size(A) is not divisible by block.
[m,n] = size(a);
mpad = rem(m,block(1)); if mpad>0, mpad = block(1)-mpad; end
npad = rem(n,block(2)); if npad>0, npad = block(2)-npad; end
aa = mkconstarray(class(a), padval, [m+mpad n+npad]);
aa(1:m,1:n) = a;
[m,n] = size(aa);
mblocks = m/block(1);
nblocks = n/block(2);
b = mkconstarray(class(a), 0, [prod(block) mblocks*nblocks]);
x = mkconstarray(class(a), 0, [prod(block) 1]);
rows = 1:block(1); cols = 1:block(2);
for i=0:mblocks-1,
for j=0:nblocks-1,
x(:) = aa(i*block(1)+rows,j*block(2)+cols);
b(:,i+j*mblocks+1) = x;
end
end
elseif strcmp(kind,'sliding')