tomwalters@0: %	Fp2toFr1
tomwalters@0: %	   derives fr1 from the given fp2
tomwalters@0: %
tomwalters@0: %	function [fr1,fp1] = Fp2toFr1(n,b1,c1,b2,c2,frat,fp2)
tomwalters@0: %
tomwalters@0: %	Author:  Masashi Unoki
tomwalters@0: %	Created:  3 July 2002
tomwalters@0: %	Updated: 15 July 2002
tomwalters@0: %	Revised:  9 Oct. 2003 (checked & renamed variables)
tomwalters@0: %	Copyright (c) 2002, AIS-Lab. JAIST
tomwalters@0: %
tomwalters@0: function [fr1,fp1] = Fp2toFr1(n,b1,c1,b2,c2,frat,fp2)
tomwalters@0: if nargin < 1; help Fp2toFr1; end;
tomwalters@0: 
tomwalters@0: SR=24000;
tomwalters@0: Nfft=1024*2;
tomwalters@0: 
tomwalters@0: %%%%%%% Coefficients: ERB(fr1)=alp1*fr1+alp0 %%%%%%%
tomwalters@0: alp1=24.7*4.37/1000;
tomwalters@0: alp0=24.7;
tomwalters@0: 
tomwalters@0: %%%%%%% Coefficients: fr2=bet1*fr2+bet0 %%%%%%%
tomwalters@0: bet1=frat*(1+c1*b1*alp1/n);
tomwalters@0: bet0=frat*c1*b1*alp0/n;
tomwalters@0: 
tomwalters@0: %%%%%%% Coefficients: ERB(fr2)=zet1*fr1+zet0 %%%%%%%
tomwalters@0: zet1=alp1*bet1;
tomwalters@0: zet0=alp1*bet0+alp0;
tomwalters@0: 
tomwalters@0: %%%%%%%	D1*fr1^3 + D2*fr1^2 + D3*fr1 + D4 = 0 %%%%%%%
tomwalters@0: 
tomwalters@0: D1=((b2^2*zet1^2+bet1^2)*(c1*b1*alp1+n) + (c2*b2*zet1)*(b1^2*alp1^2+1));
tomwalters@0: D2=((b2^2*zet1^2+bet1^2)*(c1*b1*alp0-n*fp2) ...
tomwalters@0:     + (2*b2^2*zet1*zet0-2*bet1*(fp2-bet0))*(c1*b1*alp1+n) ...
tomwalters@0:     + (c2*b2*zet1)*(2*b1^2*alp1*alp0-2*fp2) + (c2*b2*zet0)*(b1^2*alp1^2+1));
tomwalters@0: D3=((2*b2^2*zet1*zet0-2*bet1*(fp2-bet0))*(c1*b1*alp0-n*fp2) ...
tomwalters@0:     + (b2^2*zet0^2+(fp2-bet0)^2)*(c1*b1*alp1+n)...
tomwalters@0:     +(c2*b2*zet1)*(b1^2*alp0^2+fp2^2) + (c2*b2*zet0)*(2*b1^2*alp1*alp0-2*fp2) );
tomwalters@0: D4=(b2^2*zet0^2+(fp2-bet0)^2)*(c1*b1*alp0-n*fp2) ...
tomwalters@0:     + (c2*b2*zet0)*(b1^2*alp0^2+fp2^2);
tomwalters@0: 
tomwalters@0: q=roots([D1 D2 D3 D4]);
tomwalters@0: candFr1=q(imag(q)==0);
tomwalters@0: 
tomwalters@0: LenFr1=length(candFr1);
tomwalters@0: if (LenFr1 > 1)		% finding the maximum peak of |Gcc(f)|
tomwalters@0:    GccAtFp2=zeros(1,LenFr1);
tomwalters@0:    for m=1:LenFr1
tomwalters@0:       fr1m=candFr1(m);
tomwalters@0:       fr2m=bet1*fr1m+bet0;
tomwalters@0:       [GcFrsp1, freq1]=GammaChirpFrsp(fr1m,SR,n,b1,c1,0,Nfft);
tomwalters@0:       [dummy ERBw2] = Freq2ERB(fr2m);
tomwalters@0:       AsymFuncFrsp = exp(c2*atan2((freq1 - fr2m),(b2*ERBw2)));
tomwalters@0:       GcFrsp = GcFrsp1.*AsymFuncFrsp;
tomwalters@0:       GcFrsp = GcFrsp/max(GcFrsp);
tomwalters@0:       [dummy,pos]=min(abs(freq1-fp2));
tomwalters@0:       GccAtFp2(m)=GcFrsp(pos);
tomwalters@0:    end
tomwalters@0:    [dummy,pos]=max(GccAtFp2);
tomwalters@0:    fr1=candFr1(pos);
tomwalters@0: else
tomwalters@0:    fr1=candFr1;
tomwalters@0: end
tomwalters@0: 
tomwalters@0: fp1=fr1+c1*b1*(alp1*fr1+alp0)/n;
tomwalters@0: 
tomwalters@0: return