%
disp(' ');
disp(' witch.m  ver 1.1  March 23, 2011 ');
disp(' ');
disp(' by Tom Irvine ');
%
disp(' ');
disp(' This program generates a non-Gaussian white noise time history where ')
disp(' the corresponding histogram has the shape of a witches hat ');
disp(' ');
%
disp(' Select type ');
disp(' 1=  skewness =0     kurtosis ~4.7  ');
disp(' 2=  skewness ~1.4   kurtosis ~5.2  ');
%
ichoice=input(' ');
% 
clear w;
clear accel;
clear xarea;
clear area;
clear xx;
clear yy;
%
if(ichoice==1)
    w=[-10	0;
	   -3	0.1;
	    0	1;
	    3	0.1;
	   10	0];
else   
%
    w=[0	0;
	   2	1;
	   3	10;
	   20	1;
	   50	0];
%
end
%
area(1)= 0.5*abs(w(1,1)-w(2,1))*w(2,2);
area(2)=     abs(w(3,1)-w(2,1))*w(2,2); 
area(3)= 0.5*abs(w(3,1)-w(2,1))*(w(3,2)-w(2,2));
area(4)= 0.5*abs(w(4,1)-w(3,1))*(w(3,2)-w(2,2));
area(5)=     abs(w(4,1)-w(3,1))*w(4,2); 
area(6)= 0.5*abs(w(5,1)-w(4,1))*w(4,2);
%
total_area=sum(area);
%
w(:,2)=w(:,2)/total_area;
%
disp(' ');
sr=input(' Enter initial sample rate (samples/sec) ');
disp(' ');
%
dt=1/sr;
%
dur=input(' Enter duration (sec) ');
%
sdd=input(' Enter standard deviation ');
%
nt=round(dur/dt);
%
err=1.0e+90;
xmean=0;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
n=1001;
step=(w(5,1)-w(1,1))/(n-1);
%
clear xarea;
xarea=zeros(n,2);
xprevious1=0;
xprevious2=0;
xprevious3=0;
%
         x1=w(1,1);
         x2=w(2,1);
         y1=w(1,2);
         dy=w(2,2)-w(1,2);
         dx=w(2,1)-w(1,1);
         m1=dy/dx;
         b1=y1-m1*x1;
xprevious1=(m1*x2^2/2)+b1*x2 - ((m1*x1^2/2)+b1*x1);
         x1=w(2,1);
         x2=w(3,1);
         y1=w(2,2);
         dy=w(3,2)-w(2,2);
         dx=w(3,1)-w(2,1);
         m2=dy/dx;
         b2=y1-m2*x1;     
xprevious2=(m2*x2^2/2)+b2*x2 - ((m2*x1^2/2)+b2*x1);
         x1=w(3,1);
         x2=w(4,1);
         y1=w(3,2);
         dy=w(4,2)-w(3,2);
         dx=w(4,1)-w(3,1);
         m3=dy/dx;
         b3=y1-m3*x1;      
xprevious3=(m3*x2^2/2)+b3*x2 - ((m3*x1^2/2)+b3*x1);
         x1=w(4,1);
         x2=w(5,1);
         y1=w(4,2);
         dy=w(5,2)-w(4,2);
         dx=w(5,1)-w(4,1);
         m4=dy/dx;
         b4=y1-m4*x1;
%
for(i=1:n)
%    
    x=(i-1)*step+w(1,1);
%    
    if(x>=w(1,1) & x<w(2,1))
         x1=w(1,1);
         xarea(i,2)=(m1*x^2/2)+b1*x - ((m1*x1^2/2)+b1*x1);
    end
    if(x>=w(2,1) & x<w(3,1))
         x1=w(2,1);   
         xarea(i,2)=(m2*x^2/2)+b2*x - ((m2*x1^2/2)+b2*x1) +xprevious1;   
    end
    if(x>=w(3,1) & x<w(4,1))
         x1=w(3,1);     
         xarea(i,2)=(m3*x^2/2)+b3*x - ((m3*x1^2/2)+b3*x1) +xprevious1+xprevious2;   
    end
    if(x>=w(4,1) & x<=w(5,1))
         x1=w(4,1);   
         xarea(i,2)=(m4*x^2/2)+b4*x - ((m4*x1^2/2)+b4*x1) +xprevious1+xprevious2+xprevious3;    
    end    
%
    xarea(i,1)=x;
%
    if( abs(xarea(i,2)-0.5)<err)
        err=abs(xarea(i,2)-0.5);
        xmean=xarea(i,1);
    end
%
end
%
if(ichoice ==2)
    xarea(:,1)=xarea(:,1)-xmean;
