psi.m

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

function [f,df,d2f]=psi(t,c),

%                 c > 0 - use  smoothing of the absolute value function with parameter 'c'
%                              become more sharp and accurate  when c->0
%                                recommended for sparse/sepergaussian sources 
%                                (see the paper newt_ica_ica2003.pdf)
%
%                  c < 0 - use  the power function |t|^|c| (lam= -4 ... -20),  
%                                recommended for subgaussian sources


if c<0,          % the power function 
  c=-c;  
  f=t.^c;
  df=c.*t.^(c-1);
  d2f=c.*(c-1).*t.^(c-2);

else              %smoothing of the absolute value function 

  tt=(1./c).*t;
  u = abs(tt);
  f= c.*(u-log(1+u));
  df=tt./(1+u);
  d2f= (1/c)./ (1+u).^2;

end

return


%%%%%%%% Test %%%%%%%%%%%%%
figure; c=-10; x=[-1.1:0.1:1.1]';[f,df, d2f]=psi(x,c);plot(x,[f,df, d2f]);grid
eps=1e-7; [f1,df1, d2f1]=psi(x+eps,c);
dfnum=(f1-f)/eps;
d2fnum=(df1-df)/eps;
err_df=(df-dfnum)', err_d2f=(d2f-d2fnum)'
figure;plot(x,[d2f, d2fnum]);grid