function new_amp = func3d_shift_0527(PERCENT_MAX,PERCENT_MIN,STEP,file,sf,TST,TEND,n_formant);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %
% This is a program to creat a new spetrogram with %
% formant shifts and plot a 3-D figures showing    %
% the spectrogram of a wave file with one or more  %
% formants shift. The input is a 3-D spectrogram   %
% and the output is a 3-D spectrogram with F1, F2  %
% or/and F3 shift, which can be used to synthesize %
% formant-shifted vowels in STRAIGHT               %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% There are several parameters needing input or a default value in the
% program: 
% 1. Sampling frequency; 
% 2. Maximum and minimum percentage of formant shift; 
% 3. Number of steps for formant shift;
% 4. Starting and ending time (time range) for formant shift;

% Important vaiables
% F - the average formant frequencies over the whole time-frame 
% formant_freq - a 2-D variable for formant frequency at each time frame
% valley_freq - a 2-D variable for valley frequency at each time frame
% Updated by Chang Liu, 9.12.04

% Chang Liu, 7-28-1999
% Speech Psychophysics Lab, IU

% Add more comments for better understanding
% 12.11.2003

% Solve the problem for shifting a spectrum with flat valleys.
% Usually this problem is met when making a second shift following the
% first shift.
% 9.12.2004, Chang Liu

% debug the negative shift
% 2.5.05, Chang Liu

% Made function: DKP  2/10/05

% Eliminated shifted peaks from being cut off as they approach other peaks:
% ENO 5/27/05
tic;

% clear; Don't clear for coming synthesis

% Set F0, formant frequency and amplitude edges
F0MAX = 250;			  
F1MAX = 1000; %default is 1000
F2MAX = 2500;
F3MAX = 3500;
V4MAX = 4000;
AMPMAX = 25000;
fft_size = 1024;

%PERCENT_MAX = input('Input the MAX percetage for formant shift:');	% Maximum formant shift percentage
%PERCENT_MIN = input('Input the MIN percetage for formant shift:'); 	% Minimum formant shift percentage
%STEP = input('Input the number of steps for formant shift:');		% Steps of formant shift

for j = 1:STEP
   percent(j) = PERCENT_MIN + (PERCENT_MAX - PERCENT_MIN)/(STEP-1)*(j-1);
end   

%file = input('Input the file name: ', 's');
amp = load(file);

%sf = input('Input the sampling frequency:');
df = sf/fft_size;
% Straight analysis of amplitude as a funciton of frequency and time
[x y] = size(amp);
freq = (1:x)*df;
time = 1:y;
 
%TST = input ('Input the starting time of the vowel (ms):');				  % Starting time of formant shift
%TEND = input('Input the ending time of the vowel (ms):');                  % Ending time of formant shift

% Get peaks and valeys in the average spectrum over the 3-D spectrogram
meanamp = mean(amp(:, TST:TEND), 2);

for i = 1:(length(meanamp)-1)
  amp_diff(i) = meanamp(i+1) - meanamp(i);
end

number_peak = 0;
number_valey = 0;
for i = 1:(length(amp_diff) - 1)
   if ((amp_diff(i) > 0) & (amp_diff(i+1) < 0))
      number_peak = number_peak + 1;
      amp_peak(number_peak) = meanamp(i);
      index_peak(number_peak) = i;
      freq_peak(number_peak) = i*df;
   end
   
   if ((amp_diff(i) < 0) & (amp_diff(i+1) > 0))
      number_valey = number_valey + 1;
      amp_valey(number_valey) = meanamp(i);
      index_valey(number_valey) = i;
      freq_valey(number_valey) = i*df;
   end
end

% Get amplitude and frequency of F1, F2 and F3 in average
peak_f0 = 0;
peak_f1 = 0;
peak_f2 = 0;
peak_f3 = 0;
for i = 1:number_peak
   if (freq_peak(i) <=F0MAX)
      if (amp_peak(i) >= peak_f0)
         peak_f0 = amp_peak(i);
         F0 = freq_peak(i);
      end
   elseif ((freq_peak(i) > F0MAX) & (freq_peak(i) <= F1MAX))
      if (amp_peak(i) >= peak_f1)
         peak_f1 = amp_peak(i);
         F(1) = freq_peak(i);
      end
   elseif ((freq_peak(i) > F1MAX) & (freq_peak(i) <= F2MAX))
      if (amp_peak(i) >= peak_f2)
         peak_f2 = amp_peak(i);
         F(2) = freq_peak(i);
      end
   elseif ((freq_peak(i) > F2MAX) & (freq_peak(i) < F3MAX))
      if (amp_peak(i) >= peak_f3)
         peak_f3 = amp_peak(i);
         F(3) = freq_peak(i);
      end
   end
