clear all; warning off;
Fe=48000;
clear;
tlim=0.010;flim=50.0;
larg=464;
limite=10;% fraction du max du pulse servant à définir ses limites
cd('C:\Documents and Settings\JMLC\Bureau');
repertoire=uigetdir;
cd (repertoire);
[nomfichier,nomchemin]=uigetfile('*.wav','choix d''un fichier .wav');
nomchemin;nomfichier;
[sig,Fe,Nbits]=wavread([nomchemin nomfichier]);
tic;
Fe; Nbits;
ilim=tlim*Fe;
N=length(sig);
y2=abs(sig);
imax=find(y2==max(y2));
imin=imax-330000;
Ny = 2^19;
y = zeros(1,Ny);y2=y;
y =sig(imin+1 : imin+Ny); 
sig = y;
length(sig)
Fy=Fe;
fmax=23000; 
if fmax>Fe/2; fmax=Fe/2;end
%%###############################################
% definition du nombre de points sur l'axe des fréquences
Nfreqs=300;fmin=100;
if fmin<10; fmin=10;end %                      boucle 1 d & f
freq=ones(1,Nfreqs); 
%definition des frequences où sont calculés amplitude et %phase
freq(1)=fmin; dfreq=(fmax/fmin)^(1/(Nfreqs-1));
for i = 2:Nfreqs      %                        boucle 2  d
     freq(i)=freq(i-1)*dfreq;
 end                  %                        boucle 2  f
amplitude=zeros(5,Nfreqs);
%%###############################################
sigfft=fft(sig,2^19);
ih=round(200.0*Ny/48000)
    for i = 1:ih       %                        boucle 3  d
    sigfft(1+i)=0.0;
    sigfft(N-i)=0.0;
    end                %                        boucle 3  f
sig=real(ifft(sigfft,Ny));
sig=real(ifft(sigfft,Ny));
limf(1)=log10(5)*461720.46;
for i = 2 : 4
    limf(i)=max(1,round(limf(1)-(461720.46)*log10(i)));
end
limf(5)=1;
H1 = sig(limf(1):Ny);
H2 = sig(limf(2):limf(1));
H3 = sig(limf(3):limf(2));
H4 = sig(limf(4):limf(3));
H5 = sig(limf(5):limf(4));

%%###############################################

for h=1 : 5           %                        boucle 4   d
    if h==1;sig=H1;Fe=48000; end
    if h==2;sig=H2;Fe=48000/2; end
    if h==3;sig=H3;Fe=48000/3; end
    if h==4;sig=H4;Fe=48000/4; end
    if h==5;sig=H5;Fe=48000/5; end

%%###############################################
N=length(sig);
Ny=65536;
if N>Ny; Ny=N; end
y=zeros(1,Ny);y2=y;
y(1:N)=sig;
Fy=Fe;
% recentrage du signal actif dans la fenetre
temps=linspace(-((round(Ny/2)-1)/Fe),((round(Ny/2))/Fe),Ny);
ymax=abs(max(y)); ymin = abs(min(y));
if ymax<ymin; y=-y; end
y2=abs(y);
imax=find(y2==max(y2));
imin = imax;                        
dec=1+round(((imax+imin-2)-Ny)/2);  
if dec>0;y2(Ny-dec+1:Ny)=y(1:dec);y2(1:Ny-dec)=y(dec+1:Ny);end;
if dec<0;y2(-dec+1:Ny)=y(1:Ny+dec);y2(1:-dec)=y(Ny+dec+1:Ny);end
y=y2;
%%###############################################
y2=abs(y);
imax=find(y2==max(y2));
y(imax+ilim:Ny)=0.0;
y2=abs(y);
imax=find(y2==max(y2));
     y2(1)=0.0;
     for i = 2:Ny       %                        boucle 5   d
     y2(i) = y2(i-1) + y2(i);
     end                %                        boucle 5   f
