1 Star2 Stars3 Stars4 Stars5 Stars (6 votes, average: 3.67 out of 5)
Loading ... Loading ...
Print Print

Polyphase filters for interpolation

by Krishna Sankar 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 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

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.

{ 8 comments… read them below or add one }

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,

Reply

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?

Reply

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!

Reply

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.

Reply

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?

Reply

Krishna Pillai 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.

Reply

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!

Reply

Krishna Pillai June 21, 2009 at 12:20 pm

@sumo: Nice question. Let me try to answer…
With the original sampling frequency of 1000Hz, we were able to ‘see’ frequencies from [-500Hz to +500Hz). Now, with the oversampled frequency of 4000Hz, we now can ‘see’ frequencies from [-2000Hz to 2000Hz). Now as the specrum gets replicated at multiples for sampling frequency, the frequency at 200Hz is replicated at 1200Hz, -800Hz, -1800Hz. So, instead of 1 frequency, we have four frequencies at 1/4th amplitude. Hence the reduction in magnitude by 20*log10(1/4) = -12dB.

Do you agree?

Reply

Leave a Comment

Previous post:

Next post: