function lpcspect(action)
%lpcspect Demonstration of LPC spectra, poles and LPC error

% Martin Cooke, April 1999



if nargin < 1
  action='init';
  % material below is for recording user actions
  % else
  % ud=get(gcf,'UserData');
  % ud.record=strvcat(ud.record,[datestr(now,13) ' ' action]);
  % set(gcf,'UserData',ud);
end

switch action
case 'init'
  
  f=findobj('Tag','lpcspect_fig');
  if ~isempty(f)
    figure(f);
  else     
    lpcspect_gui;
    ud.record=[];
    ud.linkCursors=get(findobj('Tag','linkCursorsCheck'),'Value');
    ud.useHamming=get(findobj('Tag','useHammingCheck'),'Value'); 
    ud.fs=22050;
    ud.neverloaded=1;
    ud.lpcN=getuicontrolvalue(ud.lpcordermenu,'num');
    set(findobj(gcf,'Style','checkbox'),'Enable','off');
    set(findobj(gcf,'Style','pushbutton'),'Enable','off');
    set(gcf,'UserData',ud,'Name','lpcspect - [no data loaded]');
  end

case 'create'
  [sig,fs,name]=createsig;
  if length(sig)
    ud=get(gcf,'UserData');
    ud.sig=sig; ud.name=name; ud.fs=fs/1000;
    set(gcf,'UserData',ud);
    lpcspect 'newsig'
  end
  
case 'load'
  ud=get(gcf,'UserData');
  [name,sig,fs]=loadsig(fullfile(madroot,'data','sounds'));
  if length(sig) 
    ud.sig=sig; ud.name=name; ud.fs=fs/1000;
    set(gcf,'UserData',ud);
    lpcspect 'newsig'
  end

case 'newsig'
  ud=get(gcf,'UserData');
     if ud.neverloaded
       ud.neverloaded=0;
       set(findobj(gcf,'Style','checkbox'),'Enable','on'); 
       set(findobj(gcf,'Style','pushbutton'),'Enable','on');
       set(findobj(gcf,'Style','popupmenu'),'Enable','on');
     else
       delete([ud.r ud.l ud.seg ud.dft ud.lpcspectrum ud.errordft ud.sigline ud.poles]);
     end  
    
    ud.maxsig=max(abs(ud.sig(:)));
    if get(findobj('Tag','preemphCheck'),'Value')
       ud.sig=preemph(ud.sig);
    end
    ud.sig = ud.sig/ud.maxsig;
    ud.sig=ud.sig(:)';
    ud.samps=length(ud.sig);
    
    ud.winleft=1;
    ud.winright=ud.samps;
    set(ud.wavAxes,'NextPlot','replacechildren','Drawmode','fast','ButtonDownFcn','lpcspect play',...
       'XLim',[0 length(ud.sig)],'YLim',[-1.2 1.2]);
    axes(ud.wavAxes);
    ud.left=round(0.1*length(ud.sig));
    ud.right=min(ud.left+512,round(0.9*length(ud.sig)));   
    ud.l=line([ud.left ud.left],get(gca,'YLim'),'Color','r',...
      'LineWidth',2,'ButtonDownFcn','lpcspect left','EraseMode','xor','LineStyle','--');
    hold on;
    ud.r=line([ud.right ud.right],get(gca,'YLim'),'Color','r',...
      'LineWidth',2,'ButtonDownFcn','lpcspect right','EraseMode','xor','LineStyle','--');
    ud.sigline=plot(ud.sig);
    set(ud.sigline,'EraseMode','background' ,'ButtonDownFcn','lpcspect play');
    ud.freqs=0:0.010:(ud.fs/2);   % 20 Hz spacing 
    
    axes(ud.wavsegAxes);
    set(ud.wavsegAxes,'NextPlot','replacechildren','Drawmode','fast','YLim',[-1.7 1.2],'YTick',[],'XTick',[]);
    ud.wavseg=ud.sig(ud.left:ud.right);     
    ud.seg=plot(ud.wavseg);
    hold on
    set(ud.seg,'EraseMode','background','ButtonDownFcn','lpcspect playseg');
    ud.errorsig=zeros(size(ud.wavseg));
    ud.err=plot(ud.errorsig);
    set(ud.err,'EraseMode','background','Color','g','ButtonDownFcn','lpcspect playerror');

    % LPC and DFT spectra axes
    spect=performFFT(ud.wavseg,ud.freqs,ud.fs);    
    axes(ud.spectAxes);
    ud.dft=plot(ud.freqs,spect);
    hold on
    ud.lpcspectrum=plot(ud.freqs,spect);
    set(ud.dft,'EraseMode','background'); 
    set(ud.lpcspectrum,'EraseMode','background','Color','g','LineWidth',2); 
    set(ud.spectAxes,'Drawmode','fast','XLim',[0 ud.fs/2],'YLim',[-80 50],'YTick',[]);

    errordft=performFFT(ud.sig(ud.left:ud.right),ud.freqs,ud.fs);    
    ud.errordft=plot(ud.freqs,errordft);
    set(ud.errordft,'EraseMode','background','Color','g'); 

    xlabel('frequency (kHz)');
    ylabel('log magnitude (dB)');
    
    % pole-zero (no zeros!) axes
    axes(ud.PZAxes);
    set(ud.PZAxes,'XLim',[-1.1 1.1],'YLim',[-1.1 1.1],'ButtonDownFcn','lpcspect showpole');
      
    % define unit circle
    th = 0:pi/50:2*pi;
    xunit = cos(th);
    yunit = sin(th);
    uc=plot(xunit,yunit); hold on
    xa=plot([-1.1 1.1],[0 0]);  
    ya=plot([0 0],[-1.1 1.1]);
    set([uc xa ya],'ButtonDownFcn','lpcspect showpole','LineStyle','--','color','r');
    ud.poles=line([0 0],[0 0]);
    set(ud.poles,'EraseMode','background','Marker','x','LineStyle','none','Color','g','MarkerSize',13,...
       'ButtonDownFcn','lpcspect showpole');
    
    set(gcf,'UserData',ud,'Name',['lpcspect - ' ud.name]);
    lpcspect 'linkCursors'
    lpcspect 'recalc'
  

