tomwalters@0: 
tomwalters@0: function generate_kurtosis_noise(sig,kurtois
tomwalters@0: 
tomwalters@0: v=get(glob.f);
tomwalters@0: % variance_wanted=get(params,'Variance');
tomwalters@0: % skwednes=get(params,'Skewedness');
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % change the kurtosis input to something meaningful
tomwalters@0: % these values are measured for different inputs of x
tomwalters@0: % if raw value for "Kurtosis" is inputted as:
tomwalters@0: input_kurt=[1 0.95 0.9 0.85 .8 .7 .6 .5 .4 .3 .25 .2 .15 .125 0.1 ];
tomwalters@0: % then the Kurtosis will be:
tomwalters@0: output_kurt=[-.05 0.04 0.13 0.26 .4 .7 1.12 1.72 2.58 4 5.11 6.87 9.6 11.8 15];
tomwalters@0: 
tomwalters@0: % to input the real Kurtosis, we need to turn these values around:
tomwalters@0: % and find the "input" value for a given Kurtosis:
tomwalters@0: kurtosis=interp1(output_kurt,input_kurt,kurtosis_wanted);
tomwalters@0: 
tomwalters@0: % % and the same for the variance. changing the variance does not change the
tomwalters@0: % % kurtosis
tomwalters@0: % input_var=0:15;
tomwalters@0: % output_var=[0.97 .68 .53 .43 .35 .3 .27 .24 .22 .2 .18 .17 .15 .14 .13 .12];
tomwalters@0: % variance=interp1(output_var,input_var,variance_wanted);
tomwalters@0: 
tomwalters@0: variance=1;
tomwalters@0: 
tomwalters@0: grafix=0;
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0:     figure(3)
tomwalters@0:     clf
tomwalters@0:     hold on
tomwalters@0: end
tomwalters@0: 
tomwalters@0: nr_points=2000;
tomwalters@0: 
tomwalters@0: stepx=10/nr_points;
tomwalters@0: % generate the pdf 
tomwalters@0: x=-5:stepx:5;
tomwalters@0: 
tomwalters@0: p_pdf=pdf('Normal',x,0,1);
tomwalters@0: p_pdf=p_pdf/max(p_pdf);
tomwalters@0: 
tomwalters@0: % p_cdf=cumsum(p_pdf);
tomwalters@0: % p_cdf=p_cdf/max(p_cdf);
tomwalters@0: % if grafix
tomwalters@0: %      plot(x,p_pdf,'b.');hold on
tomwalters@0: % end
tomwalters@0: 
tomwalters@0: % scale along x to change the variance:
tomwalters@0: y=p_pdf;
tomwalters@0: 
tomwalters@0: 
tomwalters@0: new_x=x*variance;
tomwalters@0: ny=interp1(new_x,y,x);
tomwalters@0: ny(isnan(ny))=0;
tomwalters@0: 
tomwalters@0: glob.pdf_x=x; % save for plotting
tomwalters@0: glob.pdf_y=ny;
tomwalters@0: 
tomwalters@0: % % and the change in y to modify the Kurtosis
tomwalters@0: % 
tomwalters@0: % generate random numbers by reversing the cdf
tomwalters@0: new_cdf=cumsum(ny); %calculate the pdf:
tomwalters@0: new_cdf=new_cdf/max(new_cdf); % and normalise
tomwalters@0: 
tomwalters@0: for i=1:nr_points/2-1
tomwalters@0: %     dx=abs(x(nr_points/2-i));
tomwalters@0:     dy=abs(0.5-new_cdf(i));
tomwalters@0:     ndy=power((2*dy),kurtosis);
tomwalters@0:     new_cdf(i)=0.5-ndy/2;
tomwalters@0:     new_cdf(nr_points-i)=0.5+ndy/2;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: if grafix
tomwalters@0:     
tomwalters@0:     old_cdf=cumsum(y);
tomwalters@0:     old_cdf=old_cdf/max(old_cdf); % and normalise
tomwalters@0:     
tomwalters@0:     plot(x,old_cdf,'-b.');
tomwalters@0:     plot(x,new_cdf,'-r.'); hold on
tomwalters@0:     legend({'cumulative normal distribution';'modfied cdf'},2);
tomwalters@0: end
tomwalters@0: 
tomwalters@0: 
tomwalters@0: serachy=1/nr_points:1/nr_points:1;
tomwalters@0: rev_cdf=zeros(1,nr_points-1);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: 
tomwalters@0: for i=1:nr_points-1
tomwalters@0:     rev_cdf(i)=find(new_cdf>serachy(i),1,'first');
tomwalters@0: end
tomwalters@0: 
tomwalters@0: % clf
tomwalters@0: % plot(rev_cdf,'r.'); 
tomwalters@0: 
tomwalters@0: r=round(rand(length(v),1)*(nr_points-1)+0.5);
tomwalters@0: 
tomwalters@0: rev_p=rev_cdf(r);
tomwalters@0: 
tomwalters@0: % clf
tomwalters@0: % hist(rev_p,100)
tomwalters@0: 
tomwalters@0: 
tomwalters@0: % and normalise to the right output range
tomwalters@0: rev_p=(rev_p-nr_points/2)/nr_points*10;
tomwalters@0: 
tomwalters@0: % return
tomwalters@0: % v=random('norm',0,sqrt(variance),length(v),1);
tomwalters@0: 
tomwalters@0: 
tomwalters@0: glob.f=set(glob.f,rev_p');
tomwalters@0: glob.f=scaletorms(glob.f,sqrt(get(params,('rms-value'))));
tomwalters@0: 
tomwalters@0: assignin('base','f',glob.f);
tomwalters@0: 
tomwalters@0: % f=scaletomaxvalue(f,1);