Cuki_IR_gen1s.m

% Acoustic guitar IR generator "light" 1.0
% Author: Kien Phan Huy on the 14th of July 2020. Copyright 2020
% Uses the Min Phase response computation from J.O. Smith, 1982-2002
% Uses amodified versiona of Oct_spectrum from M. Buzzoni, Dec. 2018
% Uses IIR butterworth coefficient computation from Neil Robertson , 12/29/17
% Uses embbeded waitbar Yuanfei (2020). Embedding Waitbar Inside A GUI (https://www.mathworks.com/matlabcentral/fileexchange/47896-embedding-waitbar-inside-a-gui), MATLAB Central File Exchange. Retrieved July 15, 2020.
% See copyrights for subfunctions and roaotine below
% Graphic user interface updated to Octave 7.1 by pokibg
function Cuki_IR_gen1s()
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Position ui controls
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    f=[143 129 1024 768];
  Nu=[1024 768 1024 768];
  u1=[305   115   665   430]./Nu; % set(gca,'units','pixels','position',[305   115   665   430])
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  u2=[20 740 280 30]./Nu;
  u3=[300 740 300 30]./Nu;
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hh=-15; % decalage tout sauf d?but
  h2=-30; % d?calage sous interface audio
  dd=100; % d?calage droite interface audio
  u4=[450 635+hh+h2 500 20]./Nu;;
  u5=[20 710+hh 500 30]./Nu;
  u6=[50 680+hh dd 30]./Nu;
  u7=[20+dd 680+hh 500 30]./Nu;
  u8=[50 640+hh dd 30]./Nu;
  u9=[20+dd 640+hh 500 30]./Nu;
  u10=[550+dd 710+hh 250 30]./Nu; % uicontrol('style','text', 'string','3) Choose Frequency sampling:', ...
  u11=[550+dd 680+hh 250 30]./Nu;
  u12=[550+dd 640+hh 250 30]./Nu;
  u13=[20 630+hh+h2 400 30]./Nu;
  u14=[20 600+hh+h2 150 30]./Nu;
  u15=[20 550+hh+h2 220 30]./Nu;
  u16=[150 520+hh+h2 100 30]./Nu;
  u17=[20 520+hh+h2 100 30]./Nu;
  u18=[20 470+hh+h2 200 30]./Nu;
  u19=[20 440+hh+h2 200 30]./Nu;
  u20=[20 390+hh+h2 200 30]./Nu;
  u21=[20 360+hh+h2 150 30]./Nu; % uicontrol(dlg, 'style','pushbutton', 'string','Compute IR', 'position',[20 360+hh+h2 150 30]
  u22=[20 310+hh+h2 200 30]./Nu;
  u23=[20 280+hh+h2 70 30]./Nu;
  u24=[100 280+hh+h2 70 30]./Nu;
  u25=[180 280+hh+h2 70 30]./Nu;
  u26=[20 230+hh+h2 200 30]./Nu;
  u27=[20 200+hh+h2 230 30]./Nu;
  u28=[20 160+hh+h2 150 30]./Nu;
  u29=[20 160+hh+h2 72 30]./Nu;
  u30=[100 160+hh+h2 150 30]./Nu; % uicontrol(dlg, 'style','pushbutton', 'string','Save',
  u31=[20 120+hh+h2 150 30]./Nu;
  u32=[20 20 150 30]./Nu;
  u33=[200 20 150 30]./Nu;
  u34=[730 20 500 30]./Nu;
  u35=[800 740 150 30]./Nu;
##  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation figure
  dlg = figure('name','Cuki IR generator light v1.0', 'position',f, 'menubar','none', 'numbertitle','off');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation axes
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %axes1=axes(dlg,'position',[200,100,640,480])
  %h=plot(dlg,1:10,1:10,'-o')
  ax = axes('NextPlot', 'add');
  grid (ax, 'minor');
  %set(gca,'position',[0.3   0.15   0.65   0.6])
  set(gca,'position',u1)
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   uicontrol('style','text', 'string','1) Select your configuration:', ...
      'units','normalized','position',u2, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
    % Creation du menu deroulant
  icfg=uicontrol(dlg,'style','popupmenu', 'string',{'Pickup in CH1, Mic in CH2','Pickup in CH2, Mic in CH1'}, ...
      'units','normalized','position',u3);
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  hh=-15; % decalage tout sauf d?but
  h2=-30; % d?calage sous interface audio
  dd=100; % d?calage droite interface audio
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Wait bar
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ax2=axes('NextPlot', 'add');
  set(gca,'Position',u4);
  %ax2=axes(dlg,'Position',[0.45 0.8 0.4 0.03]);
  set(ax2,'Xtick',[],'Ytick',[],'Xlim',[0 1000]);
  box on;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation liste des interface audio
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','2) Choose audio interface+driver:', ...
      'units','normalized','position',u5, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  nDevices = audiodevinfo(0);
  devinfo = audiodevinfo();
  mot=['liste= {''',devinfo.input(1).Name,''''];
  for i=2:size(devinfo.input,2),
     mot=[mot,',''',devinfo.input(i).Name,''''];
  end
  mot=[mot,'};'];
  eval(mot);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','Input:', ...
      'units','normalized','position',u6, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation du menu deroulant INPUT
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ai=uicontrol(dlg,'style','popupmenu', 'string',liste, ...
      'units','normalized','position',u7, 'callback', @cbk_popupmenu);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  mot=['liste2= {''',devinfo.output(1).Name,''''];
  for i=2:size(devinfo.output,2),
     mot=[mot,',''',devinfo.output(i).Name,''''];
  end
  mot=[mot,'};'];
  eval(mot);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','Output:', ...
      'units','normalized','position',u8, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation du menu deroulant OUTPUT
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  ao=uicontrol(dlg,'style','popupmenu', 'string',liste2, ...
      'units','normalized','position',u9);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation liste des frequences d'echantillonage
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','3) Choose Frequency sampling:', ...
      'units','normalized','position',u10, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  fsliste=uicontrol(dlg,'style','popupmenu', 'string',{'44100','48000'}, ...
      'units','normalized','position',u11);
  bitdepth=uicontrol(dlg,'style','popupmenu', 'string',{'16 bits','24 bits'}, ...
      'units','normalized','position',u12,'value',2);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Record quick, test levels
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','4) Test recording and check levels (optional):', ...
      'units','normalized','position',u13, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');

  uicontrol(dlg, 'style','pushbutton', 'string','Record 10s','units','normalized', 'position',u14, 'callback',{@record20s, ax,ai,devinfo,fsliste,ax2,bitdepth });
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Record long to make IR
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','5) Record for IR generation:', ...
      'units','normalized','position',u15, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  Trec=uicontrol(dlg,'style','popupmenu', 'string',{'1 min','2 min'}, ...
      'units','normalized','position',u16);

  uicontrol(dlg, 'style','pushbutton', 'string','Record','units','normalized', 'position',u17, 'callback',{@record2mn, ax,ai,devinfo,fsliste,Trec,ax2,bitdepth});
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Creation liste des formats de sortie
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','6) IR file format:', ...
      'units','normalized','position',u18, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  listfmt={'1024 pts, 16 bits','2048 pts, 16 bits','8096 pts, 16 bits','1024 pts, 24 bits','2048 pts, 24 bits','8096 pts, 24 bits'};
  fsfmt=uicontrol(dlg,'style','popupmenu', 'string',listfmt, ...
      'units','normalized','position',u19,'value',2);
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Compute IR
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','7) Compute IR:', ...
      'units','normalized','position',u20, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  uicontrol(dlg, 'style','pushbutton', 'string','Compute IR','units','normalized', 'position',u21, 'callback',{@Go, ax,ai,devinfo,fsfmt,icfg,ax2});
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Listen
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uilisten=uicontrol('style','text', 'string','8) Listen:', ...
      'units','normalized','position',u22, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left','value',0);
  uicontrol(dlg, 'style','togglebutton', 'string','Mic','units','normalized', 'position',u23, 'callback',{@Mic,ax,ai,devinfo,fsfmt,icfg,ao,uilisten});
  uicontrol(dlg, 'style','togglebutton', 'string','Pickup','units','normalized', 'position',u24, 'callback',{@Pickup,ax,ai,devinfo,fsfmt,icfg,ao,uilisten});
  uicontrol(dlg, 'style','togglebutton', 'string','IR','units','normalized', 'position',u25, 'callback',{@IR,ax,ai,devinfo,fsfmt,icfg,ao,uilisten});
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Save IR
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol('style','text', 'string','8) Save IR file:', ...
      'units','normalized','position',u26, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
  Svfmt=uicontrol(dlg,'style','popupmenu', 'string',{'Standard','Various flavours','Feedback friendly','JF45 Copyright: Jon Fields'}, ...
      'units','normalized','position',u27);
  if ispc
      uicontrol(dlg, 'style','pushbutton', 'string','Save','units','normalized', 'position',u28, 'callback',{@Save,ax,ai,devinfo,fsfmt,icfg,listfmt,Svfmt});
  endif
  if ismac
      uicontrol('style','text', 'string','Filename:', ...
      'units','normalized','position',u29, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
      Fnam=uicontrol('style','edit', 'string','IR.wav', ...
      'units','normalized','position',u30, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
      uicontrol(dlg, 'style','pushbutton', 'string','Save','units','normalized', 'position',u31, 'callback',{@Save,ax,ai,devinfo,fsfmt,icfg,listfmt,Svfmt,Fnam});
  endif
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Close dlg
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol(dlg, 'style','pushbutton', 'string','Close', 'units','normalized', 'position',u32, 'callback',{@Close,ax,dlg});
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Donate
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  uicontrol(dlg, 'style','pushbutton', 'string','Donate','units','normalized',  'position',u33, 'callback',{@Donate,dlg});
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Paypal gift
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    uicontrol('style','text', 'string','Copyright: Kien Phan Huy, July 2020', ...
      'units','normalized', 'position',u34, ...
      'foregroundcolor', 'b', 'backgroundcolor','w', ...
      'fontangle','italic', 'horizontalalignment','left');
####  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
####  % Load file (DEBUG button)
####  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##  uicontrol(dlg, 'style','pushbutton', 'string','Load File', 'position',[800 740 150 30], 'callback',{@Load,ax,ai,devinfo,fsfmt,icfg,ao,bitdepth});
####  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
##  pause(0.1);
##  set(dlg,'position',[143 129 1024 768]); % pour Mac

end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Load y file debug
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Load(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,ao,bitdepth)
  [filename, pathname] = uigetfile({'*.wav','wav files'}, 'Choose a y file');
  %Fname='rec24_48000GtL_MicR_.wav';
  [y,FS]=audioread([pathname,filename]);

  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Set bitdepth
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  indxbits=get(bitdepth,'value');
  if indxbits==1
    nbits=16;
  else
    nbits=24;
  endif;

  handles.nbits=nbits;
  handles.y=y;
  handles.fs=FS;
  guidata(hObject,handles);
end;

function cbk_popupmenu(hObject, eventdata)
%  menu_itemms = get(h,'string');
%  disp(menu_itemms{get(hObject,'value')})
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recording 10 sec
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function record20s(hObject, eventdata,ax,ai,devinfo,fsliste,ax2,bitdepth)
%  disp('record')
  indx=get(ai,'value');
  id=devinfo.input(indx).ID;
  channels=2;
  %nbits=24;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  indxfs=get(fsliste,'value');
  if indxfs==1
    fs=44100;
    %display('44100');
  else
    fs=48000;
    %display('48000');
  endif;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Set bitdepth
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  indxbits=get(bitdepth,'value');
  if indxbits==1
    nbits=16;
  else
    nbits=24;
  endif;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  recorder = audiorecorder (fs, nbits, channels, id);
  length=10;
  record(recorder);
  frac=0;
  %h = waitbar (frac,'Recording...');
  for i=1:length
    axes(ax2)
    cla
	  rectangle('Position',[0,0,(round(1000*i/length))+1,20],'FaceColor','b');
    text(480,10,[num2str(round(100*i/length)),'%']);
     %waitbar (i/length, h);
     pause(1)
  end;
  stop(recorder)
  %close(h)

  y = getaudiodata (recorder);
  Nb=size(y,1);
  t=(1:Nb)/fs;
  nnn=1:100:Nb;

  axes(ax)
  cla
  set(ax, 'XScale', 'linear')
  xticks('auto')
  xticklabels('auto')
  xlim([0 10])
  ylim([-1 1])
  xlabel('time(s)')
  ylabel('amplitude (a.u.)')
  plot(ax,t(nnn),y(nnn,1),t(nnn),y(nnn,2));
  plot(ax,xlim,[0.1 0.1],'Color','k','linewidth',3)
  plot(ax,xlim,[-0.1 -0.1],'Color','k','linewidth',3)
  if max(max(abs(y)))<=0.1
    title('The signal is too weak (within black lines), increase the gain','Color',[1 0 0])
  elseif max(max(abs(y)))>=1
    title('Distortion decteted: decrease the gain','Color',[1 0 0])
  else
    title('Check the signal is not too weak (see black lines)','Color',[0 0 0])
  endif;
  clear recorder

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Recording 1-2 min
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function record2mn(hObject, eventdata,ax,ai,devinfo,fsliste,Trec,ax2,bitdepth)
%  disp('record')
  indx=get(ai,'value');
  id=devinfo.input(indx).ID;
  channels=2;
  %nbits=24;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Set frequency sampling
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  indxfs=get(fsliste,'value');
  if indxfs==1
    fs=44100;
    %display('44100');
  else
    fs=48000;
    %display('48000');
  endif;

  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Set bitdepth
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  indxbits=get(bitdepth,'value');
  if indxbits==1
    nbits=16;
  else
    nbits=24;
  endif;
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Set Rec time
  %%%%%%%%%%%%%%%%%%%%%%%%%%%
  T=get(Trec,'value');
  if T==1
    length=1*60;
  else
    length=2*60;
  endif;


  recorder = audiorecorder (fs, nbits, channels, id);

  record(recorder);
  frac=0;
  %h = waitbar (frac,'Recording...');
  for i=1:length
     %waitbar (i/length, h);
    axes(ax2)
    cla
	  rectangle('Position',[0,0,(round(1000*i/length))+1,20],'FaceColor','b');
    text(480,10,[num2str(round(100*i/length)),'%']);

     pause(1)
  end;
  stop(recorder)
  %close(h)
  y = getaudiodata (recorder);
  Nb=size(y,1);
  t=(1:Nb)/fs;
  nnn=1:1000:Nb;

  axes(ax)
  cla
  set(ax, 'XScale', 'linear')
  xticks('auto')
  xticklabels('auto')
  xlim([0 length])
  ylim([-1 1])
  xlabel('time(s)')
  ylabel('amplitude (a.u.)')
  plot(ax,xlim,[0.1 0.1],'Color','k','linewidth',3)
  plot(ax,xlim,[-0.1 -0.1],'Color','k','linewidth',3)
  plot(ax,t(nnn),y(nnn,1),t(nnn),y(nnn,2));
##  if max(max(abs(y)))<=0.1
##    title('The signal is too weak (within black lines), increase the gain','Color',[1 0 0])
##  elseif max(max(abs(y)))>=1
##    title('Distortion decteted: decrease the gain','Color',[1 0 0])
##  else
##    title('Check the signal is not too weak (see black lines)','Color',[0 0 0])
##  endif;
  clear recorder
  handles.y=y;
  handles.fs=fs;
  handles.nbits=nbits;
  guidata(hObject,handles);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Go: Compute IR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Go(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,ax2)
  frac=0;
  %h = waitbar (frac,'Computing...');
  axes(ax2)
  cla
  rectangle('Position',[0,0,(round(1000*0))+1,20],'FaceColor','b');
  text(480,10,[num2str(round(100*0)),'%']);
  pause(0.01)
  handles = guidata(hObject);
  y=handles.y;
  fs=handles.fs;
  nbits=handles.nbits;
  Nb=size(y,1);

  indxinput=get(icfg,'value');
  if indxinput==1
    MIC=y(:,2); %mic;
    PIC=y(:,1); %pic;
  else
    MIC=y(:,1); %mic;
    PIC=y(:,2); %pic;
  endif
  Nb=size(y,1);
  indxfmt=get(fsfmt,'value');
  if (indxfmt == 1)||(indxfmt == 4)
  NbF=1024;
  endif
  if (indxfmt == 2)||(indxfmt == 5)
  NbF=2048;
  endif
  if (indxfmt == 3)||(indxfmt == 6)
  NbF=8096;
  endif

  Nbuff=fs;
  Nbmax=floor(Nb/Nbuff)-10;
  FS=fs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% FIR calculaiton
FIR1 =zeros(NbF,1);
clear alice;

for n=1:Nbmax, % 200
        q=n;
        i=2*FS;
        i=i+n*Nbuff;

%        window = (.42-.5*cos(2*pi*(0:Nbuff-1)/(Nbuff-1))+.08*cos(4*pi*(0:Nbuff-1)/(Nbuff-1)))';
%        FIR=fft(MIC(i:i+Nbuff-1).*window,NbF)./fft(PIC(i:i+Nbuff-1).*window,NbF);
        FIR=fft(MIC(i:i+Nbuff-1),NbF)./fft(PIC(i:i+Nbuff-1),NbF);
        IR=zeros(Nbuff,1);
        IR(1:NbF)=real(ifft(FIR));
        IR=IR/max(abs(IR));

        if ((max(isinf(FIR))==1)||(max(isnan(FIR))==1))
            IR=zeros(Nbuff,1);
            IR(1)=1;
            FIR=fft(IR,NbF);
            display('NaN or Inf');
        end;


        alice(1:NbF,q)=FIR;
        FIR1=FIR1+FIR;
        if mod(round(n/Nbmax*100),10)==0
                  %round(n/Nbmax*100)
                  %waitbar (n/Nbmax/3, h);
                  axes(ax2)
                  cla
	                rectangle('Position',[0,0,(round(1000*n/Nbmax/3))+1,20],'FaceColor','b');
                  text(480,10,[num2str(round(100*n/Nbmax/3)),'%']);
                  pause(0.01)
            %waitbar(n/Nbmax,f,'Processing...');
        end;
end;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Post selection Processing
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

NN=size(alice,1);
    for i=1:NN, % POur chaque frequence
        a=alice(i,:);

        A=a(abs(abs(a)-mean(a))<2*std(a)); % On ne garde que les echantillons compris dans 2 ecarts types

        ALICE(i)=mean(A);
        if isnan(mean(A)),
            %display('Algorithm Failed: Increase the recording length and/or improve the S/N ratio');
            ALICE(i)=1;
        end
        Q(i)=length(A);
        if mod(round(i/NN*100),10)==0
                  %round(n/Nbmax*100)
                  %waitbar (1/3+i/NN/3, h);
                  axes(ax2)
                  cla
	                rectangle('Position',[0,0,(round(1000*(1/3+i/NN/3)))+1,20],'FaceColor','b');
                  text(480,10,[num2str(round(100*(1/3+i/NN/3))),'%']);
                  pause(0.01)
            %waitbar(n/Nbmax,f,'Processing...');
        end;
    end;
dnuX=FS/NbF;
nuX=(-NbF/2:NbF/2-1)*dnuX;
nn=NbF/2+1:NbF;

window = (.42-.5*cos(2*pi*(0:2*NbF-1)/(2*NbF-1))+.08*cos(4*pi*(0:2*NbF-1)/(2*NbF-1)))';
blackmanwin=window(NbF+1:end);
ir2=ifft(ALICE);
ir2=ir2.*blackmanwin';
ir2=ir2/max(abs(ir2))*0.95;
IR2(1:NbF,1)=real(ir2);% raw IR generated

%plot(IR2)
IR0=IR2;
nn=(10*FS+1:20*FS);
MS=MIC(nn,1);
PS=PIC(nn,1);

%waitbar(0.6, h);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate Two Octave spectrum (Mic & Pic convolved)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%pkg load signal;
PS=conv(PIC(nn,1),IR0,'same');

[p,cf,overall_lev,overall_levA,sfilt] = oct_spectrum2(MS/max(abs(MS)),FS,3,1,1,0);
[p2,cf2,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS/max(abs(PS)),FS,3,1,1,0);

pp=[p+100;p2+100];

##figure(10)
##c=bar(log10(cf),pp')
##%c=bar(log10(cf),p+100)
##grid on
##legend('Mic','Convolved Pickup')
## xlabel('Frequency 10^n (Hz)')
##ylabel('dB')
##figure(3)
##plot(10*log10(cf),p-p2)

g0=p-p2; % SB-SB2;
%dgain=zeros(31,1); % il est l? le probl?me
dgain=zeros(length(cf),1); % il est l? le probl?me


for qq=1:3, % iteration GEQ fit
    IRX=zeros(NbF,1);
    IRX(1)=1;
    IR1=IR0;

    for i=1:length(p) % Browse each frequency of the octave spectrum
        g=g0(i)+dgain(i); % dgain terme correctif GEQ fit
        Q=4;
        fc=cf(i);


        % IIR coefficient calculation
        g=10^(g/20);
        t0=2*pi*fc/fs;
        if g >= 1
            beta=t0/(2*Q);
        else
            beta=t0/(2*g*Q);
        endif
        a2=-0.5*(1-beta)/(1+beta);
        a1=(0.5-a2)*cos(t0);
        b0=(g-1)*(0.25+0.5*a2)+0.5;
        b1=-a1;
        b2=-(g-1)*(0.25+0.5*a2)-a2;

        % SOS Form
        b=2*[b0 b1 b2];
        a=[1 -2*a1 -2*a2];
        IRX=filter(b,a,IRX); % test
        IR1=filter(b,a,IR1); % final

        mot=['b',num2str(i,'%02d'),'=b;'];
        eval(mot);
        mot=['a',num2str(i,'%02d'),'=a;'];
        eval(mot);

    end;


    dnuX=FS/NbF;
    nuX=(-NbF/2:NbF/2-1)*dnuX;
    nn=NbF/2+1:NbF;
    s0=20*log10(abs(fftshift(fft(IRX)))); % compute the GEQ spectrum

##    figure(1)
##    plot(10*log10(nuX(nn)),s0(nn))


    dxx=(nuX(2)-nuX(1))/2;
    for i=1:length(cf)
        n=(abs(nuX-cf(i))<dxx);
        gainGEQ(i)=s0(n); % compute the GEQ spectrum at the octave frequencies
    end;

##    %cl=lines(qq+1);
##    figure(3)
##    hold on
##    %plot(10*log10(cf),gainGEQ,'o-')
##    plot(10*log10(nuX(nn)),s0(nn))
##    hold off

    dgain=dgain-(gainGEQ'-g0')/2;

    %waitbar(2/3+qq*0.1, h);
    axes(ax2)
    cla
	  rectangle('Position',[0,0,(round(1000*(2/3+qq*0.1)))+1,20],'FaceColor','b');
    text(480,10,[num2str(round(100*(2/3+qq*0.1))),'%']);
    pause(0.01)
end;

##    figure(2)
##    plot(10*log10(nuX(nn)),s0(nn))
##    hold on
##    plot(10*log10(cf),p-p2)
##    hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
axes(ax)
cla
xticks([20 50 100 200 400 800 1600 3200 6400 12000 20000])
xticklabels({'20Hz','50Hz','100Hz','200Hz','400Hz','800Hz','1.6kHz','3.2kHz','6.4kHz','12kHz','20kHz'})
xlim('auto')
ylim('auto')
xlim([20 22050])
xlabel('frequency (Hz)')
ylabel('dB')

IR1=IR1/max(abs(IR1))*0.95;
s1=20*log10(abs(fftshift(fft(IR1)))); % compute the GEQ spectrum
n2=NbF/2+2:NbF;
semilogx(ax,nuX(n2),s1(n2));
%semilogx(nuX(n2),s1(n2));

  n800=(abs(nuX)<800);
  nuX2=nuX(n800);
  s12=s1(n800);
  [pks idx] = max(s12);
  f_mx=abs(nuX2(idx))


%[pks idx] = max(s1);
%f_mx=abs(nuX(idx));
plot(ax,xlim,[0 0],'Color','k','linewidth',1,'linestyle',':')
plot(ax,f_mx,pks,'or')
%plot(f_mx,pks,'or')
text(f_mx+50,pks,[num2str(f_mx,'%3.0f'),'Hz']);

mot=['IR spectrum, feedback frequency is: ',num2str(f_mx,'%3.0f'),'Hz'];
title(mot)

  handles.IR1=IR1; % GEQ IR
  handles.IR2=IR2; % Raw IR
  guidata(hObject,handles);
  %waitbar(1, h);
  axes(ax2)
  %cla
	rectangle('Position',[0,0,(round(1000*(1)))+1,20],'FaceColor','b');
  text(480,10,[num2str(round(100*(1))),'%']);
  pause(1)
  %close(h)
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Listen to Mic
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Mic(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,ao,uilisten)

  if (get(hObject,'value')&&(get(uilisten,'value')==0))   % ?gal ? 1 => true
    set(uilisten,'value',1);
    handles = guidata(hObject);
    y=handles.y;
    fs=handles.fs;
    nbits=handles.nbits;
    Nb=size(y,1);
    n3=10*fs:20*fs;
    indxinput=get(icfg,'value');
    if indxinput==1
      MIC=y(:,2); %mic;
      PIC=y(:,1); %pic;
    else
      MIC=y(:,1); %mic;
      PIC=y(:,2); %pic;
    endif


    indx=get(ao,'value');
    id=devinfo.output(indx).ID;

    player = audioplayer (MIC(n3), fs, nbits, id);
    play(player);
    %soundsc(MIC(n3),fs,nbits);
    pause(10.5)
    clear player;
    set(hObject,'value',0);
    set(uilisten,'value',0);
  elseif (get(hObject,'value')&&(get(uilisten,'value')==1))
    set(hObject,'value',0);
    else
    set(hObject,'value',1);
  endif;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Listen to Pickup
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Pickup(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,ao,uilisten)

if (get(hObject,'value')&&(get(uilisten,'value')==0))   % ?gal ? 1 => true
    set(uilisten,'value',1);
    handles = guidata(hObject);
    y=handles.y;
    fs=handles.fs;
    nbits=handles.nbits;
    Nb=size(y,1);
    n3=10*fs:20*fs;
    indxinput=get(icfg,'value');
    if indxinput==1
      MIC=y(:,2); %mic;
      PIC=y(:,1); %pic;
    else
      MIC=y(:,1); %mic;
      PIC=y(:,2); %pic;
    endif
    %soundsc(PIC(n3),fs,nbits);
    indx=get(ao,'value');
    id=devinfo.output(indx).ID;
    player = audioplayer (PIC(n3), fs, nbits, id);
    play(player);
    pause(10.5)
    clear player;
    set(hObject,'value',0);
    set(uilisten,'value',0);
  elseif (get(hObject,'value')&&(get(uilisten,'value')==1))
    set(hObject,'value',0);
  else
    set(hObject,'value',1);
  endif;

end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Listen to IR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function IR(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,ao,uilisten)
if (get(hObject,'value')&&(get(uilisten,'value')==0))   % ?gal ? 1 => true
    set(uilisten,'value',1);
  handles = guidata(hObject);
  y=handles.y;
  fs=handles.fs;
  nbits=handles.nbits;
  IRGEQ=handles.IR1;
  Nb=size(y,1);
  n3=10*fs:20*fs;
  indxinput=get(icfg,'value');
  if indxinput==1
    MIC=y(:,2); %mic;
    PIC=y(:,1); %pic;
  else
    MIC=y(:,1); %mic;
    PIC=y(:,2); %pic;
  endif
  PS=conv(PIC(n3),IRGEQ,'same');
  PS=PS/max(abs(PS))*max(abs(PIC(n3)));
  %soundsc(PS,fs,nbits);
  indx=get(ao,'value');
  id=devinfo.output(indx).ID;
  player = audioplayer (PS, fs, nbits, id);
  play(player);
  pause(10.5)
  clear player;
    set(hObject,'value',0);
    set(uilisten,'value',0);
  elseif (get(hObject,'value')&&(get(uilisten,'value')==1))
    set(hObject,'value',0);
  else
    set(hObject,'value',1);
  endif;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Save IR file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Save(hObject, eventdata,ax,ai,devinfo,fsfmt,icfg,listfmt,Svfmt,Fnam)

  handles = guidata(hObject);
  rawIR=handles.IR2;
  IRGEQ=handles.IR1;
  fs=handles.fs;
  nbits=handles.nbits;
  Nb=size(IRGEQ,1);


  indxfmt=get(fsfmt,'value');
  if (indxfmt <= 3)
    nbits=16;
  else
    nbits=24;
  endif

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% SMooth the high end%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
FS=fs;
NbF=Nb;
dnuX1=FS/NbF;
nuX1=(-NbF/2:NbF/2-1)*dnuX1;

wi = FS/2/1e4*logspace(1,4,NbF);
w=sort([-wi,0,wi],2);
%plot(w)
%grid on

s=fftshift(fft(IRGEQ));
ss=interp1(nuX1,s,w); % On interpole sur une grille/echelle log en frequence (perte d'info dans le haut du spectre

##figure
##plot(nuX1,abs(s),w,abs(ss))
##grid on
##title('amplitude')
##
##figure
##plot(nuX1,abs(s),w,abs(ss))
##grid on
##title('phase')

snew=interp1(w,ss,nuX1); % On retourne sur une grille lineaire
%snew=interp1(w,ss,nuX1,'spline');

% On vire les NaN
for i=1:length(IRGEQ)
    if isnan(snew(i))
        snew(i)=s(i);
    end;
end;

##s1new=20*log10(abs(snew));
##
##figure(1)
##%hold on
##plot(nuX1,s1,nuX1,s1new);
##grid on
##set(gca,'Xscale','log')
##xlim([20 20e3])
##ylim([-60 40])
##xlabel('Frequency (Hz)')
##ylabel('dB')
##%hold off

yy=real(ifft(ifftshift(snew)));
%IR0old=IR0;
IRITP=yy/max(abs(yy))*0.95;


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Standard
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if ispc
    [filename, pathname] = uiputfile('*.wav', 'Save as','IR_std.wav');
    Fname=[pathname,filename];
    audiowrite(Fname,IRGEQ/max(abs(IRGEQ))*0.95,fs,'BitsPerSample',nbits);

    [filename, pathname] = uiputfile('*.wav', 'Save as','IR_itp.wav');
    Fname=[pathname,filename];
    audiowrite(Fname,IRITP/max(abs(IRITP))*0.95,fs,'BitsPerSample',nbits);
  end
  if ismac
    Fname=get(Fnam,'string');
    audiowrite(Fname,IRGEQ/max(abs(IRGEQ))*0.95,fs,'BitsPerSample',nbits);

    Fname2=[Fname(1:end-4),'_ITP.wav'];
    audiowrite(Fname2,IRITP/max(abs(IRITP))*0.95,fs,'BitsPerSample',nbits);
  end
  Fname2=[Fname(1:end-4),'.jpg'];
  print(gcf,Fname2,'-djpeg','-r300');
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % Options
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  SAVEfmt=get(Svfmt,'value');
  if SAVEfmt==2
    % Raw
    Fname2=[Fname(1:end-4),'_raw.wav'];
    audiowrite(Fname2,rawIR/max(abs(rawIR))*0.95,fs,'BitsPerSample',nbits);
    % GEQ+Blend
    IR0=zeros(Nb,1);
    IR0(1)=1;
    IRblend=(IR0+IRGEQ)/2;
    Fname2=[Fname(1:end-4),'_Bld.wav'];
    audiowrite(Fname2,IRblend/max(abs(IRblend))*0.95,fs,'BitsPerSample',nbits);
    % Minimum Phase
    s=fft(IRGEQ);
    sm = exp( fft( fold( ifft( log( clipdb(s,-100) )))));
    IRmph=real(ifft(sm));
    Fname2=[Fname(1:end-4),'_MPh.wav'];
    audiowrite(Fname2,IRmph/max(abs(IRmph))*0.95,fs,'BitsPerSample',nbits);
  end;
  if SAVEfmt==3 % Feedback friendly
    y=handles.y;
    fs=handles.fs;
    %nbits=handles.nbits;
    Nb=size(y,1);
    FS=fs;
    indxinput=get(icfg,'value');
    if indxinput==1
      MIC=y(:,2); %mic;
      PIC=y(:,1); %pic;
    else
      MIC=y(:,1); %mic;
      PIC=y(:,2); %pic;
    endif
  nn=(10*FS+1:20*FS);
  MS=MIC(nn,1);
  PS=PIC(nn,1);


  % COmpute Feedback frequency
  NbF=Nb;
  dnuX=FS/NbF;
  nuX=(-NbF/2:NbF/2-1)*dnuX;

  s1=20*log10(abs(fftshift(fft(IRGEQ)))); % compute the GEQ spectrum

  n800=(abs(nuX)<800);
  nuX2=nuX(n800);
  s12=s1(n800);
  [pks idx] = max(s12);
  f_mx=abs(nuX2(idx))

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % generate PEQ
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  g=-6; % -6dB ? la resonance
  Q=3;
  fc=f_mx;


  % IIR coefficient calculation
  g=10^(g/20);
  t0=2*pi*fc/fs;
  if g >= 1
      beta=t0/(2*Q);
  else
      beta=t0/(2*g*Q);
  endif
  a2=-0.5*(1-beta)/(1+beta);
  a1=(0.5-a2)*cos(t0);
  b0=(g-1)*(0.25+0.5*a2)+0.5;
  b1=-a1;
  b2=-(g-1)*(0.25+0.5*a2)-a2;

  % SOS Form
  b=2*[b0 b1 b2];
  a=[1 -2*a1 -2*a2];
  IRX=filter(b,a,IRGEQ); % FIltre IR

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  % generate Shelf EQ
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    g=-1; % -1dB en shelf
    fc=1000;

    g=10^(g/20);
    if g >= 1
        a1=-0.5*tan(pi*(fc/fs-0.25));
    else
        a1=-0.5*tan(pi*(fc/fs*g-0.25));
    end
    b0t=0.25+0.5*a1;
    b2=0;
    a2=0;
    b0=0.5+(g-1)*b0t;
    b1=-a1-(g-1)*b0t;

    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    IRX1=filter(b,a,IRX); % FIltre IR
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    g=-2; % -1dB en shelf
    fc=1000;

    g=10^(g/20);
    if g >= 1
        a1=-0.5*tan(pi*(fc/fs-0.25));
    else
        a1=-0.5*tan(pi*(fc/fs*g-0.25));
    end
    b0t=0.25+0.5*a1;
    b2=0;
    a2=0;
    b0=0.5+(g-1)*b0t;
    b1=-a1-(g-1)*b0t;

    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    IRX2=filter(b,a,IRX); % FIltre IR
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    g=-3; % -1dB en shelf
    fc=1000;

    g=10^(g/20);
    if g >= 1
        a1=-0.5*tan(pi*(fc/fs-0.25));
    else
        a1=-0.5*tan(pi*(fc/fs*g-0.25));
    end
    b0t=0.25+0.5*a1;
    b2=0;
    a2=0;
    b0=0.5+(g-1)*b0t;
    b1=-a1-(g-1)*b0t;

    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    IRX3=filter(b,a,IRX); % FIltre IR
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    g=-4; % -1dB en shelf
    fc=1000;

    g=10^(g/20);
    if g >= 1
        a1=-0.5*tan(pi*(fc/fs-0.25));
    else
        a1=-0.5*tan(pi*(fc/fs*g-0.25));
    end
    b0t=0.25+0.5*a1;
    b2=0;
    a2=0;
    b0=0.5+(g-1)*b0t;
    b1=-a1-(g-1)*b0t;

    % SOS Form
    b=2*[b0 b1 b2];
    a=[1 -2*a1 -2*a2];
    IRX4=filter(b,a,IRX); % FIltre IR
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Compute brightness Boost
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Generate Octave spectrum (Pic convolved & Pic convolved + EQ)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %pkg load signal;
    PS0=conv(PIC(nn,1),IRGEQ,'same');
    PS1=conv(PIC(nn,1),IRX1,'same');
    PS2=conv(PIC(nn,1),IRX2,'same');
    PS3=conv(PIC(nn,1),IRX3,'same');
    PS4=conv(PIC(nn,1),IRX4,'same');

    [p,cf,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS0/max(abs(PS0)),FS,3,1,1,0);
    [p21,cf21,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS1/max(abs(PS1)),FS,3,1,1,0);
    [p22,cf22,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS2/max(abs(PS2)),FS,3,1,1,0);
    [p23,cf23,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS3/max(abs(PS3)),FS,3,1,1,0);
    [p24,cf24,overall_lev,overall_levA,sfilt] = oct_spectrum2(PS4/max(abs(PS4)),FS,3,1,1,0);
    OSC0=sum(cf.*abs(p))/sum(abs(p));
    OSC1=sum(cf21.*abs(p21))/sum(abs(p21));
    OSC2=sum(cf22.*abs(p22))/sum(abs(p22));
    OSC3=sum(cf23.*abs(p23))/sum(abs(p23));
    OSC4=sum(cf24.*abs(p24))/sum(abs(p24));

    osc0=[OSC0 OSC0 OSC0 OSC0];
    osc=[OSC1 OSC2 OSC3 OSC4];
    [w, iw] = min(abs(osc-osc0));
    Fname2=[Fname(1:end-4),'_Fbk.wav'];
    if (iw==1)
      audiowrite(Fname2,IRX1/max(abs(IRX1))*0.95,fs,'BitsPerSample',nbits);
    end
    if (iw ==2)
      audiowrite(Fname2,IRX2/max(abs(IRX2))*0.95,fs,'BitsPerSample',nbits);
    end
    if (iw ==3)
      audiowrite(Fname2,IRX3/max(abs(IRX3))*0.95,fs,'BitsPerSample',nbits);
    end
    if (iw ==4)
      audiowrite(Fname2,IRX4/max(abs(IRX4))*0.95,fs,'BitsPerSample',nbits);
    end



end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% JF45 BEGINING OF Jon Fields Algorithm Copyright Jon Fields
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  if SAVEfmt==4 % JF45
    y=handles.y;
    fs=handles.fs;
    %nbits=handles.nbits;
    Nb=size(y,1);
    FS=fs;
    indxinput=get(icfg,'value');
    if indxinput==1
      MIC=y(:,2); %mic;
      PIC=y(:,1); %pic;
    else
      MIC=y(:,1); %mic;
      PIC=y(:,2); %pic;
    endif

    MS=MIC(:,1);
    PS=PIC(:,1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    # IR generation:
      ss = (2^17); # ~3 second segment size for analysis
      nstep = 2; # segment step ~6 seconds
      nstart = 3; # index into .wav in segments, ~6 seconds
      nlast = 17; # up to 8 segments checked
      count = 0; # initialize good segment count
      nzcount = 0; # initialize near zero count
      accirfftnnz = zeros(ss,1); # initialize no near zero IR
      # Blackman window function:
      window = (.42-.5*cos(2*pi*(0:ss-1)/(ss-1))+.08*cos(4*pi*(0:ss-1)/(ss-1)))';
      for n = nstart:nstep:nlast # process segments
        start = (((n-1) * ss) + 1);
        finish = (n * ss);
        # load segment from wav file
        %[s,fs] = audioread([input,'.wav'],[start, finish]);
        s=y(start:finish,:);
        smax = max(max(s));
        clip = (smax > 0.999); # check for clipping
        toolow = (smax < 0.178); # 15 dB down
        # calculate per segment IR
        if (clip == 0 && toolow == 0 && count < 4)
          pickup = PS(start:finish) .* window;
          mic = MS(start:finish) .* window;
          pupfft = fft(pickup);
          micfft = fft(mic);
          pupfftnnz = pupfft;
          micfftnnz = micfft;
          nearzero = 10^(-65/20)*abs(max(pupfft)); # -65dB
          for m = 1:1:ss
            # check and fix  near zeros in pickup FFT
            if (abs(pupfft(m)) < nearzero)
              nzcount = nzcount + 1; # Count near 0s
              pupfftnnz(m) = 1; # erase near zero
              micfftnnz(m) = 1; # erase near zero
            end
          end
          # IR=FFT(mic)/FFT(PUP)
          irfftnnz = micfftnnz ./ pupfftnnz;
          # accumulate IRs
          accirfftnnz = accirfftnnz + irfftnnz;
          count = count + 1; # increment segment count
        end
      end
      if(count==0)
      ['Zero Segments Processed due to Clip/Min errors']
      return;
    end
    irnnz = ifft(accirfftnnz/count); # calc IR
    ir2048nnz = irnnz(1:2048); # truncate to 2048
    avgnzcount = nzcount / (count*2);
%    ['Processed Segments = ',num2str(count)]
%    ['Average Near Zero per Segment = ',num2str(avgnzcount)]
%    filename = ['jf45ir',input,'.wav'];
%    audiowrite(filename,ir2048nnz,fs); # write IR
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% END OF Jon Fields Algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Fname2=[Fname(1:end-4),'_JF45.wav'];
    audiowrite(Fname2,ir2048nnz/max(abs(ir2048nnz))*0.95,fs,'BitsPerSample',nbits);


  end;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Close app
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Close(hObject, eventdata,ax,dlg)
  close(dlg)
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Donate
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Donate(hObject, eventdata,ax,dlg)
  if ispc
    system('start "" "https://www.paypal.com/cgi-bin/webscr?cmd=_donations&business=kien.phanhuy%40free.fr&currency_code=USD&source=url"');
  end;
  if ismac
    system('open "" "https://www.paypal.com/cgi-bin/webscr?cmd=_donations&business=kien.phanhuy%40free.fr&currency_code=USD&source=url"');
  end;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Sub-function from J. O Smith
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [rw] = fold(r)
% [rw] = fold(r)
% Fold left wing of vector in "FFT buffer format"
% onto right wing
% J.O. Smith, 1982-2002

   [m,n] = size(r);
   if m*n ~= m+n-1
     error('fold.m: input must be a vector');
   end
   flipped = 0;
   if (m > n)
     n = m;
     r = r.';
     flipped = 1;
   end
   if n < 3, rw = r; return;
   elseif mod(n,2)==1
       nt = (n+1)/2;
       rw = [ r(1), r(2:nt) + conj(r(n:-1:nt+1)), ...
             0*ones(1,n-nt) ];
   else
       nt = n/2;
       rf = [r(2:nt),0];
       rf = rf + conj(r(n:-1:nt+1));
       rw = [ r(1) , rf , 0*ones(1,n-nt-1) ];
   end

   if flipped
     rw = rw.';
   end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
   function [clipped] = clipdb(s,cutoff)
% [clipped] = clipdb(s,cutoff)

% Clip magnitude of s at its maximum + cutoff in dB.
% Example: clip(s,-100) makes sure the minimum magnitude
% of s is not more than 100dB below its maximum magnitude.
% If s is zero, nothing is done.

clipped = s;
as = abs(s);
mas = max(as(:));
if mas==0, return; end
if cutoff >= 0, return; end
thresh = mas*10^(cutoff/20); % db to linear
toosmall = find(as < thresh);
clipped = s;
clipped(toosmall) = thresh;
end

%% OCTAVE SPECTRUM
%   Computation of octave spectrum and its overall level with dB-A weighting
%   In compliance with specification for octave-band and fractional-octave band
%   analog and digital filters ANSI S1.11-2004
%
%   INPUTS
%   s = input signal
%   fs = sampling frequency
%   b = octave ratio reccomended 1 or 3
%   dbRef = standard reference for decibel scale calculation, default value: 1
%   weightFlag = A-weighting [0,1], default value: 1
%   plotFlag = generate octave diagram [0,1], default value: 1
%
%   INPUTS
%   S = octave spectrum (dB)
%   fm = midband frequency
%   overall_lev = overall level
%   overall_levA = weighted overall level
%   sfilt = third-octave filtered signals
%
% M. Buzzoni
% Dec. 2018
% Mod by Kien Phan Huy
function [S,fm,overall_lev,overall_levA,sfilt] = oct_spectrum2(s,fs,b,dbRef,weightFlag,plotFlag);

if nargin < 3
  b = 3;
  dbRef = 1;
  weightFlag = 1;
  plotFlag = 1;
end

if nargin < 4
  dbRef = 1;
  weightFlag = 1;
  plotFlag = 1;
end

if nargin < 5
  weightFlag = 1;
  plotFlag = 1;
end

if nargin < 6
  weightFlag = 1;
end

s = s(:)';
L = length(s);
G = 10^(3/10); % octave ratio, base 10
fr = 1000; % reference frequency
x = 1:43; % band number
Gstar = G.^((x - 30)/b);
fm = Gstar.*fr; % midband frequency
f2 = G.^(+1/(2*b)).*fm;  % upper bandedge freq.

ind2cut = or(f2 < 20, f2 > fs/2); % check for exceeding frequency bins
fm(ind2cut) = [];
f2(ind2cut) = [];
x(ind2cut) = [];
f1 = G.^(-1/(2*b)).*fm; % lower bandedge freq.

Br = (f2 - f1)./fm; % normalized reference bandwidth
sfilt = zeros(length(x),L); % initialization of third-octave filtered signals

overall_lev = 10.*log10((rms2(s)./dbRef).^2); % overall level

if weightFlag == 1
% weighting function
ff = (0:L-1)*fs/L;
ff = ff(1:(floor(L/2)+1));
num = (3.5041384*10^16) * ff.^8;
den = (20.598997^2 + ff.^2).^2 .* (107.65265^2 + ff.^2) .* (737.86223^2 + ff.^2) .* (12194.217^2 + ff.^2).^2;
alphaA = num./den; alphaA = [alphaA fliplr(alphaA(2:end - 1))];
s = real(ifft(fft(s) .* alphaA)); % A-weighted signal
end

% filtering
for k = 1:length(x)
  %[B, A] = butter (2,[f1(k) f2(k)]./(fs/2));
  [B,A]= bp_synth(2,(f1(k)+f2(k))/2,abs(f1(k)-f2(k)),fs);
  sfilt(k,:) = filter(B,A,s);
end

% octave spectrum
S = 10.*log10((rms2(sfilt') ./dbRef).^ 2); % band levels in dB or dBA

% weighted overall level
if weightFlag == 1
  overall_levA = 10*log10(sum(10.^(S/10)));
else
  overall_levA = 0;
end

% plot
Smin = min(S);
if plotFlag == 1
  figure('Position',[40 40 1400 500])
  h = bar([S Smin - 5 Smin - 5],'b','BaseValue',Smin - 3,'barwidth',.3);

  if weightFlag == 1
    hold on, bar(length(S)+1,overall_levA,'g','BaseValue',Smin - 3,'barwidth',.3)
  end

  hold on, bar(length(S)+2,overall_lev,'r','BaseValue',Smin - 3,'barwidth',.3)
  ylim([Smin - 3 max([overall_lev overall_levA S]) .* 1.1])
  xlim([0 length(S)+4])

  if b == 3
    title(['third-octave spectrum'])
  elseif b == 1
    title(['octave spectrum'])
  else
    title(['1/' num2str(b) ' octave spectrum'])
  end

  xlabel('midband frequency (Hz)')

  if weightFlag == 1
    ylabel('dBA')
  else
    ylabel('dB')
  end

  xticks(1:length(S) + 2)
  xticklabels(round(fm))
  box off
  set (gca, 'ygrid', 'on');

  if weightFlag == 1
    hleg = legend('octave levels','overall (A)','overall','location','eastoutside');
  else
    hleg = legend('octave levels','overall','location','eastoutside');
  end

  yticks(floor(Smin - 3 ):3:(max([overall_lev overall_levA S]) .* 1.1))
end
endfunction

% RMS function
function rms=rms2(x)
  N=length(x);
  rms=sqrt(1/N*sum(x.^2));
endfunction

%function [b,a]= bp_synth(N,fcenter,bw,fs)      12/29/17 neil robertson
% Synthesize IIR Butterworth Bandpass Filters
%
% N= order of prototype LPF
% fcenter= center frequency, Hz
% bw= -3 dB bandwidth, Hz
% fs= sample frequency, Hz
% [b,a]= vector of filter coefficients
% https://www.dsprelated.com/showarticle/1128.php
function [b,a]= bp_synth(N,fcenter,bw,fs)
f1= fcenter- bw/2;            % Hz lower -3 dB frequency
f2= fcenter+ bw/2;            % Hz upper -3 dB frequency
if f2>=fs/2;
   error('fcenter+ bw/2 must be less than fs/2')
end
if f1<=0
   error('fcenter- bw/2 must be greater than 0')
end
% find poles of butterworth lpf with Wc = 1 rad/s
k= 1:N;
theta= (2*k -1)*pi/(2*N);
p_lp= -sin(theta) + j*cos(theta);
% pre-warp f0, f1, and f2
F1= fs/pi * tan(pi*f1/fs);
F2= fs/pi * tan(pi*f2/fs);
BW_hz= F2-F1;              % Hz prewarped -3 dB bandwidth
F0= sqrt(F1*F2);           % Hz geometric mean frequency
% transform poles for bpf centered at W0
% note:  alpha and beta are vectors of length N; pa is a vector of length 2N
alpha= BW_hz/F0 * 1/2*p_lp;
beta= sqrt(1- (BW_hz/F0*p_lp/2).^2);
pa= 2*pi*F0*[alpha+j*beta  alpha-j*beta];
% find poles and zeros of digital filter
p= (1 + pa/(2*fs))./(1 - pa/(2*fs));      % bilinear transform
q= [-ones(1,N) ones(1,N)];                % N zeros at z= -1 (f= fs/2) and z= 1 (f = 0)
% convert poles and zeros to numerator and denominator polynomials
a= poly(p);
a= real(a);
b= poly(q);
% scale coeffs so that amplitude is 1.0 at f0
f0= sqrt(f1*f2);
h= freqz(b,a,[f0 f0],fs);
K= 1/abs(h(1));
b= K*b;
endfunction




##Copyright (c) 2019, Marco Buzzoni
##All rights reserved.
##
##Redistribution and use in source and binary forms, with or without
##modification, are permitted provided that the following conditions are met:
##
##* Redistributions of source code must retain the above copyright notice, this
##  list of conditions and the following disclaimer.
##
##* Redistributions in binary form must reproduce the above copyright notice,
##  this list of conditions and the following disclaimer in the documentation
##  and/or other materials provided with the distribution
##* Neither the name of Universit? degli Studi di Ferrara nor the names of its
##  contributors may be used to endorse or promote products derived from this
##  software without specific prior written permission.
##THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
##AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
##IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
##DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
##FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
##DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
##SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
##CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
##OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
##OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
##
##
##``Introduction to Digital Filters with Audio Applications'', by Julius O. Smith III, (September 2007 Edition).
##Copyright ? 2020-05-11 by Julius O. Smith III
##Center for Computer Research in Music and Acoustics (CCRMA),   Stanford University
##
##Copyright (c) 2014, Yuanfei
##All rights reserved.
##
##Redistribution and use in source and binary forms, with or without
##modification, are permitted provided that the following conditions are met:
##
##* Redistributions of source code must retain the above copyright notice, this
##  list of conditions and the following disclaimer.
##
##* Redistributions in binary form must reproduce the above copyright notice,
##  this list of conditions and the following disclaimer in the documentation
##  and/or other materials provided with the distribution
##* Neither the name of  nor the names of its
##  contributors may be used to endorse or promote products derived from this
##  software without specific prior written permission.
##THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
##AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
##IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
##DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
##FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
##DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
##SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
##CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
##OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
##OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.