ymax = y2(imax)+(0.40 * (max(y2)-y2(imax)));
ia_max= max(find(y2<ymax));
%###############################
if abs(ia_max-imax)<900;ia_max=imax+900;end
ia_min = round((1.5*imax)-(0.5*ia_max));
if ia_min < 1;  ia_min = 1;end
if ia_max > Ny; ia_max = Ny;end
imin=ia_min;
%%###############################################
dit = round(0.5+((ia_max-ia_min+1)/larg));                               %   NEW LINE
ia_max = imax+round((ia_max-imax)*larg*dit/(ia_max-ia_min+1));           %   NEW LINE
if ia_max >Ny ; ia_max =Ny;end
ia_min = ia_max-(dit*larg)+1;  
while ia_min<1;ia_min=ia_min+dit;end  
il = ia_min; ih = ia_max;
%%###############################################
colordef white;
set (0,'Units','pixels');
sc=get(0,'Screensize');
width=  600;
height = 450;
left=fix(sc(3)/20.00);
bottom=fix(sc(4)/20.00);
rect=[ left, bottom, width , height ];
set(figure(1),'Position',rect);  
outerpos=get(figure(1),'OuterPosition');
borders=outerpos-rect;
bottom= bottom + height;  %+ borders(4);
rect=[ left, bottom, width , height ];
set(figure(2),'Position',rect);  
bottom= bottom + height;  %+ borders(4);
rect=[ left, bottom, width , height ];
%set(figure(3),'Position',rect);  
height= 450;
left = left + width;
width = 600;
bottom=fix(sc(4)/20.00);
rect=[ left, bottom, width , height ];
set(figure(3),'Position',rect);  
bottom= bottom + height;
rect=[ left, bottom, width , height ];
set(figure(4),'Position',rect);  

%###############################################

%===============================================================
%paramètre de lissage des spectres
necart=12; coeff_fenetre=12;
nsmooth=10;
degsmooth=2;
%===============================================================
% necart: normal = 12, lissé = 64, brut=6
% nsmooth: normal= 10 ,lissé=20, brut=5
% degsmooth: normal=2, lissé=1; brut=3
%
% definition des fréquences limites de la fenêtre d'analyse fréquentielle
%
%
% initialisation de tableaux
Z_sin=zeros(1,Ny);
Z_cos=zeros(1,Ny);
argument=zeros(1,Nfreqs);
% ##################################
imax;
  for i=1:Nfreqs   
  %                        boucle 6   d
  omega=2*pi*freq(i);
  indice=round(Ny/(temps(Ny)*freq(i)));
      if indice==indice+rem(indice,2);indice=indice+1;end
  dind=coeff_fenetre*((indice-1)/2);
ampsin=0;ampcos=0;
  i_inf= imax-dind;i_sup=imax+dind;
  if i_inf<1;i_inf=1;end
  if i_sup>Ny;i_sup=Ny;end
         for k=i_inf:i_sup      %                        boucle 7   d
         ampsin=ampsin+y(k)*((1+cos(omega*temps(k)/coeff_fenetre))/2)*sin(omega*(temps(k)));         ampcos=ampcos+y(k)*((1+cos(omega*temps(k)/coeff_fenetre))/2)*cos(omega*(temps(k)));
         end                     %                        boucle 7
    amplitude(h,i)=20*log10(sqrt((ampsin^2)+(ampcos^2)));
     end              %                        boucle 6    f
if h==1;
ampmax=max(amplitude(1,:))-6;
ia_min=min(find(amplitude(1,:)>ampmax));
ia_max=max(find(amplitude(1,:)>ampmax));
ampmax = median(amplitude(1,ia_min:ia_max));
end  
% end
end                       %                        boucle 4   f

for h=1 : 5 
    sig= amplitude(h,:);
    sig = smooth(sig,nsmooth,'sgolay',degsmooth);
    amplitude(h,:)=sig';
