(6 votes, average: 3.67 out of 5)

# Polyphase filters for interpolation

by on May 12, 2007

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 $x(n)$ with a sampling frequency $f_s$, change the sampling frequency to $Lf_s$, where $L$ is an integer.

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

% 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 $L-1$ zeros in between the samples?

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 $N$ coefficients will require $N$ multipliers (in this example $N=41$). 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 $L$samples alone has non-zero values. Hence let us reshape h, which is now a vector of dimension $[1\ \mbox{x}\ N]$ into a matrix of dimension $[\ L\ \mbox{x}\ \lceil \frac{N}{L} \rceil \ ]$.

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 ($N$) is not an integer multiple of $L$).

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 has$L$columns. 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$N$ coefficients, we now have $L$filtersets with $\lceil \frac{N}{L} \rceil$ coefficient each. All the $L$ filtersets are fed the input sequence at the original sampling rate $f_s$. The output from each of the filterset is taken sequentially i.e. output of filtersets #1, #2,…,#$L$ , #1, #2, … and so on, are taken at the higher sampling rate $Lf_s$ . Ofcourse, in hardware we do not need to implement $L$different filtersets. One filterset with$\lceil \frac{N}{L} \rceil$ taps and dynamically loaded coefficients will do the job.

Summarizing,

1. The original implementation required hardware implementation of a filter with $N$taps. When the modified implementation, the hardware implementation of the filter requires only $\lceil \frac{N}{L} \rceil$ taps.

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

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

References:

D id you like this article? Make sure that you do not miss a new article by subscribing to RSS feed OR subscribing to e-mail newsletter. Note: Subscribing via e-mail entitles you to download the free e-Book on BER of BPSK/QPSK/16QAM/16PSK in AWGN.

Talib October 8, 2009 at 11:59 pm

Also Krishna,

Shouldnt:

diff = (y1-y)*(y1-y)’/length(y1-y)

actually be…

diff = (y1-y)*(y1-y)’/length(y1) ?

Thanks again,

Krishna Sankar October 12, 2009 at 5:36 am

@Talib: It really does not matter, no? length(y1-y) and length(y1) reports the same number. Agree?

Talib October 8, 2009 at 10:59 pm

Hi Krishna,

Great article on polyphase – question though – could you go through in more detail about EXACTLY how you made this sinc-filter? What is its sampling rate?… Why is ‘L’ an argument? … I understand that you need to low pass filter, and that a sinc in time is a rectangle in frequency, but exactly how did you select its arguments?…

Thanks!

Krishna Sankar October 12, 2009 at 5:34 am

@Talib: My replies:
1/ In the code, I assumed a sampling rate of 1kHz
2/ L is the oversampling factor. I used a small matlab code snippet to plot the frequency response
octave:12> L = 4
octave:13> h = sinc([-20:20]/L);
octave:14> hF = fft(h,1024);
octave:15> plot([-512:511]/1024,(abs(fftshift(hF))));
octave:16> xlabel(‘freq, kHz’); ylabel(‘amplitude’);

Hope this helps.

Saria July 5, 2009 at 7:33 pm

Dear Mr. Krishna, thank you for the wonderful explanation of polyphase filters. I would like to know how you chose the sinc filter time range [-20:20]/L,which will give a vector of 41 values in time domain?

Krishna Sankar July 6, 2009 at 7:25 pm

@saira: No special reason, I just wanted to define a filter which provided around 40dB attenuation outside the passband.

sumo June 17, 2009 at 4:06 pm

why after upsampling, as your example L = 4,
the magnitude of 200Hz is decreased by 20*log10(1/4), which is about -12dB.

Do you know how to explain this case ?

Thanks!