% Script file for computing and plotting Symmetric-Alpha-Stable distributions
%
%-----------------------------------------------------------
% Copyright 2000, by Dimitris G. Manolakis, Vinay K. Ingle,
% and Stephen M. Kogon.  For use with the book
% "Statistical and Adaptive Signal Processing"
% McGraw-Hill Higher Education.
%-----------------------------------------------------------


SCRUN = 'normalized'; 
SCRPOS = [0.002,0.03,0.48,0.44];
PAPUN = 'inches';
close all;
papos = [0,0,5,5];
Hf_1 = figure('units',SCRUN,'position',SCRPOS,...
   'paperunits',PAPUN,'paperposition',papos);

% (a) 1 < Alpha <= 2
alpha = 2; sigma = 1/sqrt(2); beta = 0;
xi = linspace(-50,50,2001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*sign(xi)*tan(pi*alpha/2)));
phi = [phi(1001:2001),zeros(1,2048-2001),conj(phi(1:1000))];
pdf = abs(ifft(phi));
pdf = [pdf(1025:2048),pdf(1:1025)];
x = [-1024:1024]*2*pi/100;

% Characteristic Function
xi = linspace(-5,5,1001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*sign(xi)*tan(pi*alpha/2)));
phi = abs(phi);
subplot('position',[0.05,0.75,0.4,0.2]);
plot(xi,phi); axis([-5,5,0,1.1]);
ylabel('\Phi(\xi)');%xlabel('\xi'); 
title('\alpha = 2, \sigma = 1, \beta = 0, \mu =0');
set(gca,'xtick',[-5,0,5],'ytick',[0,1],'fontsize',8);

% Density function
pdfmax = (sum(phi)*0.01)/(2*pi);
sf = pdfmax/max(pdf);
subplot('position',[0.58,0.75,0.4,0.2]);
%plot(x,pdf*sf); axis([-5,5,0,0.5]);
plot(x,pdf*sf); axis([-5,5,0,ceil(pdfmax*10)/10]);
ylabel('f(x)');%xlabel('x'); 
title('Gaussian pdf');
set(gca,'xtick',[-5,0,5],'ytick',[0,pdfmax],'fontsize',8);

% (b) Alpha = 1
alpha = 1; sigma = 1; beta = 0;
xi = linspace(-50,50,2001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*(2/pi)*sign(xi).*log(abs(xi)+eps)));
phi = [phi(1001:2001),zeros(1,2048-2001),conj(phi(1:1000))];
pdf = real(ifft(phi));
pdf = [pdf(1025:2048),pdf(1:1025)];
x = [-1024:1024]*2*pi/100;

% Characteristic Function
xi = linspace(-5,5,1001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*(2/pi)*sign(xi).*log(abs(xi)+eps)));
phi = abs(phi);
subplot('position',[0.05,0.43,0.4,0.2]);
plot(xi,phi); axis([-5,5,0,1.1]);
ylabel('\Phi(\xi)');%xlabel('\xi'); 
title('\alpha = 1, \sigma = 1, \beta = 0, \mu =0');
set(gca,'xtick',[-5,0,5],'ytick',[0,1],'fontsize',8);

% Density function
pdfmax = (sum(phi)*0.01)/(2*pi);
sf = pdfmax/max(pdf);
subplot('position',[0.58,0.43,0.4,0.2]);
%plot(x,pdf*sf); axis([-5,5,0,0.4]);
plot(x,pdf*sf); axis([-5,5,0,ceil(pdfmax*10)/10]);
ylabel('f(x)');%xlabel('x');
title('Cauchy pdf');
set(gca,'xtick',[-5,0,5],'ytick',[0,pdfmax],'fontsize',8);

% (c) 0 < Alpha < 1
alpha = 0.5; sigma = 2; beta = 0;
xi = linspace(-50,50,2001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*sign(xi)*tan(pi*alpha/2)));
phi = [phi(1001:2001),zeros(1,2048-2001),conj(phi(1:1000))];
pdf = real(ifft(phi));
pdf = [pdf(1025:2048),pdf(1:1025)];
x = [-1024:1024]*2*pi/100;

% Characteristic Function
xi = linspace(-5,5,1001);
phi = exp(-((abs(sigma*xi)).^alpha).*(1-j*beta*sign(xi)*tan(pi*alpha/2)));
phi = abs(phi);
subplot('position',[0.05,0.1,0.4,0.2]);
plot(xi,phi); axis([-5,5,0,1.1]);
ylabel('\Phi(\xi)');%xlabel('\xi'); 
title('\alpha = 0.5, \sigma = 1, \beta = 0, \mu =0');
set(gca,'xtick',[-5,0,5],'ytick',[0,1],'fontsize',8);

% Density function
pdfmax = (sum(phi)*0.01)/(2*pi);
sf = pdfmax/max(pdf);
subplot('position',[0.58,0.1,0.4,0.2]);
%plot(x,pdf*sf); axis([-5,5,0,0.5]);
plot(x,pdf*sf); axis([-5,5,0,ceil(pdfmax*10)/10]);
ylabel('f(x)'); xlabel('x');
title('Levy pdf');
set(gca,'xtick',[-5,0,5],'ytick',[0,pdfmax],'fontsize',8);