52 lines
801 B
Matlab
Executable File
52 lines
801 B
Matlab
Executable File
% pfftfilt(x,b)
|
|
function y = pfftfilt(x,h,P,S,L)
|
|
|
|
N = length(h);
|
|
Lx = length(x);
|
|
C = 2^ceil(log(L+S*L-1)/log(2));
|
|
Np = S*L;
|
|
K = Lx/L;
|
|
|
|
y = zeros(Lx,1);
|
|
Hp = zeros(C,P);
|
|
Xkps = zeros(C,S*P);
|
|
yk = zeros(C,1);
|
|
Y = zeros(C,1);
|
|
xzp = zeros(Lx+C-L,1);
|
|
pidTbl = [2:P*S];
|
|
pidTbl(P*S) = 1;
|
|
pStbl = [0:P*S-1];
|
|
pStbl(1) = 10;
|
|
|
|
pid = 1;
|
|
|
|
% h(k) => Hp(k)
|
|
for p=1:P,
|
|
Hp(1:C,p) = fft([h((p-1)*Np+1:p*Np);zeros(C-Np,1)],C);
|
|
end;
|
|
|
|
xzp(C-L+1:Lx+C-L) = x;
|
|
|
|
for k=1:K,
|
|
|
|
kL = (k-1)*L;
|
|
|
|
Xkps(1:C,pid) = fft(xzp(kL+1:kL+C),C);
|
|
|
|
pS = pid;
|
|
Y = Hp(1:C,1).*Xkps(1:C,pid);
|
|
for p=2:P,
|
|
for i = 1:S
|
|
pS = pStbl(pS);
|
|
end;
|
|
Y = Y + Hp(1:C,p).*Xkps(1:C,pS);
|
|
end;
|
|
|
|
yk = ifft(Y,C);
|
|
|
|
y(kL+1:kL+L) = yk(C-L+1:C);
|
|
|
|
pid = pidTbl(pid);
|
|
|
|
end;
|