end
amplitude=amplitude-ampmax;
for h=1:5
   Fe=48000/h/2.1;
   for i= 1:Nfreqs 
       if freq(i)<Fe; isup(h)=i+1.0;end
   end
end
   for i= 1:Nfreqs 
       if freq(i)<200; iinf=i+1.0;end
   end
 isup(1)=min(ia_max,isup(1));      
 ia_min=1; ia_max=Nfreqs;
%#################################################################""
ia_min=min(find(amplitude(1,:)>-45));ia_min=max(ia_min,1);
ia_max=max(find(amplitude(1,:)>-20));ia_max=min(ia_max,Nfreqs);
figure(1); semilogx(freq(iinf:isup(1)),amplitude(1,iinf:isup(1)),freq(iinf:isup(2)),amplitude(2,iinf:isup(2)),freq(iinf:isup(3)),amplitude(3,iinf:isup(3)),freq(iinf:isup(4)),amplitude(4,iinf:isup(4)),freq(iinf:isup(5)),amplitude(5,iinf:isup(5)),'LineWidth',2);grid;
axis([100 100000 -70 10]); grid on;
title('Frequency response curve');
xlabel('frequency  (Hz)');ylabel('L  (dB)');
legend('H1','H2','H3','H4','H5');
f1=freq(ia_min);f2=freq(ia_max);
%%###############################################
tlim=0.01;
[sig,Fe,Nbits]=wavread([nomchemin nomfichier]);
ilim=tlim*Fe;
N=65536 ; Ny=N;
y= zeros(1,N);
imax=find(abs(sig)==max(abs(sig)))
if sig(imax)<0.0;sig=-sig;end
% recentrage du signal actif dans la fenetre
y(1:65536)=sig(imax-32767:imax+32768);
Fe=48000; Fy=48000;
ih=round(20.0*Ny/48000);
temps=linspace(-((round(Ny/2)-1)/Fe),((round(Ny/2))/Fe),Ny);
%%###############################################
y2=abs(y);
imax=find(y2==max(y2));
y(imax+ilim:Ny)=0.0;
y2=abs(y);
imax=find(y2==max(y2));
y2(1)=0.0;
for i = 2:Ny
    y2(i) = y2(i-1) + y2(i);
end
ymax = y2(imax)+(0.40 * (max(y2)-y2(imax)));
ia_max= max(find(y2<ymax));
larg=464;coef=2.0;
%##########################################################################
if abs(ia_max-imax)<ceil(larg/coef);ia_max=imax+ceil(larg/coef);end  
ia_min = round((1.5*imax)-(0.5*ia_max));
if ia_min < 1;  ia_min = 1;end
if ia_max > Ny; ia_max = Ny;end
imin=ia_min;
%%###############################################
dit = round(0.5+((ia_max-ia_min+1)/larg));                              %   NEW LINE
ia_max = imax+round((ia_max-imax)*larg*dit/(ia_max-ia_min+1));          %   NEW LINE
if ia_max >Ny ; ia_max =Ny;else;end
ia_min = ia_max-(dit*larg)+1;  
while ia_min<1;ia_min=ia_min+dit;end  
il = ia_min; ih = ia_max;
%%###############################################
%paramètre de lissage des spectres
necart=12; coeff_fenetre=12;
nsmooth=10;
degsmooth=2;
%===============================================================
% necart: normal = 12, lissé = 64, brut=6
% nsmooth: normal= 10 ,lissé=20, brut=5
% degsmooth: normal=2, lissé=1; brut=3
%
% definition des fréquences limites de la fenêtre d'analyse fréquentielle
   fmax=23000; 
if fmax>Fe/2; fmax=Fe/2;end
% definition du nombre de points sur l'axe des fréquences
Nfreqs=300;fmin=10;
if fmin<10; fmin=10;end
freq=ones(1,Nfreqs); 
%definition des frequences où sont calculés amplitude et %phase
freq(1)=fmin; dfreq=(fmax/fmin)^(1/(Nfreqs-1));
for i = 2:Nfreqs
     freq(i)=freq(i-1)*dfreq;
