0
$\begingroup$

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).

enter image description here

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.