picky.m

## Picture Kymophone 1.0
## Copyright (C) 2016 Patrick Feaster
## This program is free software and may be used for any purpose; however, the copyright 
## notice must be maintained.  Patrick Feaster is not responsible for any consequences
## of using this software, which is distributed with NO WARRANTY OF ANY KIND.
function picky(output="w",Emode="0",Ipower="1",Fmat="d",corr1="1",Pc="1",wavscale="0.95",transduction="all",samplerate="44100",bitdepth="24",sourcefileinput=0);
if(columns(transduction)!=3);#Exit function if value entered for transduction is of wrong length
  printf(strcat(transduction,": Invalid value for transduction type(s).\n"));
  return;
  endif;
if(columns(output)!=1);#Exit function if value entered for output is of wrong length
  printf(strcat(output,": Invalid value for output type(s).\n"));
  return
  endif;
if(and(Fmat!="d",Fmat!="D",Fmat!="s",Fmat!="S"));
  printf(strcat(Fmat,": Invalid value for numerical format configuration.\n"));
  return
  endif;
if(output=="h");
  printf("Command-line argument format: picky w 0 1 d 1 1 0.95 all 44100 24 1 c:/sourcefile.tif\n");
  printf("All arguments are optional but must be entered in order.  Defaults are shown above.  They specify:\n");
  printf("w = type of output requested; enter 'picky o' for list of options\n");
  printf("0 = processing mode; enter 'picky m' for list of options\n");
  printf("1 = power to which source pixel intensity values will be raised before calculations are performed\n");
  printf("d = numerical format, lowercase d or s = double or single precision, simultaneous processing;\n"); 
  printf("    D or S = double and single precision, but with columns processed separately\n");
  printf("1 = adjustments of slope and DC offset; enter 'picky a' for list of options\n");
  printf("1 = sample-to-sample impulse threshold (0=0% of amplitude scale; 1=100% of amplitude scale)\n");
  printf("0.95 = decimal value to which amplitude values will be normalized for WAV output\n");
  printf("all = type(s) of transduction requested; enter 'picky t' for list of options\n");
  printf("44100 = sample rate for WAV output\n");
  printf("24 = bit depth for WAV output\n");
  printf("c:/sourcefile.tif = file path of source image, not available for batch processing modes (50+)\n");
  return;
  endif;
if(output=="a");
  printf("If present, adjustment type must be a number. Options are:\n");
  printf("2 = apply 20 Hz high pass filter to single image result\n");
  printf("3 = center endpoints of single image result on zero\n");
  printf("5 = apply 20 Hz high pass filter to concatenated batch result\n");
  printf("7 = center endpoints of concatenated batch result on zero\n");
  printf("11 = set first sample from image in batch equal to last sample from previous image\n");
  printf("To select multiple options, multiply them.\n");
  return;
  endif;
if(output=="o");
  printf("If present, output type must be a single character.  Options are:\n");
  printf("w = Wav, creates WAV file(s) for single images\n");
  printf("x = creates concatenated WAV file(s) for a batch, plus reference file with clicks at join points\n");
  printf("s = creates concatenated stereo WAV file(s) for a batch with clicks at join points in R channel\n");
  printf("v = Vector, saves vector(s) for a single image as .mat\n");
  printf("b = Batch, saves concatenated vector(s) for a batch as .mat\n"); 
  printf("m = Matrix, saves matrix/matrices for a batch as .mat\n");
  printf("h = Help, prints list of command-line arguments\n");
  printf("t = Transduction, prints list of transduction options\n");
  printf("o = Output, prints list of output options\n");
  printf("p = Processing mode, prints list of processing mode options\n");
  printf("a = Adjustments, prints list of adjustment options\n");
  return;
  endif;
if(output=="t");
  printf("If present, transduction type must be a three-character string. Options are:\n");
  printf("all = all four available types\n");
  printf("pdp = only position displacement (raw positional values)\n");
  printf("pvl = only position velocity (derivative/rate of change of positional values)\n");
  printf("idp = only intensity displacement (raw intensity values)\n");
  printf("ivl = only intensity velocity (derivative/rate of change of intensity values)\n");
  printf("pos = both position-based options (pdp and pvl)\n");
  printf("int = both intensity-based options (idp and ivl)\n");
  printf("dis = both displacement-based options (pdp and idp)\n");
  printf("vel = both velocity-based options (pvl and ivl)\n");
  printf("Any other three-character string will suppress file generation.\n");
  return;
  endif;
