Передескретизация resample

% function y = resampler(x, p, q)
close all; clear all; clc;
 
% yt = load('D:\test\_\tst_rsmpl_polyphase\out.txt');
% plot(yt); grid;
 
x = cos(2*pi*10*(0:1000)/1000) + 1i*sin(2*pi*10*(0:1000)/1000);
p = 3; q = 5;
 
% const
taps = 10;
beta = 5;
% stride = 1;
 
if p==1 && q==1
    y = x;
    return;
end
 
cd = gcd(p, q);
up = p / cd;
down = q / cd;
pqmax = max(up, down);
 
L = 2*taps*pqmax + 1;
% fc = 1/pqmax;
% h = up * firls(L-1, [0 fc fc 1], [1 1 0 0]) .* kaiser(L, beta)';
scale = min(1, up/down);
fc = min(1/up, 1/down)/2;
xx = 2*pi*fc*((0:L-1) - fix((L-1)/2));
h = zeros(1, L);
for i = 1:length(xx)
    if xx(i)==0
        h(i) = 1;
    else
        h(i) = sin(xx(i))/xx(i);
    end
end
h = scale * h .* kaiser(L, beta)';
 
Lhalf = (L-1)/2;
nz = floor(down - mod(Lhalf, down));
z = zeros(1, nz);
h = [z h(:).'];  % ensure that h is a row vector.
Lhalf = Lhalf + nz;
delay = floor(ceil(Lhalf)/down);
% Need to zero-pad so output length is exactly ceil(Lx*p/q).
nz1 = 0;
while ceil(((length(x)-1)*up+length(h)+nz1)/down)-delay < ceil(length(x)*up/down)
    nz1 = nz1+1;
end
h = [h zeros(1,nz1)];
 
%
% y = upfirdn(x, h, up, down);
 
% outLen = fix(length(x) * up / down);
% y = zeros(1, outLen);
% tapsDiv2 = taps / 2;
 
% for idx=0:outLen-1
%     j = fix(idx * down / up);
%     tmp1 = fix(mod(idx*down, up)) + fix((L + 1)/2) - 1;
%     tmp2 = fix((j + 1) * stride);
%     tmp = 0;
%     for k = -tapsDiv2:tapsDiv2-1
%         kx = tmp2 - k*stride;
%         if kx<0; kx = length(x)+kx; end
%         hpos = k * up + tmp1;
%         tmp = tmp + h(1+hpos) * x(1+mod(kx, length(x)));
%     end
%     y(idx+1) = tmp;
% end
 
%
L = length(h);
Mp = fix(L / up);
s = [zeros(1,Mp-1) x zeros(1,Mp)];
y = zeros(1,fix((length(x)+Mp-1)*up/down)+1);
ll = 1;
n = 1;
for i = 1 : length(x)+Mp-1
    for k = 1:up
        if mod(ll+down-1,down) == 0
            sm_tmp = 0;
            for m = 1 : Mp
                sm_tmp = sm_tmp + h((m-1)*up+k) * s(i+Mp-m);
            end
            y(n) = sm_tmp;
            n = n + 1;
        end
        ll = ll + 1;
    end
end
%
 
y(1:delay) = [];
y(ceil(length(x)*p/q)+1:end) = [];
 
figure; plot(real(x)); hold on; plot(real(y),'r'); grid;