-
Notifications
You must be signed in to change notification settings - Fork 1
/
AsymCmpFrspV2.m
116 lines (98 loc) · 3.46 KB
/
AsymCmpFrspV2.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
%
% AsymCmpFrsp Version 2
% Amplitude Spectrum of Asymmetric Compensation IIR filter for the gammachirp
% corresponding to MakeAsymCmpFiltersV2.m
%
%
% Toshio Irino
% Original : 14 Apr. 99 (Version1)
% Modified : 11 Jun 2004
% Modified : 7 Jul 2005 % NfrqRsl
% Modified : 26 Apr 2015 % Response at Frs (not linear-freq bin)
%
% function [ACFFrsp, freq, AsymFunc]
% = AsymCmpFrspV2(Frs,fs,b,c,NfrqRsl,NumFilt)
%
% INPUT: fs: Sampling frequency
% Frs: array of the center frequencies
% b : array or scalar of a bandwidth coefficient
% c : array or scalar of asymmetric parameters
% NfrqRsl: >64) freq. resolution for linear frequency scale
% 0) for specify resp at Frs
% NumFilt: Number of 2nd-order filters default 4
% OUTPUT: ACFFrsp: abs(Frsp of ACF) (NumCh*NfrqRsl matrix)
% freq: freq. (1*NfrqRsl vector)
% AsymFunc: Original Asymmetric Function (NumCh*NfrqRsl matrix)
%
%
function [ACFFrsp,freq,AsymFunc] = AsymCmpFrspV2(Frs,fs,b,c,NfrqRsl,NumFilt)
if nargin < 1, help(mfilename), end;
if nargin < 5, NfrqRsl = []; end;
if length(NfrqRsl) == 0, NfrqRsl = 1024; end;
if nargin < 6, NumFilt = []; end;
if length(NumFilt) == 0, NumFilt = 4; end; % default
if NumFilt ~= 4, error('NumFilter should be 4.'); end;
Frs = Frs(:);
b = b(:);
c = c(:);
NumCh = length(Frs);
if NfrqRsl >= 64,
freq = (0:NfrqRsl-1)/NfrqRsl*fs/2;
elseif NfrqRsl == 0, % added 26 Apr 2015
freq = Frs(:)';
NfrqRsl = length(freq);
else
help(mfilename);
error('Specify NfrqRsl 0) for Frs or N>=64 for linear-freq scale');
end;
%% coef
SwCoef = 0; % self consitency
%SwCoef = 1; % referece to MakeAsymCmpFiltersV2
if SwCoef == 0
% New Coefficients. NumFilter = 4; See [1]
p0 = 2;
p1 = 1.7818 .* (1 - 0.0791*b) .* (1 - 0.1655*abs(c));
p2 = 0.5689 .* (1 - 0.1620*b) .* (1 - 0.0857*abs(c));
p3 = 0.2523 .* (1 - 0.0244*b) .* (1 + 0.0574*abs(c));
p4 = 1.0724;
else
ACFcoef = MakeAsymCmpFiltersV2(fs,Frs,b,c);
end;
%% filter coef.
% No change at 26 Apr 2015
[dummy ERBw] = Freq2ERB(Frs);
ACFFrsp = ones(NumCh,NfrqRsl);
freq2 = [ones(NumCh,1)*freq, Frs];
for Nfilt = 1:NumFilt
if SwCoef == 0,
r = exp(-p1.*(p0./p4).^(Nfilt-1) .* 2 .* pi .*b .*ERBw /fs);
delfr = (p0*p4)^(Nfilt-1).*p2.*c.*b.*ERBw;
phi = 2*pi*max(Frs + delfr,0)/fs;
psy = 2*pi*max(Frs - delfr,0)/fs;
fn = Frs;
ap = [ones(NumCh,1), -2*r.*cos(phi), r.^2];
bz = [ones(NumCh,1), -2*r.*cos(psy), r.^2];
else
ap = ACFcoef.ap(:,:,Nfilt);
bz = ACFcoef.bz(:,:,Nfilt);
end;
cs1 = cos(2*pi*freq2/fs);
cs2 = cos(4*pi*freq2/fs);
bzz0 = (bz(:,1).^2 + bz(:,2).^2 + bz(:,3).^2 )*ones(1,NfrqRsl+1);
bzz1 = (2* bz(:,2).*( bz(:,1) + bz(:,3)))*ones(1,NfrqRsl+1);
bzz2 = (2* bz(:,1).*bz(:,3))*ones(1,NfrqRsl+1);
hb = bzz0 + bzz1.* cs1 + bzz2.* cs2;
app0 = (ap(:,1).^2 + ap(:,2).^2 + ap(:,3).^2)*ones(1,NfrqRsl+1);
app1 = (2* ap(:,2).*( ap(:,1) + ap(:,3)))*ones(1,NfrqRsl+1);
app2 = (2* ap(:,1).*ap(:,3))*ones(1,NfrqRsl+1);
ha = app0 + app1.*cs1 + app2.*cs2;
H = sqrt(hb./ha);
Hnorm = H(:,NfrqRsl+1)*ones(1,NfrqRsl); % Noimalization by fn value
ACFFrsp = ACFFrsp .* H(:,1:NfrqRsl)./ Hnorm;
end;
%%% original Asymmetric Function without shift centering
fd = (ones(NumCh,1)*freq - Frs*ones(1,NfrqRsl));
be = (b.*ERBw)*ones(1,NfrqRsl);
cc = (c.*ones(NumCh,1))*ones(1,NfrqRsl); % in case when c is scalar
AsymFunc = exp(cc.*atan2(fd,be));
return