case 'clearpole'
  ud=get(gcf,'UserData');
  set(ud.freqstring,'String','');
  set(gcf,'UserData',ud);
  
case 'showpole'
  ud=get(gcf,'UserData');
  currentPoint=get(gca,'CurrentPoint');
  xpos=currentPoint(1);
  ypos=currentPoint(3);  
  %logz=log(xpos+i*ypos);
  %freq=imag(logz)*(ud.fs*1000)/(2*pi);
  freq=atan2(abs(ypos),xpos)*(ud.fs*1000)/(2*pi);
  if ypos < 0
    freq=(ud.fs*1000) - freq;
  end
  %bw=real(1/logz)*(ud.fs*1000)/pi;
  set(ud.freqstring,'String',['freq=' num2str(round(freq))],'Enable','on');
  set(gcf ,'WindowButtonUpFcn','lpcspect clearpole');
  set(gcf,'UserData',ud);
  
case 'lpcorderchanged'
  ud=get(gcf,'UserData');
  ud.lpcN=getuicontrolvalue(ud.lpcordermenu,'num');
  set(gcf,'UserData',ud);
  lpcspect 'recalc'
  
case 'linkCursors'
  ud=get(gcf,'UserData');
  ud.linkCursors=get(findobj('Tag','linkCursorsCheck'),'Value');
  if ud.linkCursors
    ud.gap=ud.right-ud.left+1;
  end
  set(gcf,'UserData',ud);


case 'useHamming'
  ud=get(gcf,'UserData');
  ud.useHamming=get(findobj('Tag','useHammingCheck'),'Value'); 
  set(gcf,'UserData',ud);
  lpcspect 'recalc'

case 'unzoom'
  ud=get(gcf,'UserData');
  mid=round((ud.winright-ud.winleft)/2);
  ud.winleft=max(1,ud.winleft-mid);
  ud.winright=min(ud.samps,ud.winright+mid); 	 
  set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  set(gcf,'UserData',ud);
  
case 'showall'
  ud=get(gcf,'UserData');
  ud.winleft=1;
  ud.winright=ud.samps; 	 
  set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  set(gcf,'UserData',ud);
  