end

%F      %Remove the '%' to cause the average formant frequency used to
         %calculate shifting endpoints to be displayed in the main matlab window.
         

% Get amplitude and frequency of F1, F2 and F3 for each time frame
formant_index = zeros(3, y);        % index number of the first three formants
formant_freq = zeros(3,y);          % first three formant freuqencies initiated at 0 Hz
formant_amp = zeros(3, y);          % amplitude of the first three formants
valley_index = zeros(4, y);         % index nubmers of the first four valley 
valley_freq = zeros(4,y);           % first four valley frequeices initiated at 0 Hz
valley_amp = zeros(4, y);           % amplitude of the first four valleys

for i = TST:TEND
    formant_index(1, i) = max(find(amp(:,i) == max(amp(find(freq>=0&freq<F1MAX), i))));
    formant_freq (1, i) = freq(formant_index(1, i));
    formant_amp (1, i) = amp(formant_index(1, i), i);
    formant_index(2, i) = max(find(amp(:,i) == max(amp(find(freq>=F1MAX&freq<F2MAX), i))));
    formant_freq (2, i) = freq(formant_index(2, i));
    formant_amp(2, i) = amp(formant_index(2, i), i);
    formant_index(3, i) = max(find(amp(:,i) == max(amp(find(freq>=F2MAX&freq<F3MAX), i))));
    formant_freq (3, i) = freq(formant_index(3, i));
    formant_amp(3, i) = amp(formant_index(3, i), i);
    
    valley_index (1, i) = min(find(amp(:,i) == min(amp(find(freq>=0&freq<formant_freq(1,i)), i))));
    valley_freq (1, i) = freq(valley_index(1, i));
    valley_amp (1, i) = amp(valley_index(1, i));
    valley_index (2, i) = min(find(amp(:,i) == min(amp(find(freq>=formant_freq(1,i)&freq<formant_freq(2,i)), i))));
    valley_freq (2, i) = freq(valley_index(2, i));
    valley_amp (2, i) = amp(valley_index(2, i));
    valley_index (3, i) = min(find(amp(:,i) == min(amp(find(freq>=formant_freq(2,i)&freq<formant_freq(3,i)), i))));
    valley_freq (3, i) = freq(valley_index(3, i));
    valley_amp (3, i) = amp(valley_index(3, i));
    valley_index (4, i) = min(find(amp(:,i) == min(amp(find(freq>=formant_freq(3,i)&freq<V4MAX), i))));
    valley_freq (4, i) = freq(valley_index(4, i));
    valley_amp (4, i) = amp(valley_index(4, i));
end
    
