tomwalters@0: %
tomwalters@0: %  AsymCmpFrsp Version 2
tomwalters@0: %  Amplitude Spectrum of Asymmetric Compensation IIR filter for the gammachirp 
tomwalters@0: %  corresponding to MakeAsymCmpFiltersV2.m
tomwalters@0: %
tomwalters@0: %  Toshio Irino
tomwalters@0: %  Original : 14 Apr. 99 (Version1)
tomwalters@0: %  Modified : 11 Jun 2004
tomwalters@0: %  Modified :  7 Jul 2005 % NfrqRsl
tomwalters@0: %
tomwalters@0: %  function [ACFFrsp, freq, AsymFunc]
tomwalters@0: %            = AsymCmpFrspV2(Frs,fs,b,c,NfrqRsl,NumFilt)
tomwalters@0: %
tomwalters@0: %  INPUT:    fs: Sampling frequency
tomwalters@0: %            Frs: array of the center frequencies
tomwalters@0: %            b : array or scalar of a bandwidth coefficient 
tomwalters@0: %            c : array or scalar of asymmetric parameters 
tomwalters@0: %            NfrqRsl: freq. resolution
tomwalters@0: %            NumFilt: Number of 2nd-order filters default 4
tomwalters@0: %  OUTPUT:   ACFFrsp: abs(Frsp of ACF)  (NumCh*NfrqRsl matrix)
tomwalters@0: %            freq: freq.                    (1*NfrqRsl vector)
tomwalters@0: %            AsymFunc:    Original Asymmetric Function (NumCh*NfrqRsl matrix)
tomwalters@0: %
tomwalters@0: %
tomwalters@0: function [ACFFrsp,freq,AsymFunc] = AsymCmpFrspV2(Frs,fs,b,c,NfrqRsl,NumFilt)
tomwalters@0: 
tomwalters@0: if nargin < 1, help AsymCmpFrspV2, end;
tomwalters@0: if nargin < 5, NfrqRsl = []; end;
tomwalters@0: if length(NfrqRsl) == 0,  NfrqRsl = 1024; end;
tomwalters@0: if nargin < 6, NumFilt = []; end;
tomwalters@0: if length(NumFilt) == 0,  NumFilt = 4; end; % default
tomwalters@0: if NumFilt ~= 4, error('NumFilter should be 4.'); end;
tomwalters@0: 
tomwalters@0: Frs = Frs(:);
tomwalters@0: b = b(:);
tomwalters@0: c = c(:);
tomwalters@0: NumCh = length(Frs);
tomwalters@0: 
tomwalters@0: SwCoef = 0;  % self consitency
tomwalters@0: %SwCoef = 1; % referece to MakeAsymCmpFiltersV2
tomwalters@0: 
tomwalters@0: if SwCoef == 0
tomwalters@0: % New Coefficients. NumFilter = 4; See [1]
tomwalters@0:   p0 = 2;
tomwalters@0:   p1 = 1.7818 .* (1 - 0.0791*b) .* (1 - 0.1655*abs(c));
tomwalters@0:   p2 = 0.5689 .* (1 - 0.1620*b) .* (1 - 0.0857*abs(c));
tomwalters@0:   p3 = 0.2523 .* (1 - 0.0244*b) .* (1 + 0.0574*abs(c));
tomwalters@0:   p4 = 1.0724;
tomwalters@0: else
tomwalters@0:   
tomwalters@0:   ACFcoef = MakeAsymCmpFiltersV2(fs,Frs,b,c);
tomwalters@0: end;
tomwalters@0:   
tomwalters@0: [dummy ERBw] = Freq2ERB(Frs);
tomwalters@0: freq = (0:NfrqRsl-1)/NfrqRsl*fs/2;
tomwalters@0: ACFFrsp = ones(NumCh,NfrqRsl);
tomwalters@0: freq2 = [ones(NumCh,1)*freq, Frs];
tomwalters@0: 
tomwalters@0: 
tomwalters@0: for Nfilt = 1:NumFilt
tomwalters@0: 
tomwalters@0:   if SwCoef == 0,
tomwalters@0:     r  = exp(-p1.*(p0./p4).^(Nfilt-1) .* 2 .* pi .*b .*ERBw /fs); 
tomwalters@0:     delfr = (p0*p4)^(Nfilt-1).*p2.*c.*b.*ERBw; 
tomwalters@0:     phi = 2*pi*max(Frs + delfr,0)/fs;
tomwalters@0:     psy = 2*pi*max(Frs - delfr,0)/fs;
tomwalters@0:     fn = Frs;    
tomwalters@0:     ap = [ones(NumCh,1), -2*r.*cos(phi),  r.^2];
tomwalters@0:     bz = [ones(NumCh,1), -2*r.*cos(psy),  r.^2]; 
tomwalters@0:   else
tomwalters@0:     ap = ACFcoef.ap(:,:,Nfilt);
tomwalters@0:     bz = ACFcoef.bz(:,:,Nfilt);
tomwalters@0:   end;
tomwalters@0: 
tomwalters@0:   cs1 = cos(2*pi*freq2/fs);
tomwalters@0:   cs2 = cos(4*pi*freq2/fs);
tomwalters@0: 
tomwalters@0:   bzz0 = (bz(:,1).^2 + bz(:,2).^2 + bz(:,3).^2 )*ones(1,NfrqRsl+1);
tomwalters@0:   bzz1 = (2* bz(:,2).*( bz(:,1) + bz(:,3)))*ones(1,NfrqRsl+1);
tomwalters@0:   bzz2 = (2* bz(:,1).*bz(:,3))*ones(1,NfrqRsl+1);
tomwalters@0:   hb = bzz0 + bzz1.* cs1 + bzz2.* cs2;
tomwalters@0: 
tomwalters@0:   app0 = (ap(:,1).^2 + ap(:,2).^2 + ap(:,3).^2)*ones(1,NfrqRsl+1);
tomwalters@0:   app1 = (2* ap(:,2).*( ap(:,1) + ap(:,3)))*ones(1,NfrqRsl+1);
tomwalters@0:   app2 = (2* ap(:,1).*ap(:,3))*ones(1,NfrqRsl+1);
tomwalters@0:   ha = app0 + app1.*cs1 + app2.*cs2;
tomwalters@0: 
tomwalters@0:   H = sqrt(hb./ha);
tomwalters@0:   Hnorm = H(:,NfrqRsl+1)*ones(1,NfrqRsl);    % Noimalization by fn value
tomwalters@0: 
tomwalters@0:   ACFFrsp = ACFFrsp .* H(:,1:NfrqRsl)./ Hnorm;
tomwalters@0: 
tomwalters@0: end;
tomwalters@0: 
tomwalters@0: %%% original Asymmetric Function without shift centering 
tomwalters@0: fd = (ones(NumCh,1)*freq - Frs*ones(1,NfrqRsl));
tomwalters@0: be = (b.*ERBw)*ones(1,NfrqRsl);
tomwalters@0: cc = (c.*ones(NumCh,1))*ones(1,NfrqRsl); % in case when c is scalar
tomwalters@0: AsymFunc = exp(cc.*atan2(fd,be));
tomwalters@0: 
tomwalters@0: return