[249cf9]: comp / comp_pfilt.m Maximize Restore History

Download this file

comp_pfilt.m    69 lines (48 with data), 1.4 kB

 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
function h=comp_pfilt(f,g,a,do_time);
%COMP_PFILT Compute filtering
%
% If *do_time* is 1, the routine will use the time-side algorithm for
% FIR filters, otherwise it will always do the multiplication in the
% frequency domain.
[L,W]=size(f);
if numel(a)==1
N=L/a;
else
N=L/a(1)*a(2);
afrac=a(1)/a(2);
end;
h=zeros(N,W,assert_classname(f));
l=(0:L-1).'/L;
realoutput = isreal(f) && isfield(g,'h') && isreal(g.h) && g.fc==0;
if isfield(g,'h') && do_time
% Use a direct algorithm
g_time=circshift(postpad(g.h,L),g.offset).*...
exp(2*pi*1i*round(g.fc*L/2)*l);
if g.realonly
g_time=real(g_time);
end;
g_time=conj(involute(g_time));
for n=0:N-1
h(n+1,:)=sum(bsxfun(@times,f,circshift(g_time,a*n)));
end;
else
% Zero-extend and use a full length fft algorithm
% This case can be further optimized
G=comp_transferfunction(g,L);
if numel(a)==1
for w=1:W
F=fft(f(:,w));
h(:,w)=ifft(sum(reshape(F.*G,N,a),2))/a;
end;
else
Llarge=ceil(L/N+1)*N;
amod=Llarge/N;
for w=1:W
F=fft(f(:,w));
h(:,w)=ifft(sum(reshape(fir2long(F.*G,Llarge),N,amod),2))/afrac;
end;
end;
end;
if realoutput
h=real(h);
end;