# Polyphase resampling delay

I want to implement polyphase resampling. For this, I wrote a script in matlab.

The script performs the following sequence of actions:

%% ---- Memory release, closing figures, clear command window

clear; close all; clc;

%% ---- Initialization of algorithm parameters

input = sin(2*pi*50*(0:1/5000:0.3-1/5000)); % Signal generation for resampling
input_length = length(input);

up_rate   = 101; % Upsampling rate
down_rate = 100; % Downsampling rate

max_rate = max(up_rate, down_rate); % Maximum sampling rate

% Filter generation
% Filter order increases linearly with an increase in sampling rate (ord = 20*max_rate),
% because narrower filter required.
% The cutoff frequency is approximately equal to the sampling rate to prevent aliasing
% or suppressing spectrum copies.
hfilter = fir1(20*max_rate, 0.9/max_rate);

filter_length = length(hfilter);
ord = filter_length - 1;

%% ---- Implementation of the logic of the function
%% The generation of the polyphase structure of the filter

% The length of the polyphase filter is equal to the nearest bigger number to the length of the
% original filter that is a multiple of up_rate.
% Missing values will be replaced with zeros
% Thus, the number of phases is integer.
polyphase_filter_length = filter_length;
if rem(filter_length, up_rate)
polyphase_filter_length = (fix(filter_length / up_rate) + 1) * up_rate;
end

polyphase_filter = zeros(1, polyphase_filter_length);
% The length of one phase of the polyphase filter
phase_length = polyphase_filter_length / up_rate;
for i = 0 : up_rate - 1 % Polyphase filter phase cycle
for j = 0 : phase_length - 1 % Cycle along the length of each phase
if (i + j*up_rate < filter_length) % If the length of the original filter is sufficient
% We write down the initial coefficients with decimation into the phase of the
% polyphase filter.
% Recording is done in reverse order to release convolution.
% Decimation coefficients is necessary to simulate the insertion of zeros in the
% original signal.
polyphase_filter(i*phase_length + (phase_length-1-j) + 1) = hfilter(i + j*up_rate + 1);
else  % If the source filter is not long enough
% Write in phase zeros
polyphase_filter(i*phase_length + (phase_length-1-j) + 1) = 0.0;
end
end
end

%% Zero padding the input signal

% Pad the input signal so that the output has from the first to the last valid sample of
% convolution of the signal with the filter.
% So convolving for the first and last samples of the output, the filter intersects with the
% signals in one sample.
for i = 0 : phase_length-1 - 1
end
for i = 0 : input_length-1
padded_input(phase_length-1 + i + 1) = input(i + 1);
end

%% Creating variables for filtering

out_length = fix(up_rate * (input_length + phase_length-1) / down_rate);
if rem(up_rate * (input_length + phase_length-1), down_rate)
out_length = out_length + 1;
end
out = zeros(1, out_length);
acc       = 0;
phase_num = 0;
out_indx  = 0;
i         = 0;
shift     = 0;

%% Filtering

while(i < input_length + phase_length-1)
for j = 0 : phase_length-1
acc = acc + padded_input(i + j + 1) * polyphase_filter(phase_num*phase_length + j + 1);
end
out(out_indx + 1) = acc;
acc = 0;
out_indx = out_indx + 1;

shift = phase_num + down_rate;
i = i + fix(shift / up_rate);
phase_num = rem(shift, up_rate);
end

%% Isolation of the valid part of the filtered signal

% Calculate the length of the zero insert after resampling in front of the original signals
% except for the filter delay.
% By removing this part of the resampling signal, I hope to get the necessary part of the signal
invalid_out_length = fix((up_rate*(phase_length-1) - (ord/2)) / down_rate);
if rem((up_rate*(phase_length-1) - (ord/2)), down_rate)
invalid_out_length = invalid_out_length + 1;
end

res_length = fix(input_length * up_rate / down_rate);
if rem(input_length * up_rate, down_rate)
res_length = res_length + 1;
end
res = zeros(1, res_length);
for i = 0 : res_length-1
res(i + 1) = up_rate * out(i + invalid_out_length + 1);
end

%% ---- Processing the results

res_ref = resample(input, up_rate, down_rate);

subplot(2,1,1);
plot(res);
hold on;
plot(res_ref);
legend('Marlab', 'Script');
title('Resampling signals');
text(5, -0.5, ['uprate = ', num2str(up_rate)]);
text(5, -0.7, ['downrate = ', num2str(down_rate)]);

subplot(2,1,2);
plot(res - res_ref);
title('Resampling error');



But it seems to me that there is a logical error in the script. Maybe in Isolation of the valid part of the filtered signal section. I get two typical oversampling error values relative to the built-in matlab function (in the pictures).

Thus, the question is how to correctly calculate the delay of the polyphase filter after oversampling.