--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/m-toolbox/test/test_matrix_MultiChannelNoise_2.m	Wed Nov 23 19:22:13 2011 +0100
@@ -0,0 +1,613 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% TEST matrix/MultiChannelNoise
+% 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% HISTORY:     23-04-2009 L Ferraioli
+%                 Creation
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%% VERSION
+
+'$Id: test_matrix_fromCSD.m,v 1.3 2009/04/23 16:03:55 luigi Exp $';
+
+%% Loading spectra
+
+fs = 10;
+f = [logspace(-6,log10(fs/2),2000)]';
+
+% currPath = cd;
+% cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
+% [TF,CSD] = mdc1_tf_models(plist('f',f,'fs',fs));
+% cd(currPath)
+
+% get noise psd shapes
+% parameters names
+parnames = {'DH','S11','S1D','SD1','SDD','Dh1','Dh2','w1','w2',...
+  'TH','Th1','Th2'};
+
+% nominal params values
+parvalues = {1,1,0,-1e-4,1,1,1,-1.3e-6,-2.0e-6,0.35,0.25,0.28};
+
+% Noise model
+Ns = matrix(plist('built-in','mdc3_ifo2ifo_noise'));
+
+% evalueat model
+[rw,cl]=size(Ns.objs);
+for ii=1:rw
+  for jj=1:cl
+    Ns.objs(ii,jj).setParams(parnames,parvalues);
+    Ns.objs(ii,jj).setXvals(f);
+    CSD(ii,jj)=eval(Ns.objs(ii,jj));
+    if ii==jj
+      CSD(ii,jj)= abs(CSD(ii,jj));
+    end
+  end
+end
+
+%% Plot
+
+plpl = plist('Linecolors', {'k', 'r', 'b'},...
+  'LineStyles', {'-', '--', ':'},...
+  'LineWidths', {3, 3, 3},...
+  'Legends', {'S_{o_1 o_1}', 'S_{o_1 o_\Delta}', 'S_{o_\Delta o_\Delta}'});
+
+% iplot(CSD(1,1),CSD(1,2),CSD(2,1),CSD(2,2))
+iplot(CSD(1,1),CSD(1,2),CSD(2,2),plpl)
+
+%%
+%--------------------------------------------------------------------------
+% 2 Dimensional Noisegen Filter
+%--------------------------------------------------------------------------
+
+
+%% 2 dim noisegen filter in z domain - symbolic precision All pass
+
+% plist for noise generation
+plns = plist(...
+  'csd',       CSD, ...
+  'targetobj',  'miir', ...
+  'fs',        fs, ...
+  'MaxIter',   70, ...
+  'PoleType',  3, ...
+  'MinOrder',  25, ...
+  'MaxOrder',  55, ...
+  'Weights',   2, ...
+  'Plot',      false,...
+  'MSEVARTOL', 1e-1,...
+  'FITTOL',    1e-4,...
+  'UseSym',    'on',...
+  'iunits',    'm',...
+  'ounits',    'm');
+
+na = matrix(plns);
+
+%% Check response of CSD
+
+na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
+na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
+na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
+na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
+
+mna = matrix([na11 na12;na21 na22]);
+
+ECSD = mna*conj(transpose(mna));
+
+for kk = 1:numel(CSD)
+  CSD(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
+  CSD(kk).setXunits(unit('Hz'));
+  ECSD.objs(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
+  ECSD.objs(kk).setXunits(unit('Hz'));
+end
+
+% calculate cross-coherence
+COH = CSD(1,2)./sqrt(CSD(1,1).*CSD(2,2));
+COH.simplifyYunits;
+RCOH = real(COH);
+ICOH = imag(COH);
+
+ECOH = ECSD.objs(1,2)./sqrt(ECSD.objs(1,1).*ECSD.objs(2,2));
+ECOH.simplifyYunits;
+RECOH = real(ECOH);
+IECOH = imag(ECOH);
+
+% set plist for plotting
+
+plpl = plist('Linecolors', {'k', 'r', 'b'},...
+  'LineStyles', {'-', '--', ':'},...
+  'LineWidths', {3, 3, 3},...
+  'Legends', {'Model', 'Fit', 'Residuals'});
+
+iplot(CSD(1,1),ECSD.objs(1,1),abs(CSD(1,1)-ECSD.objs(1,1)),plpl)
+iplot(CSD(2,2),ECSD.objs(2,2),abs(CSD(2,2)-ECSD.objs(2,2)),plpl)
+
+plpl = plist('Linecolors', {'k', 'r', 'b'},...
+  'LineStyles', {'-', '--', ':'},...
+  'LineWidths', {3, 3, 3},...
+  'Legends', {'Model', 'Fit', 'Residuals'},...
+  'Yscales', {'All', 'lin'},...
+  'Ylabels', {'All', 'Re(coh)'});
+
+iplot(RCOH,RECOH,abs(RCOH-RECOH),plpl)
+
+plpl = plist('Linecolors', {'k', 'r', 'b'},...
+  'LineStyles', {'-', '--', ':'},...
+  'LineWidths', {3, 3, 3},...
+  'Legends', {'Model', 'Fit', 'Residuals'},...
+  'Yscales', {'All', 'lin'},...
+  'Ylabels', {'All', 'Im(coh)'});
+
+iplot(ICOH,IECOH,abs(ICOH-IECOH),plpl)
+
+%% Check response of Filters
+
+% na11 = resp(na.objs(1,1).filters,plist('f',f,'bank','parallel'));
+% na12 = resp(na.objs(1,2).filters,plist('f',f,'bank','parallel'));
+% na21 = resp(na.objs(2,1).filters,plist('f',f,'bank','parallel'));
+% na22 = resp(na.objs(2,2).filters,plist('f',f,'bank','parallel'));
+
+
+% iplot(TF(1,1),na11,abs(TF(1,1) - na11))
+% iplot(TF(1,2),na12,abs(TF(1,2) - na12))
+% iplot(TF(2,1),na21,abs(TF(2,1) - na21))
+% iplot(TF(2,2),na22,abs(TF(2,2) - na22))
+
+iplot(na11,na12,na21,na22)
+
+%%
+%--------------------------------------------------------------------------
+% 2 Dimensional Noise Generation
+%--------------------------------------------------------------------------
+nsecs = 3e5;
+plng = plist(...
+  'model',   na,...
+  'Nsecs',   nsecs,...
+  'fs',      fs,...
+  'yunits',  'm');
+
+no = matrix(plng);
+
+%% Check psd
+
+no1 = no.objs(1);
+no2 = no.objs(2);
+
+no1xx = no1.psd;
+no2xx = no2.psd;
+
+iplot(no1xx,csd11,no2xx,csd22)
+
+%% Check cross-coherence
+
+no1 = no.objs(1);
+no2 = no.objs(2);
+nfft = 6e4;
+
+no12xx = cpsd(no1,no2,plist('nfft',nfft));
+no1xx = no1.psd(plist('nfft',nfft));
+no2xx = no2.psd(plist('nfft',nfft));
+%%
+
+ecoh = no12xx./(sqrt(no1xx.*no2xx));
+recoh = real(ecoh);
+iecoh = imag(ecoh);
+
+iplot(recoh,RECOH,plist('Yscales',{'All','lin'}))
+iplot(iecoh,IECOH,plist('Yscales',{'All','lin'}))
+
+%% Get averaged noise spectral curves
+
+N = 500;
+nsecs = 3e5;
+nfft = 6e4;
+
+Mnxx1 = 0; Mn2xx1 = 0;
+Mnxx2 = 0; Mn2xx2 = 0;
+
+Mnxy1 = 0; Mn2xy1 = 0;
+Mnxy2 = 0; Mn2xy2 = 0;
+
+tic
+for ii = 1:N
+  
+  fprintf('===== iter %s =====\n',num2str(ii))
+  
+  % noise generation
+  plng = plist(...
+    'model',   na,...
+    'Nsecs',   nsecs,...
+    'fs',      fs,...
+    'yunits',  'm');
+
+  no = matrix(plng);
+  
+  % periodogram calculation
+  no1 = no.objs(1);
+  no2 = no.objs(2);
+
+  no1xx = no1.psd;
+  no2xx = no2.psd;
+  f1 = no1xx.x;
+  
+  % running average and variance for periodogram
+  if ii == 1
+    Mnxx1 = no1xx.y;
+    Mnxx2 = no2xx.y;
+  else
+%     Qxx1 = no1xx.y - Mnxx1;
+    Mnxx1 = Mnxx1 + (no1xx.y - Mnxx1)./ii;
+    Mn2xx1 = Mn2xx1 + (no1xx.y - Mnxx1).*(no1xx.y - Mnxx1);
+    
+%     Qxx2 = no2xx.y - Mnxx2;
+    Mnxx2 = Mnxx2 + (no2xx.y - Mnxx2)./ii;
+    Mn2xx2 = Mn2xx2 + (no2xx.y - Mnxx2).*(no2xx.y - Mnxx2);
+  end
+  clear no1xx no2xx
+  
+  
+  % cross coherence calculation
+  no12xx = cpsd(no1,no2,plist('nfft',nfft));
+  no1xx2 = no1.psd(plist('nfft',nfft));
+  no2xx2 = no2.psd(plist('nfft',nfft));
+  f2 = no1xx2.x;
+
+  ecoh = (no12xx.y)./(sqrt((no1xx2.y).*(no2xx2.y)));
+  recoh = real(ecoh);
+  iecoh = imag(ecoh);
+  
+  clear no12xx no1xx2 no2xx2  
+  % running average and variance for cross coherence
+  if ii == 1
+    Mnxy1 = recoh;
+    Mnxy2 = iecoh;
+  else
+%     Qxy1 = recoh - Mnxy1;
+    Mnxy1 = Mnxy1 + (recoh - Mnxy1)./ii;
+    Mn2xy1 = Mn2xy1 + (recoh - Mnxy1).*(recoh - Mnxy1);
+    
+%     Qxy2 = iecoh - Mnxy2;
+    Mnxy2 = Mnxy2 + (iecoh - Mnxy2)./ii;
+    Mn2xy2 = Mn2xy2 + (iecoh - Mnxy2).*(iecoh - Mnxy2);
+  end
+  
+end
+toc
+
+% get mean and variance for psd 1
+Sxx1 = Mnxx1; % mean
+if N == 1
+  Svxx1 = Sxx1;
+else
+  Svxx1 = (Mn2xx1/(N-1))/N; % variance
+end
+% get mean and variance for psd 2
+Sxx2 = Mnxx2; % mean
+if N == 1
+  Svxx2 = Sxx2;
+else
+  Svxx2 = (Mn2xx2/(N-1))/N; % variance
+end
+% get mean and variance for real part of cross coherence
+Sxy1 = Mnxy1; % mean
+if N == 1
+  Svxy1 = Sxy1;
+else
+  Svxy1 = (Mn2xy1/(N-1))/N; % variance
+end
+% get mean and variance for imaginary part of cross coherence
+Sxy2 = Mnxy2; % mean
+if N == 1
+  Svxy2 = Sxy2;
+else
+  Svxy2 = (Mn2xy2/(N-1))/N; % variance
+end
+
+% save results
+
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat','Sxx1')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat','Svxx1')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat','Sxx2')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat','Svxx2')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy1.mat','Sxy1')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy1.mat','Svxy1')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy2.mat','Sxy2')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy2.mat','Svxy2')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f1.mat','f1')
+save('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f2.mat','f2')
+save(na,'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\na.mat')
+
+%% Load saved data
+
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxy2.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxy2.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\f2.mat'
+
+%% Plot spectra
+
+figure()
+
+p22 = loglog(f1(1:end-1),Sxx1(1:end-1) + sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
+hold on
+grid on
+p23 = loglog(f1(1:end-1),Sxx1(1:end-1) - sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
+p21 = loglog(f1(1:end-1),Sxx1(1:end-1),'r','LineWidth',3);
+
+p11 = loglog(ECSD.objs(1,1).x,ECSD.objs(1,1).y,'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+figure()
+
+p22 = loglog(f1(1:end-1),Sxx2(1:end-1) + sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
+hold on
+grid on
+p23 = loglog(f1(1:end-1),Sxx2(1:end-1) - sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
+p21 = loglog(f1(1:end-1),Sxx2(1:end-1),'r','LineWidth',3);
+
+p11 = loglog(ECSD.objs(2,2).x,ECSD.objs(2,2).y,'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+
+%% Plot coherence
+
+figure()
+
+p22 = semilogx(f2,Sxy1 + sqrt(Svxy1),'b','LineWidth',3);
+hold on
+grid on
+p23 = semilogx(f2,Sxy1 - sqrt(Svxy1),'b','LineWidth',3);
+p21 = semilogx(f2,Sxy1,'r','LineWidth',3);
+
+p11 = semilogx(RECOH.x,RECOH.y,'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('Re(coh) []','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+figure()
+
+p22 = semilogx(f2,Sxy2 + sqrt(Svxy2),'b','LineWidth',3);
+hold on
+grid on
+p23 = semilogx(f2,Sxy2 - sqrt(Svxy2),'b','LineWidth',3);
+p21 = semilogx(f2,Sxy2,'r','LineWidth',3);
+
+p11 = semilogx(IECOH.x,IECOH.y,'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('Im(coh) []','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+%% Check response of CSD
+
+currPath = cd;
+cd MDC1 % assuming to start from ..\ltpda\software\m-toolbox\test
+[TF,CSD] = mdc1_tf_models(plist('f',f1,'fs',fs));
+cd(currPath)
+
+na11 = resp(na.objs(1,1).filters,plist('f',f1,'bank','parallel'));
+na12 = resp(na.objs(1,2).filters,plist('f',f1,'bank','parallel'));
+na21 = resp(na.objs(2,1).filters,plist('f',f1,'bank','parallel'));
+na22 = resp(na.objs(2,2).filters,plist('f',f1,'bank','parallel'));
+
+mna = matrix([na11 na12;na21 na22]);
+
+ECSD = mna*conj(transpose(mna));
+
+for kk = 1:numel(CSD)
+  CSD(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
+  CSD(kk).setXunits(unit('Hz'));
+  ECSD.objs(kk).setYunits((unit('m')^2)*(unit('Hz')^-1));
+  ECSD.objs(kk).setXunits(unit('Hz'));
+end
+
+plpl = plist('Linecolors', {'k', 'r', 'b'},...
+  'LineStyles', {'-', '--', ':'},...
+  'LineWidths', {3, 3, 3},...
+  'Legends', {'Model', 'Fit', 'Residuals'});
+
+iplot(CSD(1,1),ECSD.objs(1,1),abs(CSD(1,1)-ECSD.objs(1,1)),plpl)
+iplot(CSD(2,2),ECSD.objs(2,2),abs(CSD(2,2)-ECSD.objs(2,2)),plpl)
+
+%% get expected response + window
+
+% we should considere that the window has a spectral response that can add
+% some systematics inside our spectral estimations. Since the window is a
+% product in time domain with the time series, that it corresponds to a
+% convolution in frequency domain. We can ifft the frequency response of
+% the filters, multiply in time domain for the window and then come back to
+% freq domain by fft.
+na = matrix('C:\Users\Luigi\Lavori\ltp_noisegen\matlab\na.mat');
+%%
+% filter responses
+fs = 10;
+nsecs = 3e5;
+nfft = nsecs.*fs;
+ff1 = (0:floor(nfft/2)).*fs./nfft;
+
+rsp11 = resp(na.objs(1,1).filters,plist('f',ff1,'bank','parallel'));
+rsp11 = rsp11.y;
+% rsp11 = [rsp11(1:end-1);flipud(conj(rsp11(2:end)))];
+rsp12 = resp(na.objs(1,2).filters,plist('f',ff1,'bank','parallel'));
+rsp12 = rsp12.y;
+% rsp12 = [rsp12(1:end-1);flipud(conj(rsp12(2:end)))];
+rsp21 = resp(na.objs(2,1).filters,plist('f',ff1,'bank','parallel'));
+rsp21 = rsp21.y;
+% rsp21 = [rsp21(1:end-1);flipud(conj(rsp21(2:end)))];
+rsp22 = resp(na.objs(2,2).filters,plist('f',ff1,'bank','parallel'));
+rsp22 = rsp22.y;
+% rsp22 = [rsp22(1:end-1);flipud(conj(rsp22(2:end)))];
+clear na
+
+% psd calc (npsd should contain nfft samples)
+npsd11 = rsp11.*conj(rsp11) + rsp12.*conj(rsp12);
+% npsd11 = [npsd11;flipud(npsd11(2:end-1))];
+npsd22 = rsp21.*conj(rsp21) + rsp22.*conj(rsp22);
+% npsd22 = [npsd22;flipud(npsd22(2:end-1))];
+
+%%
+
+% get window frequency response
+
+% BH92 in frequency domain
+W = BH92(2.*pi.*ff1./fs,nfft);
+W = W./max(abs(W));
+W = W.*conj(W);
+
+%% ifft dello spettro
+
+% npsd11 = [npsd11;flipud(npsd11(2:end-1))]./2;
+% npsd22 = [npsd22;flipud(npsd22(2:end-1))]./2;
+
+% g11 = ifft(npsd11);
+% g22 = ifft(npsd22);
+
+W = blackmanharris(nfft);
+if size(W,2)>1
+  W = W.';
+end
+W = W.^2;
+W = ifftshift(W);
+
+WS11 = fft(W.*g11);
+WS22 = fft(W.*g22);
+
+WS11 = 2.*WS11(1:numel(ff1));
+WS22 = 2.*WS22(1:numel(ff1));
+
+%% get convolution
+
+% WS11 = conv(W,npsd11);
+% WS11 = WS11(1:numel(ff1));
+
+
+%% Plot spectra
+
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx1.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Sxx2.mat'
+load 'C:\Users\Luigi\Lavori\ltp_noisegen\matlab\Svxx2.mat'
+%%
+figure()
+
+p22 = loglog(ff1(1:end-1),Sxx1(1:end-1) + sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
+hold on
+grid on
+p23 = loglog(ff1(1:end-1),Sxx1(1:end-1) - sqrt(Svxx1(1:end-1)),'b','LineWidth',3);
+p21 = loglog(ff1(1:end-1),Sxx1(1:end-1),'r','LineWidth',3);
+
+p11 = loglog(ff1(1:end),WS11(1:end),'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+%%
+figure()
+
+p22 = loglog(ff1(1:end-1),Sxx2(1:end-1) + sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
+hold on
+grid on
+p23 = loglog(ff1(1:end-1),Sxx2(1:end-1) - sqrt(Svxx2(1:end-1)),'b','LineWidth',3);
+p21 = loglog(ff1(1:end-1),Sxx2(1:end-1),'r','LineWidth',3);
+
+p11 = loglog(ff1(1:end-1),WS22(1:end-1),'k','LineWidth',3);
+
+
+xlabel('Frequency [Hz]','fontsize',18)
+ylabel('amplitude [m^2 Hz^{-1}]','fontsize',18)
+legend([p11(1) p21(1) p22(1)],'Fit Model', 'Data', 'Error \pm \sigma')
+
+%% get statistics
+
+Sn11 = (npsd11(14:end-1)-Sxx1(14:end-1))./npsd11(14:end-1);
+Sn11 = ao(cdata(Sn11));
+Sn22 = (npsd22(14:end-1)-Sxx2(14:end-1))./npsd22(14:end-1);
+Sn22 = ao(cdata(Sn22));
+
+% get equivalent normal distribution and histogram
+hpl1 = plist('N', 50);
+hpl2 = plist('N', 200);
+
+dS11h = hist(Sn11,hpl1);
+dS11h = dS11h./max(dS11h);
+dS11n = normdist(Sn11,hpl2);
+dS11n = dS11n./max(dS11n);
+
+dS22h = hist(Sn22,hpl1);
+dS22h = dS22h./max(dS22h);
+dS22n = normdist(Sn22,hpl2);
+dS22n = dS22n./max(dS22n);
+
+%% check
+
+iplot(dS11h,dS11n)
+iplot(dS22h,dS22n)
+
+%% Plot results of the statistics analysis
+
+figure()
+
+p1 = bar(dS11h.x,dS11h.y,'b');
+hold on
+grid on
+p2 = plot(dS11n.x,dS11n.y,'r','LineWidth',3);
+ 
+xlabel('Normalized Deviation [\Deltax / x]','fontsize',18)
+ylabel('Normalized Counts','fontsize',18)
+legend([p1(1) p2(1)],'Histogram', 'Equivalent Norm. Distr.')
+
+figure()
+
+p1 = bar(dS22h.x,dS22h.y,'b');
+hold on
+grid on
+p2 = plot(dS22n.x,dS22n.y,'r','LineWidth',3);
+ 
+xlabel('Normalized Deviation [\Deltax / x]','fontsize',18)
+ylabel('Normalized Counts','fontsize',18)
+legend([p1(1) p2(1)],'Histogram', 'Equivalent Norm. Distr.')
+
+
+
+%%
+m1 = mean(Sn11.y);
+m2 = mean(Sn22.y);
+s1 = std(Sn11.y);
+s2 = std(Sn22.y);
+
+
+
+
+
+%%
+
+nfft = 1e2;
+fs = 10;
+frq = linspace(0,pi,10000);
+% D=DirichletKernel(frq,nfft);
+% D=D./max(D);
+% figure;
+% plot(frq,20.*log10(D))
+% grid on
+W=BH92(frq,nfft);
+% W=W./max(W);
+figure;
+plot(frq./pi,20.*log10(W))
+grid on
+