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 sourcesif 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;endreturn%%%%%%%% Test %%%%%%%%%%%%%figure; c=-10; x=[-1.1:0.1:1.1]';[f,df, d2f]=psi(x,c);plot(x,[f,df, d2f]);grideps=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