end
% initialisation de tableaux
Z_sin=zeros(1,Ny);
Z_cos=zeros(1,Ny);
amplit=zeros(1,Nfreqs);
argument=zeros(1,Nfreqs);
% boucle principale
%imax
  for i=1:Nfreqs
  omega=2*pi*freq(i);
  indice=round(Ny/(temps(Ny)*freq(i)));
      if indice==indice+rem(indice,2);indice=indice+1;end
  dind=coeff_fenetre*((indice-1)/2);
ampsin=0;ampcos=0;
  i_inf= imax-dind;i_sup=imax+dind;
  if i_inf<1;i_inf=1;end
  if i_sup>Ny;i_sup=Ny;end
         for k=i_inf:i_sup
         ampsin=ampsin+y(k)*((1+cos(omega*temps(k)/coeff_fenetre))/2)*sin(omega*(temps(k)));
         ampcos=ampcos+y(k)*((1+cos(omega*temps(k)/coeff_fenetre))/2)*cos(omega*(temps(k)));
         end
    amplit(i)=20*log10(sqrt((ampsin^2)+(ampcos^2)));
    argument(i)=atan2(ampsin,ampcos);
     end
  %  depliement de la phase, normalisation module
  phase = -(180/pi)*argument';
phase_depliee=(180/pi)*unwrap(argument,pi);
ampmax=max(amplit)-6;
ia_min=min(find(amplit>ampmax));
ia_max=max(find(amplit>ampmax));
ampmax = median(amplit(ia_min:ia_max));
amplit=amplit-ampmax;
amplit=smooth(amplit,nsmooth,'sgolay',degsmooth);
ia_min=1; ia_max=Nfreqs;
phase_depliee=smooth(phase_depliee,20,'sgolay',1);
retard=zeros(1,Nfreqs);
 vit=344/360;
for i =2:Nfreqs-1
    retard(i)=vit*(phase_depliee(i+1)-phase_depliee(i-1))/(freq(i+1)-freq(i-1));
end
retard(1)=2*retard(2)-retard(3);
retard(Nfreqs)=2*retard(Nfreqs-1)-retard(Nfreqs-2);
ia_min=min(find(amplit>-15));ia_min=max(ia_min,1);
ia_max=max(find(amplit>-20));ia_max=min(ia_max,Nfreqs);
figure(2); semilogx(freq(ia_min:ia_max),retard(ia_min:ia_max),'LineWidth',2);% grid;
axis([100 100000 -0.1  0.5]); grid on;
title('Group delay');
xlabel('frequency  (Hz)');ylabel('equivalent distance to the GD (m)');
f1=freq(ia_min);f2=freq(ia_max);
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
%  pause
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
% = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =
 %===============================================================
%      limits of frequency of the spectrogram
%===============================================================
%      largeur 464,  hauteur = 366
fmin=200;
fmax=20000; 
if fmax>Fe/2; fmax=Fe/2;end
Nfreqs=366;
if f1<fmin; f1=fmin;end
freq=ones(1,Nfreqs); sfen=ones(1,Nfreqs);logfreq=zeros(1,Nfreqs);
freq(1)=fmin; dfreq=(fmax/fmin)^(1/(Nfreqs-1));
for i = 2:Nfreqs
     freq(i)=freq(i-1)*dfreq;
end
logfreq=log10(freq);  % for vertical axis of the spectrogram
%===============================================================
t = temps(il:ih);           
y = y(il:ih);
N=length(t);
%************************************************************************
%   Spectrogram calculation
%************************************************************************
   spectro=zeros(larg,Nfreqs);
   ts=zeros(larg,1);
    for i=1:larg
    ts(i)=t(1)+i*dit/Fe;   %new time scale of the spectrogram with 481 steps
    end
   
   spectro=spectro-40.0;
   specm=0.0;