if(output=="p");
  printf("If present, processing mode type must be a number. Options are:\n");
  printf("0 = Process single image independently\n");
  printf("1 = Process single image independently; save amplitude range values\n");
  printf("2 = Process single image independently; overwrite any saved amplitude range values that are exceeded\n");
  printf("3 = Process single image using previously saved amplitude range values\n");
  printf("50 = Process multiple images independently from one another\n");
  printf("51 = Process multiple images to consistent amplitude range auto-detected for group\n");
  printf("52 = Process multiple images using previously saved amplitude range values\n");
  return;
  endif;
##Convert input variables to numbers if necessary
if isnumeric(Ipower)==0;
  Ipowerstr=Ipower;
  Ipower=str2num(Ipower);
  endif;
if isnumeric(Pc)==0;
  Pcstr=Pc;
  Pc=str2num(Pc);
  endif;
if isnumeric(wavscale)==0;
  wavscalestr=wavscale;
  wavscale=str2num(wavscale);
  endif;
if isnumeric(samplerate)==0;
  sampleratestr=samplerate;
  samplerate=str2num(samplerate);
  endif;
if isnumeric(bitdepth)==0;
  bitdepthstr=bitdepth;
  bitdepth=str2num(bitdepth);
  endif;
if isnumeric(Emode)==0;
  Emode=str2num(Emode);
  endif;
if isnumeric(corr1)==0;
  corr1str=corr1;
  corr1=str2num(corr1);
  endif;
##Check validity of numerical values
if(isempty(wavscale));
  printf(strcat(wavscalestr,": Invalid value for amplitude normalization.\n"));
  return;
  endif;
if(isempty(Ipower));
  printf(strcat(Ipowerstr,": Invalid value for power to which pixel intensity values will be raised.\n"));
  return;
  endif;
if(isempty(Pc));
  printf(strcat(Pcstr,": Invalid value for impulse rejection threshold.\n"));
  return;
  endif;
if(isempty(corr1));
  printf(strcat(corr1str,": Invalid value for slope and DC offset adjustment.\n"));
  return;
  endif;
if(or(isempty(samplerate),(samplerate<=0)));
  printf(strcat(sampleratestr,": Invalid value for sample rate.\n"));
  return;
  endif;
if(or(isempty(bitdepth),(samplerate<=0)));
  printf(strcat(bitdepthstr,": Invalid value for bit depth.\n"));
  return;
  endif;
