Multirate systems - downsampling

Question

As my teacher explained the process in class, I wanted to write a MATLAB code for learning purposes that downsamples a signal with a downsampling factor of L. Here's the code I wrote and the output that followed:

clear;
clc;

f0 = 5;   % sin() frequency [Hz]
fc1 = 100;   % sampling rate [Hz]
T1 = 1/fc1; 
L = 2;     % downsampling factor
fc2 = fc1/L;    % lower sampling rate [Hz] 
T2 = 1/fc2;
B = fc2/2;  % anti-aliasing filter's band [Hz]
W = 0.1;    % transition band width
N = 100;    % sin() number of samples

% Let's design the anti-aliasing low-pass FIR filter with Parks McClellan
% method
delta = 10^(-7);

amp = [1,0]; 
freq =[B-W*B B]; 
deviaz =[delta, delta];

[N2,normal_freq,amp_resp,W] = firpmord(freq,amp,deviaz,fc1);
hlp=firpm(N2,normal_freq,amp_resp,W);
delay = mean(grpdelay(hlp));    % filter group delay
%fvtool(hlp,1);

% first sampling, at frequency fc1
times1 = (0:N-1)*T1;
samples1 = sin(2*pi*f0*times1);

% filtering the input signal (in this case a sin()) with the anti-aliasing
% filter and compensate for filter delay
samples1_filt = filter(hlp,1,[samples1 zeros(1,delay)]);
samples1_filt = samples1_filt(delay+1:end);


% downsampling by taking one sample each L samples (discarding L-1 samples
% each L samples)
times2 = (0:(N/L)-1)*T2;
samples2 = zeros(1,N/L);

for i = 0:N/L-1
    samples2(i+1) = samples1_filt((i*L)+1);
end


subplot(3,2,1);
stem(times1,samples1,Color='#0072BD');
title(sprintf("fc1 = %d Hz",fc1));


subplot(3,2,2);
stem(times2 ,samples2, Color='#D95319');
title(sprintf("fc2 = %d Hz, L = %d",fc2,L));

subplot(3,2,[3,4]);
hold on
stem(times1 ,samples1);
stem(times2 ,samples2);
hold off
title(sprintf("fc1 (%d Hz) vs fc2 (%d Hz)",fc1,fc2));
legend("fc1","fc2");

subplot(3,2,[5,6]);
hold on
stem(times2 ,samples2, Color='#D95319');
samples2 = decimate(samples1,L,'fir');   % downsampling using built-in function
stem(times2 ,samples2, Color="#77AC30");
hold off
title("comparison with built-in function");
legend("fc2","fc2 built-in");

The program seems to be working, but upon closer inspection, there is a slight difference between the downsampled samples produced by the built-in function and the one that comes from my code. Furthermore, my samples don't start from zero.

So now, I have some questions:

1. Is my code correct? If yes, why do I get these differences with the built-in function? Could it be because of the filter?
2. How should an anti-aliasing filter be designed to achieve the best results? I understand that the filter needs to filter up to $\frac{fc2}{2}$, but does this mean that the band-pass should extend to $\frac{fc2}{2}$ (so the transition band will exceed this boundary), or should both the band-pass and the transition band be within $\frac{fc2}{2}$ like in my code above?
3. How should I handle a non-integer group delay?

Thanks for the help

EDIT 1
With "How should I handle a non-integer group delay?" I mean, for example, choosing delta = 10^(-6); in the filter project I get a low-pass FIR with a group delay of delay = 146.5. Since zeros() only takes integer numbers I can't compensate for the delay introduced by the filter with:

samples1_filt = filter(hlp,1,[samples1 zeros(1,delay)]);
samples1_filt = samples1_filt(delay+1:end);

how should I handle this situation?

I tried to run the same program on a low-pass least-squared FIR filter generated by the following code:

hlp = firls(400,[0 0.45 0.5 1],[1 1 0 0]);

even playing with the filter order and adding some weights I get approximately the same error compared with the built-in function. Why the biggest difference between the built-in function and mine is only on the first sample?

Answer 1

The OP's code appears to be correct and shows a very close match to the built-in function. The differences will exist based on specific filter used, and the initial samples will have errors depending on the initial transient in the filter.

As for "handling" non-integer delays, the delay is based on the filter used and will be what it is. For a linear phase filter (which the OP's filter algorithm used will produce) the delay in samples is $(N-1)/2$ where $N$ is the number of coefficients in the filter. So if a non-integer delay is desired, then an odd-number of filter coefficients should be used. When post processing, a non-causal "zero-phase" forward-reverse filter can be used which is convenient in that no further time alignment is required (although as I described above, if an integer delay is used the alignment is trivial to do). These zero-phase filters are available as filtfilt in MATLAB, Octave and Python's scipy.signal. The filter coefficients are used twice (both in the forward and backward direction), so the filter used is effectively applied twice.

As for understanding the optimal design of decimation filters, decimation involves aliasing just as occurs with A/D conversion. Decimating is the combination of anti-alias filtering and down-sampling. Decimating is digital resampling from a higher rate to a lower rate; A/D conversion is resampling from an infinitely high rate to a lower rate). Just as an anti-alias filter must be used prior to A/D conversion to avoid aliasing distortions, a similar digital anti-alias filter must be used prior to down-sampling. The ideal filter is one that selects the desired signal in a passband with no distortion, and completely rejects all the possible alias regions. This ideal is not achievable, but we can approach it to minimize the distortion to any level depending on the filter complexity and delay.

For example, with decimate by 2 of a real waveform occupying the lower end of the frequency spectrum (low pass signal) from DC to an upper range given as $f_p$, the region in the first Nyquist zone of $f_s/2-f_p$ to $f_s$ prior to down-sampling will alias into the primary spectrum after down-sampling. (Where $f_s$ is the sampling rate prior to down-sampling).

The ideal filter would pass the passband without any distortion and completely reject the stopband, as depicted below. Note how everything else is a "don't care" region. If there was any signal here that could be a distortion, we still have an opportunity down-stream to filter that out, so does not restrict the decimation process.

Where it gets interesting is with higher decimation rates, where efficient multi-band filters can be used (for which the tools readily support with the Parks McClellan method the OP used). These filters will concentrate the rejection only where needed, thus maximizing rejection and minimizing passband ripple for any given filter order. Below shows where the alias regions would be for the decimate by 4 case:

And below further details the passband and multi-stop bands for the decimate by 4:

Note as the decimation ratio increases, the number of regions that can alias into the primary passband increases. The Parks-McClellan filter design has a "equi-ripple" stopband, so the aliasing from each of the stopbands is the same, and thus the rejected noise if flat across the band (white noise) would increase by $10\log_{10}(D)$ in the passband. For this reason, the least-squares optimal filter design algorithm for linear phase filters (firls in MATLAB, Octave and Python's scipy.signal) is preferred as the go-to design for multi-band resampling filters, and will result in the minimum distortion for a filter of the same complexity (number of coefficients) for a signal in the presence of white noise (such as quantization noise or thermal noise).

Note: this topic and many other resampling approaches, is part of an on-line live course I teach "DSP for Wireless Communications" offered several times a year. More details are here.