case 'zoomToCursors' 
  ud=get(gcf,'UserData'); 
  if ud.linkCursors
    ud.winleft=max(1,round(ud.left-0.1*(ud.right-ud.left)));
    ud.winright = min(ud.samps,round(ud.right+0.1*(ud.right-ud.left)));
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
    set(gcf,'UserData',ud);
  else	 
    ud.winleft = ud.left;
    ud.winright = ud.right;
    ud.left = round(ud.winleft+0.1*(ud.winright-ud.winleft));
    ud.right = round(ud.winright-0.1*(ud.winright-ud.winleft));
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
    set(ud.l,'XData',[ud.left ud.left]);
    set(ud.r,'XData',[ud.right ud.right]);   
    set(gcf,'UserData',ud);
    lpcspect 'recalc'
  end 

case 'zoom'
  ud=get(gcf,'UserData');
  if ud.right > ud.left+5
    
    cmid=round((ud.left+ud.right)/2);
    if ud.linkCursors
      ud.winleft = min(ud.left-5,round((ud.winleft+cmid)/2));
      ud.winright = max(ud.right+5,round((ud.winright+cmid)/2));
      set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
      set(gcf,'UserData',ud);
    else
      ud.winleft = round((ud.winleft+cmid)/2);
      ud.winright = round((ud.winright+cmid)/2);
      set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
      redo=0;
      if ud.left < ud.winleft+5
        ud.left = ud.winleft+5;
        set(ud.l,'XData',[ud.left ud.left]);
        redo=1;
      end
      if ud.right > ud.winright-5
        ud.right = ud.winright-5;
        set(ud.r,'XData',[ud.right ud.right]);
        redo=1;
      end
      set(gcf,'UserData',ud);
      if redo
        lpcspect 'recalc'
      end	 
    end	
  end 

case 'left'
  set(gcf,'WindowButtonMotionFcn','lpcspect moveleft');
  set(gcf,'WindowButtonUpFcn','lpcspect release');
  
case 'right'
  set(gcf,'WindowButtonMotionFcn','lpcspect moveright');
  set(gcf,'WindowButtonUpFcn','lpcspect release');
  
case 'moveleft'
  ud=get(gcf,'UserData');
  currentPoint=get(gca,'CurrentPoint');
  ud.left=min(ud.right-1,max(1,ceil(currentPoint(1))));
  if ud.linkCursors
    step=50;
  else	
    step=round(0.2*(ud.winright-ud.winleft));
  end	 
  if ud.left <= ud.winleft
    ud.winleft=max(1,ud.winleft-step);
    if ud.winleft > 1
      ud.winright=ud.winright-step;
    end	  
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  end
  if ud.linkCursors
    ud.right = ud.left + ud.gap;
    if ud.right >= ud.samps
      ud.right=ud.samps;
      ud.left=ud.samps-ud.gap;
    end
  end	  
  if ud.right >= ud.winright
    ud.winright=min(ud.samps,ud.winright+step);
    if ud.winright < ud.samps
      ud.winleft=ud.winleft+step;
    end	  
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  end    

  set(ud.l,'XData',[ud.left ud.left]);
  set(ud.r,'XData',[ud.right ud.right]);
  set(ud.gapstring,'String',['selection: ' num2str(round((ud.right-ud.left+1)/ud.fs)) 'ms']);
  set(gcf,'UserData',ud);

case 'moveright'
  ud=get(gcf,'UserData');
  currentPoint=get(gca,'CurrentPoint');
  ud.right=min(ud.samps,max(ud.left+1,ceil(currentPoint(1))));
  if ud.linkCursors
    step=50;
  else	
    step=round(0.2*(ud.winright-ud.winleft));
  end	
  if ud.right >= ud.winright
    ud.winright=min(ud.samps,ud.winright+step);
    if ud.winright < ud.samps
      ud.winleft=ud.winleft+step;
    end	  
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  end  
  if ud.linkCursors
    ud.left=ud.right-ud.gap;
    if ud.left < 0
      ud.left=1;
      ud.right=ud.left+ud.gap-1;
    end  
  end	
  if ud.left <= ud.winleft
    ud.winleft=max(1,ud.winleft-step);
    if ud.winleft > 1
      ud.winright=ud.winright-step;
    end	  
    set(ud.wavAxes,'XLim',[ud.winleft ud.winright]);
  end 
  set(ud.l,'XData',[ud.left ud.left]);
  set(ud.r,'XData',[ud.right ud.right]);  
  set(ud.gapstring,'String',['selection: ' num2str(round((ud.right-ud.left+1)/ud.fs)) 'ms']);
  set(gcf,'UserData',ud);