##Convert corr to its prime factorization
corr=factor(corr1);
if (Emode<50);
  ##SINGLE-IMAGE PROCESSING ROUTINE
  if isnumeric(sourcefileinput);#Prompt user to select a file if none specified in command
    [A,B]=uigetfile({"*.gif;*.bmp;*.tif;*.jpg","Select an image file"},"Select an image file");
    if(A==0);
      printf("No source image selected.\n");
      return;
      endif;
    F=imread(strcat(B,A));
  else;
    [B,D,C]=fileparts(sourcefileinput);#Parse sourcefileinput into directory B and filename A and file extension C
    A=strcat(D,C);
    C=lower(C);#ensure file extension is lowercase
    ##Exit function if specified file extension is not supported
    if((strcmp(C,".gif").+strcmp(C,".bmp").+strcmp(C,".tif").+strcmp(C,".jpg"))==0);
      printf(strcat(A,": Specified file extension is not supported.\n"));
      return;
    else;
      F=imread(strcat(B,"/",A));#Read selected file as matrix F
      endif;
    endif;
  if(ndims(F)>3); #Exit function if source image has more than three dimensions
    printf("Source image has too many dimensions.\n");
    return;
    endif;
  if(ndims(F)==3); #Collapse 3-channel (color) source image to 2-channel 
    F=sum(F,3);
    printf(strcat("Multi-channel source image (",A,") summed to single channel\n"));
    endif;
  if(rows(F)>columns(F));#If source image is taller than it is wide, rotate ninety degrees counterclockwise
    F=rot90(F,1);
    endif;  
  if(or(Fmat=="D",Fmat=="S"));
    ##Column-by-column processing option for larger arrays
    for i=1:columns(F);
      Fcol=F(:,i);
      ##Convert to requested numerical format
      if(Fmat=="S");
        Fcol=single(Fcol);
      else;
        Fcol=double(Fcol);
        endif;
      Fcol=Fcol.^Ipower;#Raise source pixel intensity values to selected power
      Scol=sum(Fcol,1);#Create vector S of sums of pixel intensities by column
      S(i)=Scol;
      if(Scol==0);
        Qcol=0;#Sidestepping NaN values here to avoid division-by-zero warnings
      else;    
        Qcol=(Fcol./Scol);#Create vector Q of pixel intensities as fractions of total column intensity
        endif;
      if(i==1);
        Zcol(1:(rows(Fcol)))=1:rows(Fcol);#Create row vector Z of descending row numbers
        Zcol=rot90(Zcol,1);#Convert Z to column vector
        endif;
      Ycol=(Zcol.*Qcol);#Create vector Y of row numbers multiplied by pixel intensity
      Wcol=sum(Ycol);#Create vector W as sum of weighted row number values
      W(i)=Wcol;
      endfor;
  else;
    ##Convert to requested numerical format
    if(Fmat=="s");
      F=single(F);
    else;
      F=double(F);
      endif;
    F=F.^Ipower;#Raise source pixel intensity values to selected power
    S=sum(F,1);#Create vector S of sums of pixel intensities by column
    Q=(F./S);#Create matrix Q of pixel intensities as fractions of total column intensity
    Q(isnan(Q))=0;#replaces attempted divisions of zero by zero with 0
    Z(1:(rows(F)))=1:rows(F);#Create row vector Z of descending row numbers
    Z=rot90(Z,1);#Convert Z to column vector
    Y=(Z.*Q);#Create vector Y of row numbers multiplied by pixel intensity
    W=sum(Y,1);#Create vector W as sum of weighted row number values
    endif;
  ##Create derivative vectors
  T=diff(W);#Create vector T as derivative of W (converts displacement to velocity)
  N=diff(S);#Create vector N as derivative of S (converts displacement to velocity)
  ##Run impulse noise attenuation algorithm if requested
  if(and((Pc<1),(or(transduction=="all",transduction=="pdp",transduction=="pvl",transduction=="pos",transduction=="dis",transduction=="vel"))));
    if(abs(min(T))>max(T));#Sets maximum at initial level throughout processing
      Tmaxcl=abs(min(T));
    else;
      Tmaxcl=max(T);
      endif;
    if(Emode==3);
      load tdif.mat;
      if((Tdif/2)>Tmaxcl);
        Tmaxcl=(Tdif/2);
        endif;
      endif;
    for i=1:length(T);
      impset=1;
      if(abs(T(i))>Pc*Tmaxcl);#If a sample in vector T exceeds detection threshold
        if(i==1);
          T(i)=0;#if the first sample, set it equal to zero
        else;
          for j=1:10;#Searches up to 10 samples ahead for one that falls under threshold
            if((i+j)>length(T));
              break;
              endif;
            if(abs(T(i+j))<Pc*Tmaxcl);
              impset=2;
              impend=j;
              break;
              endif;
            endfor;
          if(impset==2);
            T(i)=((T(i-1)+T(i+impend))./2);#sets to average of nearest adjacent pair of accepted samples
          else;
            T(i)=T(i-1);#sets to value of preceding sample if no subsequent sample accepted
            endif;
          endif;
        endif;
      endfor;
    for i=2:(length(T)+1);#Re-integrate W
      W(i)=W(i-1)+T(i-1);
      endfor;
    endif;
  if(and((Pc<1),(or(transduction=="all",transduction=="idp",transduction=="ivl",transduction=="int",transduction=="dis",transduction=="vel"))));
    if(abs(min(N))>max(N));
      Nmaxcl=abs(min(N));
    else;
      Nmaxcl=max(N);
      endif;
    if(Emode==3);
      load ndif.mat;
      if((Ndif/2)>Nmaxcl);
        Nmaxcl=(Ndif/2);
        endif;
      endif;
    for i=1:length(N);
      impset=1;
      if(abs(N(i))>Pc*Nmaxcl);
        if(i==1);
          N(i)=0;
        else;
          for j=1:10;
            if((i+j)>length(N));
              break;
              endif;
            if(abs(N(i+j))<Pc*Nmaxcl);
              impset=2;
              impend=j;
              break;
              endif;
            endfor;
          if(impset==2);
            N(i)=((N(i-1)+N(i+impend))./2);
          else;
            N(i)=N(i-1);
            endif;
          endif;
        endif;
      endfor;
    for i=2:(length(N)+1);#Re-integrate S
      S(i)=S(i-1)+N(i-1);
      endfor;
    endif;
  ##Set median value to zero
  W=W-(median(W));
  S=S-(median(S));
  T=T-(median(T));
  N=N-(median(N));
  ##Set endpoints equal to zero if requested
  if(any(corr==3));
    Wstart=0-W(1);
    Wstop=(0-W(end))-Wstart;
    for i=1:length(W);
      W(i)=W(i)+Wstart+(Wstop.*(i/length(W)));
      endfor;
    Tstart=0-T(1);
    Tstop=(0-T(end))-Tstart;
    for i=1:length(T);
      T(i)=T(i)+Tstart+(Tstop.*(i/length(T)));
      endfor;
    Sstart=0-S(1);
    Sstop=(0-S(end))-Sstart;
    for i=1:length(S);
      S(i)=S(i)+Sstart+(Sstop.*(i/length(S)));
      endfor;
    Tstart=0-T(1);
    Tstop=(0-T(end))-Tstart;
    for i=1:length(T);
      T(i)=T(i)+Tstart+(Tstop.*(i/length(T)));
      endfor;
    endif;
  ##Apply Butterworth filter at single-image level if requested
  if(any(corr==2));
    Nyq=(samplerate/2);
    Order=2;#This will yield a 2nd-order Butterworth filter
    Hzcutoff=20;#Cut-off frequency set to 20 Hz.
    Cutoff=(Hzcutoff/Nyq);
    pkg load signal;
    [c,e]=butter(Order,Cutoff,"high");
    W=filter(c,e,W);
    S=filter(c,e,S);
    endif;
  ##Get maximum and minimum values
  Wmax=max(W);
  Wmin=min(W);
  Tmax=max(T);
  Tmin=min(T);
  Smax=max(S);
  Smin=min(S);
  Nmax=max(N);
  Nmin=min(N);
  if(Emode==1);#Save current range values
    Wdif=(Wmax-Wmin);
    save wdif.mat Wdif;
    Sdif=(Smax-Smin);
    save sdif.mat Sdif;
    Tdif=(Tmax-Tmin);
    save tdif.mat Tdif;
    Ndif=(Nmax-Nmin);
    save ndif.mat Ndif;
    endif;
  if(Emode==2);#If any current range values are greater than previously saved values, overwrite saved values
    load wdif.mat;
    load sdif.mat;
    load tdif.mat;
    load ndif.mat;
    if((Wmax-Wmin)>Wdif);
      Wdif=(Wmax-Wmin);
      save wdif.mat Wdif;
      endif;
    if((Smax-Smin)>Sdif);
      Sdif=(Smax-Smin);
      save sdif.mat Sdif;
      endif;
    if((Tmax-Tmin)>Tdif);
      Tdif=(Tmax-Tmin);
      save tdif.mat Tdif;
      endif;
    if((Nmax-Nmin)>Ndif);
      Ndif=(Nmax-Nmin);
      save ndif.mat Ndif;
      endif;
    endif;
  if(Emode==3);#Reset range to maximum with excess spread evenly between top and bottom
    load wdif.mat;
    load sdif.mat;
    load tdif.mat;
    load ndif.mat;
    if(Wdif>(Wmax-Wmin));
      Wmargin=(Wdif-(Wmax-Wmin));
      Wmax=(Wmax+(0.5.*Wmargin));
      Wmin=(Wmin-(0.5.*Wmargin));
      endif;
    if(Sdif>(Smax-Smin));
      Smargin=(Sdif-(Smax-Smin));
      Smax=(Smax+(0.5.*Smargin));
      Smin=(Smin-(0.5.*Smargin));
      endif;
    if(Tdif>(Tmax-Tmin));
      Tmargin=(Tdif-(Tmax-Tmin));
      Tmax=(Tmax+(0.5.*Tmargin));
      Tmin=(Tmin-(0.5.*Tmargin));
      endif;
    if(Ndif>(Nmax-Nmin));
      Nmargin=(Ndif-(Nmax-Nmin));
      Nmax=(Nmax+(0.5.*Nmargin));
      Nmin=(Nmin-(0.5.*Nmargin));
      endif;
    endif;
  ##Create scaled vectors from W, T, S, and N in -1 to +1 range (or 0 if no range at all)
  if(abs(Wmin)>Wmax);
    Wabsmax=abs(Wmin);
  else;
    Wabsmax=Wmax;
    endif;
  if(Wabsmax==0);
    Wscaled(1:length(W))=0;
  else;  
    Wscaled=(W./Wabsmax)*wavscale;
    endif;
  if(abs(Tmin)>Tmax);
    Tabsmax=abs(Tmin);
  else;
    Tabsmax=Tmax;
    endif;
  if(Tabsmax==0);
    Tscaled(1:length(T))=0;
  else;  
    Tscaled=(T./Tabsmax)*wavscale;
    endif;
  if(abs(Smin)>Smax);
    Sabsmax=abs(Smin);
  else;
    Sabsmax=Smax;
    endif;
  if(Sabsmax==0);
    Sscaled(1:length(S))=0;
  else;  
    Sscaled=(S./Sabsmax)*wavscale;
    endif;
  if(abs(Nmin)>Nmax);
    Nabsmax=abs(Nmin);
  else;
    Nabsmax=Nmax;
    endif;
  if(Nabsmax==0);
    Nscaled(1:length(N))=0;
  else;  
    Nscaled=(N./Nabsmax)*wavscale;
    endif;
  Timestring=strtok(num2str(time),".");#Get time value to add to WAV filenames
  ##Output selected scaled vectors as WAV or bumped vectors as .mat
  if or(transduction=="all",transduction=="pos",transduction=="dis",transduction=="pdp");
    if (output=="w");
      wavwrite(Wscaled,samplerate,bitdepth,strcat(strtok(A,"."),"_",Timestring,"_pdp.wav"));
      endif;
    if (output=="v");
      save w.mat W;
      endif;
    endif;
  if or(transduction=="all",transduction=="pos",transduction=="vel",transduction=="pvl");
    if (output=="w");
      wavwrite(Tscaled,samplerate,bitdepth,strcat(strtok(A,"."),"_",Timestring,"_pvl.wav"));
      endif;
    if (output=="v");
      save t.mat T;
      endif;
    endif;
  if or(transduction=="all",transduction=="int",transduction=="dis",transduction=="idp");
    if (output=="w");
      wavwrite(Sscaled,samplerate,bitdepth,strcat(strtok(A,"."),"_",Timestring,"_idp.wav"));
      endif;
    if (output=="v");
      save s.mat S;
      endif;
    endif;
  if or(transduction=="all",transduction=="int",transduction=="vel",transduction=="ivl");
    if (output=="w");
      wavwrite(Nscaled,samplerate,bitdepth,strcat(strtok(A,"."),"_",Timestring,"_ivl.wav"));
      endif;
    if (output=="v");
      save n.mat N;
      endif;
    endif;
  ##Notify of completion
  printf(strcat(A," processed mode_",num2str(Emode),"_",output,"_",transduction," ID=",Timestring,"\n"));
  if(range(W)==0);
    printf(strcat("Output value range for (",A,") is zero; recommend retry with different numerical format setting\n"));
    endif;
