%
disp(' ');
disp(' white_kurtosis.m  ver 1.2  April 17, 2011 ');
disp(' ');
disp(' by Tom Irvine  Email: tomirvine@aol.com');
disp(' ');
disp(' This program synthesizes a white noise time history ');
disp(' where the user specifies the kurtosis ');
disp(' ');
%
clear w;
clear accel;
clear xarea;
clear area;
clear xx;
clear yy;
clear ar;
clear br;
clear store;
clear table;
clear x;
clear X;
clear Z;
%
clear xx;
clear yy;
%
krec=0;
srec=0;
%
max_error=1.0e+90;
%
nv=1;
%
disp(' ');
disp(' Enter the standard deviation ');
sdd=input(' ');
%
disp(' ');
disp(' Enter the kurtosis ');
kde=input(' ');
%
disp(' ');
disp(' Enter sample rate (samples/sec)');
sr=input(' ');
dt=1/sr;
%
disp(' Enter the duration (sec) ');
dur=input(' ');
%
disp(' Enter the number of trials ');
trials=input(' ');
%
nt=round(dur/dt);
%
accel=zeros(nt,2);
%
np=1001;
%
for(i=1:nt)    
%
  tt=dt*(i-1);
  accel(i,1)=tt;
end
%
xx=linspace(accel(1,1),accel(nt,1),10*nt);
x=zeros(np,1);
x=linspace(-10,10,np);
x=x';
dx=x(2)-x(1);
%
kmax=1.0e+80;
%
disp(' ');
disp(' Trial  Skewness   Kurtosis    nex    kex ');
%
progressbar % Create figure and set starting time
num=trials;
%
for(ijk=1:num)
    progressbar(ijk/num);
ar=zeros(np,1);
br=zeros(np,1);
br=linspace(0,1,np);
pr=zeros(np,1);
%
mu=0;
sigma=1;
den=2*sigma^2;
%
nex=10*rand();
kex=1.5*rand();
%
if(ijk>5)
    qqq=rand();
    if(qqq>0.5 & qqq<0.75 )
        nex=nex_rec*(0.04*rand()+0.98);
        kex=kex_rec*(0.04*rand()+0.98);   
    end
%
    if(qqq>=0.75)
        nex=nex_rec;
        kex=kex_rec;   
    end
end
%
if(ijk==1)
    nex=1;
    kex=1;
end
%
for(i=1:np)
    if(abs(x(i))>0)
        sss=x(i)/abs(x(i));
    else
       sss=0;
    end
    xy=(abs(x(i))^kex)*sss;
    arg=-(xy)^2/den;   % mu=0
    pr(i)=exp(arg);
    pr(i)=pr(i)^nex;
end
pr=pr/sqrt(2*pi*sigma^2);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
    ar=zeros(np,1);
    ar(1)=pr(1)*dx/2;
%
    for(i=2:(np-1))
        ar(i)=ar(i-1)+pr(i)*dx;
    end
    ar(np)=ar(np-1)+pr(np)*dx/2;
%
    ar=ar/ar(np);
%    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
xarea=[x ar];
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
nx=max(size(xarea));
%
for(i=1:nt)    
%
a=1.0*rand();
%
    for(j=1:nx-1)
        if( a>=xarea(j,2) & a<=xarea(j+1,2))
            len=xarea(j+1,1)-xarea(j,1);
            xa=(a-xarea(j,2));
            c2=xa/len;
            c1=1-c2;
            b=c1*xarea(j,1)+c2*xarea(j+1,1);
            break;
        end
    end
%
    accel(i,2)=b;
%
end
%
%% disp(' ');
%% disp(' Spline fit ');
%
accel(:,2)=accel(:,2)-mean(accel(:,2));
yy = spline(accel(:,1),accel(:,2),xx);
%
yy=yy*sdd/std(yy);
%
clear amp;
%
    amp=yy-mean(yy);
    n = length(amp);
%
    mu=mean(amp);
    sd=std(amp);
    kt=0.;
    sk=0.;
%
    for i=1:n
        kt=kt+(amp(i)-mu)^4;
        sk=sk+(amp(i)-mu)^3;
    end      
    kt=kt/(n*sd^4);
    sk=sk/(n*sd^3);
%
    store(ijk,1)=sk;
    store(ijk,2)=kt;
%
   sde=0;
   error=abs(kt-kde) +abs(sk-sde);
%
    if(error < max_error)
%
        max_error=error;
        clear accel_best;
        accel_best=[xx' yy'];
        kmax=abs(kt-kde);
        krec=kt;
        srec=sk;
        nex_rec=nex;
        kex_rec=kex;
%
        out1=sprintf('   %d  %8.4g  %8.4g   %8.4g  %8.4g ',ijk,sk,kt,nex_rec,kex_rec);
        disp(out1);
%
    end
    table(nv,:)=[kt,nex,kex];
    nv=nv+1;
%
end
%
    disp(' ');
    disp(' The time history is:  accel_best ');
%
fignum=1;
figure(1);
fignum=fignum+1;
mu=mean(accel_best(:,2));
sd=std(accel_best(:,2));
plot(accel_best(:,1),accel_best(:,2));
if(mu<sd/500)
    mu=0;
end
out1=sprintf(' Time History  mean=%6.3g std dev=%6.3g skew=%6.3g kurtosis=%6.3g',mu,sd,srec,krec);
title(out1);
ylabel('Amplitude');
xlabel('Time(sec)');
grid on;
%
clear amp;
amp=accel_best(:,2);
%
    disp(' ')
    disp(' Plot histogram? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
        while(1)
            disp(' ');
            nbar=input(' Enter the number of bars ');
            fignum=fignum+1;
            figure(fignum);    
            hist(amp,nbar)
            title(' Histogram ');
            ylabel(' Counts ');
            xlabel('Amplitude');
            disp(' ');
            disp(' Re-plot with different number of bars? ');
            disp(' 1=yes  2=no ');
            irp=input(' ');
            if(irp==2)
                break;
            end
        end
%
        disp(' ')
        disp(' Plot cumulative histogram? ')
            cum_choice = input(' 1=yes 2=no '); 
% 
            if(cum_choice==1)
                [nnn,xout]=hist(amp,nbar);
                c_elements = cumsum(nnn);
                fignum=fignum+1;
                figure(fignum);
                bar(xout,c_elements);
                title(' Cumulative Histogram ');
                ylabel(' Counts ');
                xlabel('Amplitude');
                grid on;
            end
%
    end