gld_momentfit.m

这是盲信号的代码 都已经通过编译了 做这方面的同仁可以参考一下 我觉得蛮惯用的

function lambda = gld_momentfit(alpha3,alpha4,alpha1,alpha2)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% GLD_MOMENTFIT - Estimates the parameters of the Generalized 
% Lambda Distribution using the first four moments (the method of moments). 
%
% Input
%   alpha3 = sample 3rd moment
%   alpha4 = sample 4th moment
%   alpha1 = sample mean       (optional)
%   alpha2 = sample variance   (optional)
%  
% Output
%   lambda = the distribution parameters vector  
%            [lambda(1),lambda(2),lambda(3),lambda(4)]
%
% Parameters are linearly interpolated between the values 
% provided by the GLD_moment_table. The table is adapted from
%   E.J. Dudewicz, Z.A. Karian: "The Extended Generalized Lambda
%   Distribution (EGLD) System for Fitting Distributions to Data with 
%   Moments, II: Tables", American J. of Math. and Management
%   Sciences, 1996, Vol. 16, NOS. 3& 4, 271--332
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Jan Eriksson, Juha Karvanen, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.
%
% Last modified: 5.9.2000
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Constant values defining the size of the GLD table. Should be
% changed only if the table range is changed.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
table_min_base  = 1.8;
table_min_slope = 1.75;
table_max_base  = 14.7;
table_max_slope = 1.8;
  
if nargin<2,  disp('Usege: gld_momentfit(alpha3,alpha4) or');
              disp('gld_momentfit(alpha3,alpha4,alpha1,alpha2)');
              error('Not enough input arguments.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% If alpha1 and alpha2  are not given, the default values alpha1=0
% and alpha2=1 are substituted.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if nargin==2, alpha1=0;
              alpha2=1;
end
global GLD_moment_table;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check if the moment values are covered by the table.
% If not, the closest point is used.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (table_min_base+table_min_slope*alpha3^2>alpha4),
              disp(' ');
              disp(['Warning: The moment values, alpha3=',num2str(alpha3),...
		    ', alpha4=',num2str(alpha4),',']);
	      disp('         are not covered by the GLD distribution.');
	      disp('         The closest point is used.');
end
if (table_max_base+table_max_slope*alpha3^2<alpha4),
              disp(' ');
              disp(['Warning: The moment values, alpha3=',num2str(alpha3),...
		    ', alpha4=',num2str(alpha4),',']);
	      disp('         are not covered by the GLD tables.');
	      disp('         The closest point is used.');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The closest points corresponding  to <= alpha3, >= alpha3,
% <=alpha4, and >=alpha4 are searced from the table. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ta3=GLD_moment_table(:,1);

a3max=max(ta3);
a3abs=min(abs(alpha3),a3max);

a3high=min(ta3(find(ta3>=a3abs)));
a3hplace=find(ta3==a3high);
a3low=max(ta3(find(ta3<=a3abs)));
a3lplace=find(ta3==a3low);

ta4h=GLD_moment_table(a3hplace,2:6);

a4hlow= max(ta4h(find(ta4h(:,1)<=alpha4),1));
if isempty(a4hlow)
  a4hlow=min(ta4h(:,1));
end
a4hhigh=min(ta4h(find(ta4h(:,1)>=alpha4),1));
if isempty(a4hhigh)
  a4hhigh=max(ta4h(:,1));
end

lambda_a4hlow=ta4h(find(a4hlow==ta4h(:,1)),2:5);
lambda_a4hhigh=ta4h(find(a4hhigh==ta4h(:,1)),2:5);
if a4hlow==a4hhigh
     a4hlambda=lambda_a4hlow;
else
  a4hlambda=((a4hhigh-alpha4)*lambda_a4hlow+(alpha4-a4hlow)*lambda_a4hhigh)/(a4hhigh-a4hlow);
end

ta4l=GLD_moment_table(a3lplace,2:6);

a4llow= max(ta4l(find(ta4l(:,1)<=alpha4),1));
if isempty(a4llow)
  a4llow=min(ta4l(:,1));
end
a4lhigh=min(ta4l(find(ta4l(:,1)>=alpha4),1));
if isempty(a4lhigh)
  a4lhigh=max(ta4l(:,1));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The parameters vector is obtained by linearily interpolating between
% the lambda values of the points obtained above. This vector is
% for the standarized data (mean=0, variance=1).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lambda_a4llow=ta4l(find(a4llow==ta4l(:,1)),2:5);
lambda_a4lhigh=ta4l(find(a4lhigh==ta4l(:,1)),2:5);
if a4llow==a4lhigh
     a4llambda=lambda_a4llow;
else
  a4llambda=((a4lhigh-alpha4)*lambda_a4llow+(alpha4-a4llow)*lambda_a4lhigh)/(a4lhigh-a4llow);
end

if a3low==a3high 
  lambda=a4llambda;
else
  lambda=((a3high-a3abs)*a4llambda+(a3abs-a3low)*a4hlambda)/ ...
	 (a3high-a3low);
end
  
if sign(alpha3)==-1
  lambda=[-lambda(1) lambda(2) lambda(4) lambda(3)];
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now the mean and the variance are taken into account.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lambda=[lambda(1)*alpha2+alpha1,lambda(2)/alpha2,lambda(3),lambda(4)];

% The end of the function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%