%#################    MAIN LOOP   #########################"
for j =1 : Nfreqs
    spsin=zeros(N,1);spcos=zeros(N,1);
    Fr = freq(j);
    omega= 2*pi*Fr;
    Nw=round(Fe/Fr);
    if Nw > N; Nw=N-2;end
            spsin(1)=0.0; spcos(1)=0.0;          
            for i = 2 : N
            spsin(i)=spsin(i-1)+(y(i)*sin(omega*t(i)));
            spcos(i)=spcos(i-1)+(y(i)*cos(omega*t(i)));
            end
      for i = 2 : N-Nw
      spsin(i)= spsin(i-1)-spsin(i)+spsin(i+Nw);
      spcos(i)= spcos(i-1)-spcos(i)+spcos(i+Nw);
      end
          
           spsin(1)=0.0; spcos(1)=0.0;
           for i = N-Nw : N
           spsin(i)= spsin(i-1);
           spcos(i)= spcos(i-1);
           end
     
      for i = 2 : N-Nw
      spsin(i)= spsin(i-1)-spsin(i)+spsin(i+Nw);
      spcos(i)= spcos(i-1)-spcos(i)+spcos(i+Nw);
      end
      
           spsin(1)=0.0; spcos(1)=0.0;
           for i = N-Nw : N
           spsin(i)= spsin(i-1);
           spcos(i)= spcos(i-1);
           end
   
     spsin(1)=0.0; spcos(1)=0.0;
     for i = 1: larg
          ia= 1+((i-1)*dit)-round((1.5*Nw))  ;
          ib= ia + Nw  ;
     if ia<1; ia=1;end
     if ia>N; ia=N;end
     if ib<1; ib=1;end
     if ib>N; ib=N;end
    spec = 20*log10((1/(Nw*Nw))*sqrt(((spsin(ib)-spsin(ia)))^2+((spcos(ib)-spcos(ia))^2)));
 
if spec>specm; specm=spec;end
     spectro(i,j)= spec;
     end
 end
%#################  END OF MAIN LOOP   #########################"

%===============================================================
%  normalization inside the interval [0 , -40db]
       for j =1 : Nfreqs
       for i = 1 : larg
       spectro (i,j) = spectro (i,j) - specm;
       end
       end
  for j =1 : Nfreqs
  for i = 1 : larg
  if spectro (i,j)<-40; spectro (i,j) = -40.0; end
  end
  end
%===============================================================
%===============================================================
%           graphical output of the spectrogram 
%===============================================================
t = 1000.00*t;
YT=[40 60 80 100 200 400 600 800 1000 2000 4000 6000 8000 10000];
YT= log10(YT);
Ymax= .1*fix(1+max(y)*10.0); Ymin= .1*fix(-1+min(y)*10.0);Ylength =Ymax-Ymin;
Yt=linspace(Ymin,Ymax,fix(1+(Ymax-Ymin)/0.1));    
figure(3); plot (t,y,'LineWidth',2);  %set(gca,'YTick',Yt);
axis([t(1) t(N) Ymin Ymax]);
grid on;
title('Impulse response');
xlabel('time (ms)');ylabel('signal');

figure(4);clf;orient tall;zoom on
imagesc(t,logfreq,spectro'); axis xy;colormap(jet(256));
set(gca,'YTick',YT);
set(gca,'YTickLabel',{'4','6','8','100','2','4','6','8','1000','2','4','6','8','10000'})
grid on;
title('Spectrogram');
xlabel('time (ms)');ylabel('frequency'); % colorbar
%===============================================================
%===============================================================
toc;                    % clock stop
%************************************************************************
%                             END
%************************************************************************ 
cd('C:\MATLAB6p5p1\mes_fonctions');