case 'recalc'
  ud=get(gcf,'UserData');
  % recompute and display lpc spectrum, error signal and pole plot
  if ud.useHamming
     seg=window(ud.right-ud.left+1,'hamming')'.*ud.sig(ud.left:ud.right);
  else
     seg=ud.sig(ud.left:ud.right);
  end

  set(ud.wavsegAxes,'XLim',[1 length(seg)]);
  set(ud.seg,'YData',seg,'XData',1:length(seg));  
  spect=performFFT(seg,ud.freqs,ud.fs);                   % DFT
  [a,g]=lpc(seg,ud.lpcN);                                        % LPC
  lpcspectrum = 20*log10(freqz(g,a,ud.freqs,ud.fs));  % LPC spectrum
  
  % compute error signal
  a(1)=0;
  f=filter(-a,1,seg);
  err=seg-f;
  set(ud.err,'YData',err-1,'XData',1:length(err)); 
  ud.errsig=err;
  
  % compute MSE as a percentage of signal MSE
  sigmse=sqrt(mean(seg.^2));
  errmse=sqrt(mean(err.^2));
  set(ud.totalerror,'String',['mse: ' num2str(100*errmse/sigmse) '%']);
  errdft=performFFT(err,ud.freqs,ud.fs);
  a(1)=1; poles=roots(a);
  set(ud.poles,'XData',real(poles),'YData',imag(poles));
  
  % display them
  set(ud.lpcspectrum,'YData',lpcspectrum);
  set(ud.dft,'YData',spect);
  set(ud.errordft,'YData',errdft-30);
  set(gcf,'UserData',ud);
  
case 'release'
  ud=get(gcf,'UserData');
  set(gcf,'WindowButtonMotionFcn','');
  set(gcf,'WindowButtonUpFcn','');
  if  get(findobj('Tag','playCheck'),'Value'); 
    ud=get(gcf,'UserData');
    sound(ud.sig(ud.left:ud.right),ud.fs*1000);
  end	
  lpcspect 'recalc'

case 'preemphChecked'
  ud=get(gcf,'UserData');
  %ud.preemph=get(findobj('Tag','premphCheck'),'Value')
  ud.sig=ud.maxsig*ud.sig;
  if get(findobj('Tag','preemphCheck'),'Value')
     ud.sig=preemph(ud.sig);
  else
     ud.sig=deemph(ud.sig);
  end
  ud.maxsig=max(abs(ud.sig));
  ud.sig = ud.sig/ud.maxsig;
  set(ud.sigline,'YData',ud.sig);
  set(gcf,'UserData',ud);
  lpcspect 'recalc'
  
case 'play'
  ud=get(gcf,'UserData');
  currentPoint=get(gca,'CurrentPoint');
  pos=currentPoint(1);
  if and(pos>ud.left,pos<ud.right)
    soundsc(ud.sig(ud.left:ud.right),ud.fs*1000);
  else
    soundsc(ud.sig(ud.winleft:ud.winright),ud.fs*1000);
  end

case 'playerror'
  ud=get(gcf,'UserData');
  soundsc(ud.errsig,ud.fs*1000);
  
case 'playseg'
  ud=get(gcf,'UserData');
  soundsc(ud.sig(ud.left:ud.right),ud.fs*1000);
   
otherwise
  maduievent('lpcspect',action);
  
end

function s=performFFT(sig,fvals,fs)
  % compute FFT of signal and return result at specified frequencies
  lsig=length(sig);
  lf=length(fvals);
  % find nearest power of 2 to signal length
  f=fft(sig,max(64,power(2,round(log2(lsig)))));
  f=20*log10(1e-10+abs(f(2:round(length(f)/2))));
  s=interp1(linspace(0,fs/2,length(f)),f,fvals);