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. padded_input_length = input_length + 2*(phase_length-1); padded_input = zeros(1, padded_input_length); for i = 0 : phase_length-1 - 1 padded_input(i + 1) = 0.0; padded_input(padded_input_length-1 - i + 1) = 0.0; 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.