%ExampleStat.m sample program
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                               %
%    Extraction of the weibull parameters of the rectangular    %
%         components of the E-field over M observations         %
%                                                               %
%        by E. Amador (emmanuel.amador@insa-rennes.fr)          %
%                         IETR/DGA                              %
%                                                               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear

global Lt c R N POS

tic
c = 299792458;%
Lt = 1e-6; %Time-window length in seconds WARNING: for unbiased results, one should simulate at least 3*tau (see model companion)
nbre_elements = 1; %number of radiating elements

%loading the Position matrix from the image generator
filename = sprintf('%delem_%dns1s8.mat',nbre_elements,round(Lt/(1e-9)));
load(filename)

f_l = 300e6; %lowest frequency
f0 = 2e9; %max frequency

N = round(3*Lt*f0) %number of points for the chosen time-window (Lt)
t = 0:Lt/(N-1):Lt; %time scale


%Loss coefficient
R = 0.998;

%Dimensions of the chamber
l = 8.7;
p = 3.7;
h = 2.9;


M=150; % Number of reception points equiv. number of stirrer positions

%A distance of at least lambda/2 at the lowest frequency f_l
%should be respected between each reception points and with the walls

lambda=c/f_l;
X=[(l-lambda)*rand+lambda/2];
Y=[(p-lambda)*rand+lambda/2];
Z=[(h-lambda)*rand+lambda/2];

while (length(X)<M)
    X_1=(l-lambda)*rand+lambda/2;
    Y_1=(p-lambda)*rand+lambda/2;
    Z_1=(h-lambda)*rand+lambda/2;
    D=sqrt((X_1-X).^2+(Y_1-Y).^2+(Z_1-Z).^2);
    if min(D)>lambda/2
        X=[X;X_1];
        Y=[Y;Y_1];
        Z=[Z;Z_1];
    end
end

for i=1:M
    [Sx,Sy,Sz] = CIR(X(i),Y(i),Z(i));
    
    Fs = N/Lt; %sampling frequency
    T = 1/Fs;                     % sample length
    L = N;                  %Number of points
    tt = (0:L-1)*T;                % time...
    
    NFFT = 2^nextpow2(N); % Next power of 2 from length of y
    Yx = fft((Sx),NFFT);
    Yy = fft((Sy),NFFT);
    Yz = fft((Sz),NFFT);
    f = Fs/2*linspace(0,1,NFFT/2);
    
    % Three axis frequency response
    FFTx(i,:) = Yx(1:NFFT/2);
    FFTy(i,:) = Yy(1:NFFT/2);
    FFTz(i,:) = Yz(1:NFFT/2);
    freq = f;
    
    disp([num2str(round((i/M)*1000)/10),' %'])
    
end

%extraction of the weibull parameters
%definition used of the Weibull:
%Weibull pdf is give by f(x)=a*x^(b-1)*exp(-a*x^b/b)
%the statistical toolbox is needed if not, a mean square etimation of the
%parameters is possible and given ine the litterature

disp('Weibull parameters extraction')
for k=1:length(freq)
   paramx = wblfit(abs(FFTx(:,k))/mean(abs(FFTx(:,k))));
   paramy = wblfit(abs(FFTy(:,k))/mean(abs(FFTy(:,k))));
   paramz = wblfit(abs(FFTz(:,k))/mean(abs(FFTz(:,k))));
   
   bx(k)=paramx(2);
   ax(k)=1/(paramx(1).^bx(k));
   by(k)=paramy(2);
   ay(k)=1/(paramy(1).^by(k));
   bz(k)=paramz(2);
   az(k)=1/(paramz(1).^bz(k));
    
end

%filtering of the results over p frequencies
p=10;

axf=filter(ones(p,1)/p,1,ax);
ayf=filter(ones(p,1)/p,1,ay);
azf=filter(ones(p,1)/p,1,az);

bxf=filter(ones(p,1)/p,1,bx);
byf=filter(ones(p,1)/p,1,by);
bzf=filter(ones(p,1)/p,1,bz);


figure(1)
subplot(2,1,1)
plot(freq/1e6,axf,freq/1e6,ayf,freq/1e6,azf)
title('Scale parameter $a$','Interpreter','Latex')
xlabel('Frequency [MHz]','Interpreter','Latex')
grid on
h=legend('$E_x$','$E_y$','$E_z$');
set(h,'Interpreter','Latex')
grid on
xlim([100 1000])
ylim([.7 .9])

subplot(2,1,2)
plot(freq/1e6,bxf,freq/1e6,byf,freq/1e6,bzf)
title('Shape parameter $b$','Interpreter','Latex')
xlabel('Frequency [MHz]','Interpreter','Latex')
h2=legend('$E_x$','$E_y$','$E_z$');
set(h2,'Interpreter','Latex')
grid on
xlim([100 1000])
ylim([1.3 2.2])
toc