% TMDF_EXMAPLE (S) Calculates the second harmonic of kantele using the
% downsampled fundamental according to the TMDF formulation.

% CE, 05.2002, for the JASA Kantele Paper

% NEEDS: SIGNAL PROCESSING TOOLBOX 
% PRODUCES: FIGURE 10

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% PARAMETERS
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 E = 2*10^9;
 r = (0.35/2)*10^(-3);
 S = pi * r^2;
 L = 0.4015;
 
 TMDF = pi^2*E*S/8/L^2;
 
 a  = sqrt(2)*2*10^(-3);
 as = 2 * a^2;
 
 fs = 44100;
 idx_beg = 1429; % PRECALCULATED ONSET
 f0 = 392.7265;
 
 % LONGITUDINAL ADMITTANCE MAGNITUDE AT 2f0, PRECALCULATED
 LON_ADD = 9.2450;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Read the measurement
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 wavedata = wavread('tmdf_example.wav');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Index of the onset =  1429
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 wavedata = wavedata(idx_beg:length(wavedata));
 len_sig = length(wavedata);
 [val_a,idx_a]=max(wavedata);

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %% Get the hopsize, make an even envelope size
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 period_f0 = fs/f0;
 hop_size = fix(period_f0 * 2.5 / 4);
 env_size = fix(len_sig/hop_size)-...
            (~iseven(fix(len_sig/hop_size)));
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % TUNE BUTTERWORTHs FOR 1st and 2nd harmonics
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 W1 = [0.95*f0 1.05*f0]/22050;
 W2 = [1.95*f0 2.05*f0]/22050;
  
 [B1,A1] = BUTTER(2,W1);
 [B2,A2] = BUTTER(2,W2);
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
 % ZERO-PHASE FILTERING
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 harm1 = filtfilt(B1,A1,wavedata);
 harm2 = filtfilt(B2,A2,wavedata);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % ENVELOPE EXTRACTION
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 harm1 = harm1(1:hop_size*env_size);
 harm2 = harm2(1:hop_size*env_size);

 ENV1 = reshape(harm1,hop_size,env_size);
 ENV2 = reshape(harm2,hop_size,env_size);
 
 ENV1 = max(abs(ENV1));
 ENV2 = max(abs(ENV2));
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % DOWNSAMPLE THE FIRST HARMONIC
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 harm1_d = resample(harm1,1,2);
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % SCALE THE DS FUNDAMENTAL TO UNITY AMPLITUDE
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 HMAX1 = max(abs(harm1));
 HMAX2 = max(abs(harm2));
 HMAX3 = max(abs(harm1_d));

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % CALCULATE THE GENERATED SECOND HARMONIC
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 ad_harm2 = TMDF * as * harm1_d/HMAX3;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % CALCULATE THE ENVELOPE OF GENERATED SECOND HARMONIC
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 ENV2_D = reshape(ad_harm2,hop_size,env_size/2);
 ENV2_D = max(abs(ENV2_D));
 
 ENV1 = reshape(harm1,hop_size,env_size);
 ENV2 = reshape(harm2,hop_size,env_size);
 
 ENV1 = max(abs(ENV1));
 ENV2 = max(abs(ENV2));
 
 idx_0 = (0:env_size-1); idx_0 = idx_0 * hop_size/44100;
 idx_2 = (0:env_size/2-1); 
 idx_2 = idx_2 * hop_size/44100;
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % ADMITTANCE NORMALIZATION; 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 AMP_NOR = -max(20*log10(ENV1)); %dB
 PHA_NOR = 0.03; %s
    
 ENV1 = 20*log10(ENV1)+AMP_NOR;
 ENV2 = 20*log10(ENV2)+AMP_NOR;
 ENV2_D = 20*log10(ENV2_D)+AMP_NOR+LON_ADD; % ADMITTANCE COMPENSATION

 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % FIGURE
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 figure(1);clf;
  plot(idx_0,ENV1,'m-.',...
       idx_0,ENV2,'r--',...
       idx_2+PHA_NOR,ENV2_D,'b'); 

  set(gca,'FontName','TimesNewRoman');
  set(gca,'FontSize',16);
 
  xlabel('Time (s)');
  ylabel('Normalized Velocity (dB)');
  grid on,axis('tight')
  set(gca,'XLim',[0 1.5])
  set(gca,'YLim',[-40 0])