else;
  ##MULTIPLE-IMAGE PROCESSING ROUTINE
  ##Get multiple files
  [A,B]=uigetfile({"*.*","Select MULTIPLE image files"},"Select MULTIPLE image files.","MultiSelect","on");
  if(iscell(A)==0);
    printf(strcat("Requested processing mode (",num2str(Emode),") requires multiple image selection.\n")); 
    return;
    endif;
  ##Pad beginning of concatenated file with 3000 samples of silence if one will be created
  if(or((output=="x"),(output=="b"),(output=="s")));
    Wconc(1:3000)=0;
    Sconc(1:3000)=0;
    Tconc(1:3000)=0;
    Nconc(1:3000)=0;
    clickfile(1:3000)=0;
    endif;
  ##Preliminary batch processing loop for setting consistent amplitude levels, if requested
  if (Emode==51);
    picky("w",1,Ipower,Fmat,corr1,Pc,wavscale,"000",samplerate,bitdepth,strcat(B,A{1,1}));
    for i=2:columns(A);
      picky("w",2,Ipower,Fmat,corr1,Pc,wavscale,"000",samplerate,bitdepth,strcat(B,A{1,i}));
      endfor;
    printf("Amplitude levels set for batch\n");
    endif;
  ##Main batch processing loop
  for i=1:columns(A);
    ##Begin subroutine for batch result concatenation and matrix creation, if requested
    if(or((output=="x"),(output=="b"),(output=="s"),(output=="m")));
      ##Process each image using mode 0 for batch mode 50 or mode 3 for batch modes 51 and 52
      if (or((Emode==51),(Emode==52)));
        picky("v",3,Ipower,Fmat,corr1,Pc,wavscale,transduction,samplerate,bitdepth,strcat(B,A{1,i}));
      elseif (Emode==50);
        picky("v",0,Ipower,Fmat,corr1,Pc,wavscale,transduction,samplerate,bitdepth,strcat(B,A{1,i}));
        endif;
      ##Load the saved vectors
      load w.mat;
      load s.mat;
      load t.mat;
      load n.mat;
      ##Set first sample from new image equal to last sample from previous image if requested
      if(and(any(corr==11),i>1));
        Wchg=Wlast-W(1);
        Schg=Slast-S(1);
        Tchg=Tlast-T(1);
        Nchg=Nlast-N(1);
        W=W+Wchg;
        S=S+Schg;
        T=T+Tchg;
        N=N+Nchg;
        endif;
      ##Create matrix rows from loaded vectors
      Wmx(i,1:length(W))=W;
      Smx(i,1:length(S))=S;
      Tmx(i,1:length(T))=T;
      Nmx(i,1:length(N))=N;
      ##Append loaded vectors to concatenation vectors
      Wconc((end+1):(end+columns(W)))=W(1:columns(W));
      Sconc((end+1):(end+columns(W)))=S(1:columns(W));
      if (i==1);#Double first velocity sample of batch to keep in sync with clickfile
        Tconc(end+1)=T(1);
        Nconc(end+1)=N(1);
      else;#Calculate transitional velocity samples between image pairs
        Tconc(end+1)=Wlast-W(1);
        Nconc(end+1)=Slast-S(1);
        endif;
      Tconc((end+1):(end+columns(W)-1))=T(1:columns(W)-1);
      Nconc((end+1):(end+columns(W)-1))=N(1:columns(W)-1);
      ##Create reference vector 'clickfile' with clicks at image joins
      clickfile(end+1)=1;
      clickfile((end+1):(end+columns(W)-2))=0;
      clickfile(end+1)=(-1);
      ##Save last values of loaded vectors before next loop iteration
      Wlast=W(end);
      Slast=S(end);
      Tlast=T(end);
      Nlast=N(end);
      printf("Result concatenated\n");
    else; #If concatenation not chosen, then handle images separately
      if (Emode==50);
        picky(output,0,Ipower,Fmat,corr1,Pc,wavscale,transduction,samplerate,bitdepth,strcat(B,A{1,i}));
      elseif (or((Emode==51),(Emode==52))==1);
        picky(output,3,Ipower,Fmat,corr1,Pc,wavscale,transduction,samplerate,bitdepth,strcat(B,A{1,i}));
        endif;
      endif;
    endfor;
  ##Pad end of concatenation files with 3000 samples of silence if one is being created
  if(or((output=="x"),(output=="b"),(output=="s")));
    Wconc((end+1):(end+3000))=0;
    Sconc((end+1):(end+3000))=0;
    Tconc((end+1):(end+3000))=0;
    Nconc((end+1):(end+3000))=0;
    clickfile((end+1):(end+3000))=0;
    endif;
  ##Get time value to add to output filenames
  Timestring=strtok(num2str(time),".");
  ##Output any requested matrices for batch
  if (output=="m");
    if or(transduction=="all",transduction=="pos",transduction=="dis",transduction=="pdp");
      save wmatrix.mat Wmx;
      endif;
    if or(transduction=="all",transduction=="pos",transduction=="vel",transduction=="pvl");
      save tmatrix.mat Tmx;
      endif;
    if or(transduction=="all",transduction=="int",transduction=="dis",transduction=="idp");
      save smatrix.mat Smx;
      endif;
    if or(transduction=="all",transduction=="int",transduction=="vel",transduction=="ivl");
      save nmatrix.mat Nmx;
      endif;
    endif;
  ##Prepare and write other concatenated results if requested
  if(or((output=="x"),(output=="b"),(output=="s"))==1);
    ##Set endpoints equal to zero if requested
    if(any(corr==7));
      Wstart=0-Wconc(1);
      Wstop=(0-Wconc(end))-Wstart;
      for i=1:length(Wconc);
        Wconc(i)=Wconc(i)+Wstart+(Wstop.*(i/length(Wconc)));
        endfor;
      Tstart=0-Tconc(1);
      Tstop=(0-Tconc(end))-Tstart;
      for i=1:length(Tconc);
        Tconc(i)=Tconc(i)+Tstart+(Tstop.*(i/length(Tconc)));
        endfor;
      Sstart=0-Sconc(1);
      Sstop=(0-Sconc(end))-Sstart;
      for i=1:length(Sconc);
        Sconc(i)=Sconc(i)+Sstart+(Sstop.*(i/length(Sconc)));
        endfor;
      Nstart=0-Nconc(1);
      Nstop=(0-Nconc(end))-Nstart;
      for i=1:length(Nconc);
        Nconc(i)=Nconc(i)+Nstart+(Nstop.*(i/length(Nconc)));
        endfor;
      endif;
    ##Apply Butterworth filter at concatenated batch level if requested
    if(any(corr==5));
      Nyq=(samplerate/2);
      Order=2;#This will yield a 2nd-order Butterworth filter
      Hzcutoff=20;#Cut-off frequency set to 20 Hz.
      Cutoff=(Hzcutoff/Nyq);
      pkg load signal;
      [c,e]=butter(Order,Cutoff,"high");
      Wconc=Wconc-(0.5*(max(Wconc)));
      Sconc=Sconc-(0.5*(max(Sconc)));
      Tconc=Tconc-(0.5*(max(Tconc)));
      Nconc=Nconc-(0.5*(max(Nconc)));
      Wconc=filter(c,e,Wconc);
      Sconc=filter(c,e,Sconc);
      Tconc=filter(c,e,Tconc);
      Nconc=filter(c,e,Nconc);
      ##Reset initial and terminal 3000-sample silences to zero
      Wconc(1:3000)=0;
      Wconc((end-2999):end)=0;
      Sconc(1:3000)=0;
      Sconc((end-2999):end)=0;
      Tconc(1:3000)=0;
      Tconc((end-2999):end)=0;
      Nconc(1:3000)=0;
      Nconc((end-2999):end)=0;
      endif;
    ##Get maximum and minimum values
    Wmax=max(Wconc);
    Wmin=min(Wconc);
    Tmax=max(Tconc);
    Tmin=min(Tconc);
    Smax=max(Sconc);
    Smin=min(Sconc);
    Nmax=max(Nconc);
    Nmin=min(Nconc);
    ##Create scaled vectors from Wconc, Tconc, Sconc, and Nconc in -1 to +1 range (or 0 if no range at all)
    if(abs(Wmin)>Wmax);
      Wabsmax=abs(Wmin);
    else;
      Wabsmax=Wmax;
      endif;
    if(Wabsmax==0);
      Wcscaled(1:length(Wconc))=0;
    else;  
      Wcscaled=(Wconc./Wabsmax)*wavscale;
      endif;
    if(abs(Tmin)>Tmax);
      Tabsmax=abs(Tmin);
    else;
      Tabsmax=Tmax;
      endif;
    if(Tabsmax==0);
      Tcscaled(1:length(Tconc))=0;
    else;  
      Tcscaled=(Tconc./Tabsmax)*wavscale;
      endif;
    if(abs(Smin)>Smax);
      Sabsmax=abs(Smin);
    else;
      Sabsmax=Smax;
      endif;
    if(Sabsmax==0);
      Scscaled(1:length(Sconc))=0;
    else;  
      Scscaled=(Sconc./Sabsmax)*wavscale;
      endif;
    if(abs(Nmin)>Nmax);
      Nabsmax=abs(Nmin);
    else;
      Nabsmax=Nmax;
      endif;
    if(Nabsmax==0);
      Ncscaled(1:length(Nconc))=0;
    else;  
      Ncscaled=(Nconc./Nabsmax)*wavscale;
      endif;
    Timestring=strtok(num2str(time),".");#Get time value to add to WAV filenames
    ##Output concatenated batch WAV or vector if requested
    if(output=="x");
      wavwrite(clickfile,samplerate,bitdepth,strcat("Batch_",Timestring,"_clk.wav"));
      endif;
    if(output=="b");
      save clk.mat clickfile;
      endif;
    if or(transduction=="all",transduction=="pos",transduction=="dis",transduction=="pdp");
      if(output=="x");
        wavwrite(Wcscaled,samplerate,bitdepth,strcat("Batch_",Timestring,"_pdp.wav"));
        endif;
      if(output=="s");
        Wstereo=[clickfile;Wcscaled];
        Wstereo=rot90(Wstereo,-1);
        wavwrite(Wstereo,samplerate,bitdepth,strcat("Batch_",Timestring,"_pdpst.wav"));
        endif;
      if(output=="b");
        save wconc.mat Wconc;
        endif;
      endif;
    if or(transduction=="all",transduction=="pos",transduction=="vel",transduction=="pvl");
      if(output=="x");
        wavwrite(Tcscaled,samplerate,bitdepth,strcat("Batch_",Timestring,"_pvl.wav"));
        endif;
      if(output=="s");
        Tstereo=[clickfile;Tcscaled];
        Tstereo=rot90(Tstereo,-1);
        wavwrite(Tstereo,samplerate,bitdepth,strcat("Batch_",Timestring,"_pvlst.wav"));
        endif;
      if(output=="b");
        save tconc.mat Tconc;
        endif;
      endif;
    if or(transduction=="all",transduction=="int",transduction=="dis",transduction=="idp");
      if(output=="x");
        wavwrite(Scscaled,samplerate,bitdepth,strcat("Batch_",Timestring,"_idp.wav"));
        endif;
      if(output=="s");
        Sstereo=[clickfile;Scscaled];
        Sstereo=rot90(Sstereo,-1);
        wavwrite(Sstereo,samplerate,bitdepth,strcat("Batch_",Timestring,"_idpst.wav"));
        endif;
      if(output=="b");
        save sconc.mat Sconc;
        endif;
      endif;
    if or(transduction=="all",transduction=="int",transduction=="vel",transduction=="ivl");
      if(output=="x");
        wavwrite(Ncscaled,samplerate,bitdepth,strcat("Batch_",Timestring,"_ivl.wav"));
        endif;
      if(output=="s");
        Nstereo=[clickfile;Ncscaled];
        Nstereo=rot90(Nstereo,-1);
        wavwrite(Nstereo,samplerate,bitdepth,strcat("Batch_",Timestring,"_ivlst.wav"));
        endif;
      if(output=="b");
        save nconc.mat Nconc;
        endif;
      endif;
     endif;
  ##Notify of completion
  if(or(output=="x",output=="b",output=="m",output=="s")==1);
    printf(strcat("Concatenated batch results exported mode_",num2str(Emode),"_",output,"_",transduction," ID=",Timestring,"\n"));
    endif
  printf("Requested batch processing complete\n");
  endif;
endfunction;