Polyphase filters for interpolation

In typical digital signal processing applications, there arises need to increase the sampling frequency of a signal sequence, where the higher sampling frequency is an integer multiple of the original sampling frequency i.e for a signal sequence with a sampling frequency , change the sampling frequency to , where is an integer.

Typically, upsampling by a factor of is done by inserting zeros between the original samples and changing the sampling frequency to .

% original sequence
fs = 1000; % sampling frequency of 1kHz
x = cos(2*pi*200*[0:1/fs:10]) + j*sin(2*pi*200*[0:1/fs:10]); % complex sinusoidal with frequency 200Hz
figure
pwelch(x,[],[],[],fs,’twosided’) % psd plot

L = 4; % interpolation factor
xU = [x; zeros(L-1,length(x))];
xU = xU(:).’; % upsampled sequence with zero padding
figure;
pwelch(xU,[],[],[],L*fs,’twosided’) % spectrum of the upsampled sequence

As can be seen from the second figure, the process of inserting zeros in the time domain caused the original spectrum to get repeated i.e. the original sinusoidal present in 200Hz is now present in 200Hz, 1kHz + 200Hz, 2kHz + 200Hz etc…If it is desired to remove the ‘extra aliases of the original spectrum’, we can do so by having a low pass filter following the upsampled sequence.

% lowpass filtering by sinc() filter
h = sinc([-20:20]/L);
y = conv(xU,h);
figure;
pwelch(y,[],[],[],L*fs,’twosided’) % spectrum of the filtered upsampled sequence

With a simple sinc() shaped low pass filter, we cut down the aliases of the orignal spectrum to be 40dB down the desired spectrum.

Having removed the aliases the original spectrum by the low pass filter, comes the important question:

Is it possible to reduce the hardware complexity of the low pass filter implementation considering that the input sequence has zeros in between the samples?

The answer is YES.

To understand how, let us first replace convolution operation by a matrix multiplication where the input sequence is changed to toeplitz form (see previous post).

xUM = toeplitz([xU(1) zeros(1,size(h,2)-1) ], [xU zeros(1,size(h,2)-1) ]);
y1 = h*xUM;
% mean square error
diff = (y1-y)*(y1-y)’/length(y1-y)

This typical direct form filter implementation with coefficients will require multipliers (in this example ). It is reasonably obvious that, due to the presence of the zeros in the input sequence, not all of the multiplier outputs are used. At a particaly time instant, the taps seperated by samples alone has non-zero values. Hence let us reshape h, which is now a vector of dimension into a matrix of dimension .

h = [h zeros(1,L*ceil(length(h)/L) – length(h))]; % padding zeroes
hM = reshape(h,L,length(h)/L);

(We need to pad zeros if the number of coefficients in h () is not an integer multiple of ).

Now, instead of using the toeplitz representation of the upsampled input sequence, let us form the toeplitz representation of the original sequence.

xM = toeplitz([x(1) zeros(1,size(hM,2)-1) ], [x zeros(1,size(hM,2)-1) ]); % toeplitz representation

The matrix multiplication for implementing convolution can be written as,

y2 = hM*xM;

The output y2 hascolumns. Reshaping (i.e reading the first column, second column, third column and so on…) the output y2 for comparison

y2 = y2(:).’;

% mean square error
diff = (y2-y)*(y2-y)’/length(y2-y)

As can be observed, the output y2 is identical to the output of the direct form implementation (y, y1) obtained prior. Instead of single filter with coefficients, we now have filtersets with coefficient each. All the filtersets are fed the input sequence at the original sampling rate . The output from each of the filterset is taken sequentially i.e. output of filtersets #1, #2,…,# , #1, #2, … and so on, are taken at the higher sampling rate . Ofcourse, in hardware we do not need to implement different filtersets. One filterset with taps and dynamically loaded coefficients will do the job.

Summarizing,

1. The original implementation required hardware implementation of a filter with taps. When the modified implementation, the hardware implementation of the filter requires only taps.

2. The original implementation required the input samples to the filter to be clocked at the higher sampling rate of . With the modified implementation, the input samples to the filter is clocked at the original sampling rate of . However, we need to have additional circuitry to dynamically load the coefficients and latch the output at the upsampled frequency .

These filtersets are called polyphase filters. Details about polyphase filters in sufficient mathematical detail is described in Chapter 9.XX in [1]

References:

[1] Digital Signal Processing – Principles, Algorithms and Applications, John G. Proakis, Dimitris G. Manolakis