% Get shift for one formant
   for j = 1:STEP
      for i = 1:length(freq)
         new_freq(i,j) = freq(i);
         new_amp(i,:,j) = amp(i,:);
      end
   end
   
  % Get the frequence shift values (Hz)for all the percent shifts
 % n_formant = input('Enter which formant you want to shift (1, 2 or 3): ');
  shift_freq = F(n_formant)*percent;
  index_change = floor((F(n_formant) + F(n_formant)*percent)/df) - floor(F(n_formant)/df);
  
 
  % calculate new frequency and new amplitude; 
  % Fill in points on low-frequency valley, no change high-frequency valley
  for j = 1:STEP
      for k = TST:TEND            % time range for formant shift      
          if (shift_freq(j) >= 0) % formant shifts to the right (positively)
              for i = (valley_index(n_formant, k) + index_change(j)) : (valley_index(n_formant+1, k) - 2)   %CHANGE ENO
                  if (i-index_change(j)) <= 1
                      new_amp(i,k,j) = amp(i-index_change(j), k);
                  else
                      x = [freq(i - index_change(j) - 1) freq(i - index_change(j))]';
                      y = [amp(i-index_change(j)-1, k) amp(i-index_change(j), k)]';
   	                  new_amp(i,k,j) = interp1q(x, y, freq(i)-shift_freq(j));
                  end
        	  end

              for i = (valley_index(n_formant+1, k) - 1):(valley_index(n_formant+1, k) - 1 + index_change(j)) %CHANGE ENO
                      x = [freq(i - index_change(j) - 1) freq(i - index_change(j))]';
                      y = [amp(i-index_change(j)-1, k) amp(i-index_change(j), k)]';
   	                  
                     if new_amp(i,k,j) > interp1q(x, y, freq(i)-shift_freq(j))  %CHANGE ENO
                          %nothing happens
                     else
                         new_amp(i,k,j) = interp1q(x, y, freq(i)-shift_freq(j)); %CHANGE ENO
                     end           %CHANGE ENO
              end
              
              for i = valley_index(n_formant, k) : (valley_index(n_formant, k) + index_change(j) - 1)
                  new_amp(i,k,j) = amp(valley_index(n_formant, k), k);
              end
          else 	% formant shifts to the left
               for i = (valley_index(n_formant, k) + 1):(valley_index(n_formant+1, k)+index_change(j)) %CHANGE ENO
                      x = [freq(i - index_change(j) - 1) freq(i - index_change(j))]';
                      y = [amp(i-index_change(j)-1, k) amp(i-index_change(j), k)]';
                      new_amp(i,k,j) = interp1q(x, y, freq(i)-shift_freq(j));
               end
              
              for i = (valley_index(n_formant, k) + index_change(j)):valley_index(n_formant, k) %CHANGE ENO
                      x = [freq(i - index_change(j) - 1) freq(i - index_change(j))]';
                      y = [amp(i-index_change(j)-1, k) amp(i-index_change(j), k)]';
   	                  
                     if new_amp(i,k,j) > interp1q(x, y, freq(i)-shift_freq(j))  %CHANGE ENO
                          %nothing happens
                     else
                         new_amp(i,k,j) = interp1q(x, y, freq(i)-shift_freq(j)); 
                     end           
              end
            
	          for i = (valley_index(n_formant + 1, k) + index_change(j) + 1) : (valley_index(n_formant + 1, k)) 
   		          new_amp(i, k, j) = amp(valley_index(n_formant + 1, k), k);
	          end
          end
      end
  end
 
    
  % Get a 3-D graph of amplitude as a function of time and frequency. 
  
  % First, set up the frequency range for the 3-D graph, depending on which
  % formant is shifted. 
  
  if (n_formant == 1)
      f_start = 200;   
      f_end = 1200;	 	
  elseif (n_formant == 2)
      f_start = 1000;   
      f_end = 3000;	 	
  elseif (n_formant == 3)
      f_start = 2000;   
      f_end = 4000;	 	
  else
      f_start = 200;   
      f_end = 4000;	 	
  end    


  % starting and ending frequency
  f_startPoint = floor(f_start/df);
  f_endPoint = floor(f_end/df);
  az = 75;			% azimuth angle to frequency axis
  el = 30;			% elevation angle 

  
%for j = 1:STEP      % plot all the graphs with each formant-shift step
% Just plot STEPth figure
figure;
j=STEP;
  subplot(1,2,1);       % standard vowel
  mesh(time, freq(f_startPoint:f_endPoint), amp(f_startPoint:f_endPoint,:));
  view([az el]);
  axis([0 length(time) f_start f_end  0 AMPMAX]);
  xlabel('Time (ms)');
  ylabel('Frequency (Hz)');
  zlabel('Amplitdue (dB)');
  title('Original Spectrogram');

  subplot(1,2,2);       % vowel with formant shift
  mesh(time, new_freq(f_startPoint:f_endPoint,j), new_amp(f_startPoint:f_endPoint,:,j));
  view([az el]);
  axis([0 length(time) f_start f_end  0 AMPMAX]);
  xlabel('Time (ms)');
  ylabel('Frequency (Hz)');
  zlabel('Amplitdue (dB)');
  title('Spectrogram with formant shift');
  
  %end

toc