end    
%
figure(1);
plot(w(:,1),w(:,2));
title(' Probability Density Function');
xlabel('Amplitude');
ylabel('Probability');
grid on;
%
figure(2);
plot(xarea(:,1),xarea(:,2));
title('Cumulative Probability');
ylabel('Area');
xlabel('Amplitude');
grid on;
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
clear accel;
accel=zeros(nt,2);
%
nx=max(size(xarea));
%
progressbar % Create figure and set starting time
%
for(i=1:nt)
%
progressbar(i/nt);    
%
tt=dt*(i-1);
%
a=1.0*rand();
%
    for(j=1:nx-1)
        if( a>=xarea(j,2) & a<=xarea(j+1,2))
            len=dx;
            x=(a-xarea(j,2));
            c2=x/len;
            c1=1-c2;
            b=c1*xarea(j,1)+c2*xarea(j+1,1);
            break;
        end
    end
%
    accel(i,1)=tt;
    accel(i,2)=b;
%
end
%
clear xx;
clear yy;
%
disp(' ');
disp(' Spline fit ');
%
xx=linspace(accel(1,1),accel(nt,1),10*nt);
%
%% figure(2);
%% plot(accel(:,1),accel(:,2));
%
accel(:,2)=accel(:,2)-mean(accel(:,2));
yy = spline(accel(:,1),accel(:,2),xx);
%
yy=yy*sdd/std(yy);
%
clear amp;
clear tim;
%
    amp=yy;
    tim=xx;
    n = length(amp);
%disp(' mean values ')
    mu=mean(amp);
    sd=std(amp);
    mx=max(amp);
    mi=min(amp);
    rms=sqrt(sd^2+mu^2);
    kt=0.;
    sk=0.;
%1
    for i=1:n
        if( amp(i)==mx)
            tt_max=tim(i);
        end
        if( amp(i)==mi)
            tt_min=tim(i);
        end
        kt=kt+(amp(i)-mu)^4;
        sk=sk+(amp(i)-mu)^3;
    end      
    kt=kt/(n*sd^4);
    sk=sk/(n*sd^3);
%1
    disp(' ')
    disp(' Time Stats ')
    tmx=max(tim);
    tmi=min(tim);
    duration=(tmx-tmi);
% 
    out3 = sprintf('   start  = %g sec    end = %g sec   duration = %g sec  ',tmi,tmx,duration);
    disp(out3)
%
    dt=duration/(n-1);
    out5 = sprintf('\n number of samples = %d  ',n);
    disp(out5)
%
    clear t;
    t=tim;
    disp(' ')
    disp(' Time Step ');
    clear difft;
    difft=diff(t);
    dtmin=min(difft);
    dtmax=max(difft);
%
    out4 = sprintf(' dtmin  = %8.4g sec  ',dtmin);
    out5 = sprintf(' dt     = %8.4g sec  ',dt);
    out6 = sprintf(' dtmax  = %8.4g sec  ',dtmax);
    disp(out4)
    disp(out5)
    disp(out6)
%
    disp(' ')
    disp(' Sample Rate ');
    out4 = sprintf(' srmin  = %8.4g samples/sec  ',1/dtmax);
    out5 = sprintf(' sr     = %8.4g samples/sec  ',1/dt);
    out6 = sprintf(' srmax  = %8.4g samples/sec  \n',1/dtmin);
    disp(out4)
    disp(out5)
    disp(out6)
    clear t;
%
    disp(' ')
    disp(' Amplitude Stats ')
     out0 = sprintf(' number of points = %d ',n);
    out1 = sprintf(' mean = %8.4g    std = %8.4g    rms = %8.4g ',mu,sd,rms);
    out2a = sprintf(' max  = %9.4g  at  = %8.4g sec             ',mx,tt_max);
    out2b = sprintf(' min  = %9.4g  at  = %8.4g sec            \n',mi,tt_min);
    disp(out0);
    disp(out1);
    disp(out2a);
    disp(out2b);
%
    amax=abs(mx);
    amin=abs(mi);
    if(amax<amin)
        amax=amin;
    end
    crest=amax/rms;
    out4 = sprintf(' crest factor = %8.4g ',crest);
    out5 = sprintf(' skewness = %8.4g  kurtosis  = %8.4g ',sk,kt);
    disp(out4);
    disp(out5);
%
clear accel;
accel=[xx' yy'];
disp(' ');
disp(' The time history is:  accel ');
%
figure(3);
plot(xx,yy);
if(mu<sd/500)
    mu=0;
end
out1=sprintf('Time History  mean=%4.2g  std dev=%4.3g skew=%4.2g kurtosis=%4.2g',mu,sd,sk,kt);
title(out1);
ylabel('Amplitude');
xlabel('Time(sec)');
grid on;
%
    disp(' ')
    disp(' Plot histogram? ')
    choice = input(' 1=yes 2=no ');
%
    if(choice == 1)
        while(1)
            disp(' ');
            nbar=input(' Enter the number of bars ');
            figure(5);    
            hist(yy,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
    end