-
Notifications
You must be signed in to change notification settings - Fork 13
/
Copy pathtst_cifa.m
100 lines (95 loc) · 1.86 KB
/
tst_cifa.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
% test gammatone filterbank analysis
function tst_cifa
td=4;
load('test/tst_cifa')
y=real(y);
t=1000*(0:(length(y)-1))/rate;
figure(1);clf
reset_color(1);
plot(t,y)
title('complex IIR filterbank impulse response')
y_lim=[min(min(y)) max(max(y))]* 1.05;
t_lim=[-0.5 15.5];
axis([t_lim y_lim])
xlabel('time (ms)')
figure(2);clf
H=ffa(y);
f=linspace(0,rate/2000,length(H))';
d=td*ones(size(f));
fm=max(max(f));
m_lim=[0.05 fm -50 10];
d_lim=[0.05 fm 0 8];
M=db(H);
D=gd(H,f);
D(M<-40)=NaN;
subplot(2,1,1)
reset_color(1);
semilogx(f,M)
axis(m_lim)
ylabel('magnitude (dB)')
title('complex IIR filterbank transfer functions')
subplot(2,1,2)
reset_color(1);
semilogx(f,D,f,d,':k')
axis(d_lim)
xlabel('frequency (kHz)')
ylabel('delay (ms)')
% plot filter-sum transfer function
figure(3);clf
H=ffa(sum(y,2));
f=linspace(0,rate/2000,length(H))';
d=td*ones(size(f));
fm=max(max(f));
m_lim=[0.05 fm -5 5];
d_lim=[0.05 fm 0 8];
M=db(H);
D=gd(H,f);
D(M<-40)=NaN;
subplot(2,1,1)
reset_color(1);
semilogx(f,M)
axis(m_lim)
ylabel('magnitude (dB)')
title('complex IIR filterbank combined transfer function')
subplot(2,1,2)
semilogx(f,D,f,d,':k')
axis(d_lim)
xlabel('frequency (kHz)')
ylabel('delay (ms)')
return
function reset_color(i)
ax=gca;
if (isobject(ax))
ax.ColorOrderIndex=i;
end
return
function y=db(x)
y=20*log10(max(eps,abs(x)));
return
function d=gd(H,f)
[nf,nc]=size(H);
d=zeros(nf,nc);
for k=1:nc
p=unwrap(angle(H(:,k)))/(2*pi);
d(:,k)=-cdif(p)./cdif(f);
end
return
function dx=cdif(x)
n=length(x);
dx=zeros(size(x));
dx(1)=x(2)-x(1);
dx(2:(n-1))=(x(3:n)-x(1:(n-2)))/2;
dx(n)=x(n)-x(n-1);
return
% ffa - fast Fourier analyze real signal
% usage: H=ffa(h)
% h - impulse response
% H - transfer function
function H=ffa(h)
H=fft(real(h));
n=length(H);
m=1+n/2; % assume n is even
H(1,:)=real(H(1,:));
H(m,:)=real(H(m,:));
H((m+1):n,:)=[]; % remove upper frequencies
return