Decimation and Interpolation

Demonstration of noble identities and polyphase decomposition using MATLAB.

1. Department of Digital Signal
Processing
Master of Science in Electronics
Multirate Systems
Homework 1
Decimation and interpolation
Dr. Gordana Jovanovic Dolecek
Ojeda Loredo Fernando
June/15/2015
Sta. Ma. Tonantzintla, Puebla

2. Introduction
The decimator is a device that reduces the sampling rate by an integer factor of M,
whereas the interpolator is used to increase the rate by L.
In many applications the sampling rate of a system needs to be changed for a lower or
higher sampling rate for an appropriate processing of the signal, for example in a digital
mobile receiver system, the received signal often presents a high frequency which cannot be
digitized by the ADC converter , so it needs to be rst converted to a lower frequency using
an Rf analog downconversion, the resultant frequency is commonly known as intermedia-
te frequency, after ADC a FIR lter is used to select a desired bandwidth and then it is
downsampled or upsampled to an appropriate sampling rate. In digital audio applications a
number of sampling rates coexist, for example, the sampling rate for studio work is 48KHz,
whereas for cd production is 44.1KHz, and the broadcast rate is 32KHz.All this rate con-
versions can be made using decimators and interpolators devices, systems that use dierent
sampling rates are known as multirate systems.[3]
Downsampling
The process of downsampling consist in reducing the sampling rate of a signal that
is already sampled, by an integer factor M, basically we are resampling the signal, but
because we are reducing the sampling rate then we call this as subsampling. When a signal
is downsampled we reduce the amount of data by taking only every M-th sample of the
signal and discarding all others as we can see in g 2.
Figure 1: Downsampling.Figure 2: Downsampling in time.
1

3. The downsamplig process can be described in time and frequency in two steps:
step 1. we obtain a signalx
0
(n) fromx(n) by setting all samples whose indexes are not
integer multiples of M to zero. We can considerx
0
(n) as a multiplication ofx(n) andCM(n)
where M denotes the downsamplig factor, in this case we useM= 2 to explain the process.
x
0
(n) =x(n)CM(n) (1)
We express this signal in matlab as follows:
N = 26;
n = 0:N-1;
x = 0.8*sin(2*pi*0.0625*n)+0.3*sin(2*pi*0.2*n);
cm = [zeros(1,26)];
for k = 1:2:26
cm(k) = 1;
end
x2 = x.*cm;
We can see the result of rst step in the gure 3 where only every 2nsamples ofx(n) are
taken and the others are set to zero.
Figure 3: Step 1, time.
This operation does not change the content of the signal x'(n), the spectrum undergoes
a replica every 2=Mas we can see in g. 4, as M = 2 then the spectrum replicas appears
every, where= 1 in normalized frequency. Because we removed samples in time, in
frequency the spectrum is scaled by a factor of 2, we can see this scaling in g. 4.b. Later in
this document we will present the eects of aliasing due to the overlapping of the frequency
bands and how we can avoid this undesired eect that distorts the signal.
2

4. Figure 4: Step 1, frequency.
step 2. The next step is to remove the zeros inx
0
(n) and thus obtainy(m):
y = downsample(x,2);
This operation introduce time scaling by a factor of 1/M. Due to te principle of duality
in time-frequency, in time we get a compression so we expect an expansion in frequency as
we can see in g 6. The frequency is scaled by a factor of 2.
Figure 5: Step 2, scaling in time.Figure 6: Scaling in frequency.
3

5. Aliasing
When a signal is downsampled, the replicas generated by the rst step can overlap, this
occurs if the original signal is not bandlimited to=M. This overlapping eect is called
aliasing, to avoid this we use a low-pass digital lter that should precede the down-sampler,
in g 7 we can see this eect.
Figure 7: Aliasing with M = 4.
Properties of downsampling
A system is said to be linear if it complies whit the principle of superposition
Tfx1(n) +x2(n)g=Tx1(n) +Tx2(n) (2)
We demonstrate this principle using MATLAB as follows:
N = 16;
n = 0:N1;
a = 2;
b = 3;
x1 = 0.8sin (2pi0.0625n ) ;
x2 = 0.3sin (2pi0.2n ) ;
xs = a .x1 + b .x2 ;
xsd = downsample ( xs , 2 ) ; % Structure 1
xd1 = downsample (x1 , 2 ) ;
xd2 = downsample (x2 , 2 ) ;
x12d = a .xd1 + b .xd2 ; % Structure 2
4

6. We observe in g. 8 that xsd and x12 are equal, thus we prove the linearity of the down-
sampling process.
Figure 8: linear property.
We can also demonstrate this principle using Simulink.
Figure 9: Structure 1.Figure 10: Structure 2.
If an input signal undergoes a transformation through a system we get a transformed
output signal, also if the same input signal is shifted in time k samples then we expect the
same transformed output signal but shifted k samples in time, if this is accomplished then
we say that the system is time invariant. we will demonstrate the time variance or time
invariance of a downsampling process as follows:
N = 32;
n = 0:N-1;
k = 4;
x = 0.8*sin(2*pi*0.0625*n);
xs = 0.8*sin(2*pi*0.0625*(n - k));
xd = downsample(x,2);
xsd = downsample(xs,2);
5

7. We observe in g. 12 that the downsampled shifted signal xsd does not present the same
shift as the input signal xs, thus we conclude downsampling is a time varying operation.
Figure 11: Delayed signal by 4 samples.Figure 12: Downsampled delayed signal.
When the relation K/M is not an integer the downsampled delayed signal will have
a fractional delay, this means that the downsampled signal with no input delay and the
downsampled delayed signal will have a dierent shape as we can see in g 13.
Figure 13: Downsampled signals with fractional delay.
6

8. Downsampling identities
First identity
This identity follows directly from the principle of superposition, as previously discussed.
Second identity
Establishes that the delay of M samples before the downsampler is equivalent to a delay
of one sample after the downsampler. We will show this second identity using Simulink, as
follows:
Figure 14: Structure 1.Figure 15: Structure 2.
Comparing the two graphics in g. 16 and g 17 it follows that structure 1 and structure
2 are equivalent.
Figure 16: Structure 1 plot.Figure 17: Structure 2 plot.
7

9. Third identity
This identity is related to cascade connection of a linear time invariant systemH(Z) and
a downsampler. Filtering withH(Z
M
) and downsampling byMis equal to the downsampling
byMand ltering withH(Z).
Figure 18: Third identity structures.
The third identity is demonstrated using a lter expanded with M = 4 which precedes
the downsampler as can be seen in structure 1 in gure 18. Another lter is designed to
implement structure 2, and nally the graphs in gure 19, show this identity, the Matlab
implementation is as follows:
n = 0 : 6 0 ; % Time index
x = cos (2pi0.05n ) ; % Generating the o r i g i n a l s i g n a l
h = f i r 1 ( 1 0 , 0 . 5 ) ; % Designing the f i l t e r t r a n s f e r function H( z )
hu = upsample (h , 4 ) ; % Transfer function H( z^M)
y1 = f i l t e r (hu ,1 , x ) ; % F i l t e r i n g
y = downsample (y1 , 4 ) ; % Downsampling
m = 0: length (y)1; % Time index
f i g u r e (1)
subplot (3 ,2 ,1) , stem (n , x ) , ylabel ( ' x [ n ] ' )
t i t l e ( ' Figure 18( a ) ')
subplot (3 ,2 ,3) , stem (n , y1 ) , ylabel ( ' y1 [ n ] ' )
subplot (3 ,2 ,5) , stem (m, y ) , ylabel ( ' y [m] ' )
xlabel ( ' Time index ' )
y2 = downsample (x , 4 ) ; % down sampling
y = f i l t e r (h ,1 , y2 ) ; % f i l t e r i n g
subplot (3 ,2 ,2) , stem (n , x , ' r ' ) , ylabel ( ' x [ n ] ' )
t i t l e ( ' Figure 18(b ) ')
subplot (3 ,2 ,4) , stem (m, y2 , ' r ' ) , ylabel ( ' y2 [m] ' )
subplot (3 ,2 ,6) , stem (m, y , ' r ' ) , ylabel ( ' y [m] ' )
xlabel ( ' Time index ' )
In gure 19 the left-hand side shows the signals for structure 1, and the right-hand side
presents the signals for structure 2. This results shown in g. 19 demonstrate the equivalence
of the cascade connections dened by the third identity.
8

10. Figure 19: Illustration of the third identity.
Polyphase decimation
The decimation structure consists of two block as can be seen in gure 20, a low-pass
lter which discard all frequencies above=Mto avoid aliasing, and the downsamplign block
which reduce the sampling rate of the signal.
Figure 20: Decimation structure.
UsuallyH(z) is a FIR lter, which consist of N coecients. If N is a integer multiple of
M we can obtain M dierent discretely sampled components of the impulse response. That
means that the corresponding Z-transform ca be partitioned into M sub-signals
H(z) =
M1
X
k=0
z
k
Hk(z
M
) (3)
where
Hk(z
M
) =
N=M1
X
n=0
h(nM+k)(z
M
)
n
) (4)
9

11. The lter shown in the structure from gure 20, can be decomposed as we can observe in
gure 21.a this is a direct application of the rst identity, now in gure 21.b we can see an
application of the third identity. With this implementation we obtain a more ecient version
in which both the number of lter operations as well the amount of memory required are
reduced by a factor of M.
Figure 21: Polyphase decimation.
Now we present an example in Matlab with a FIR lter which consist of N = 64 coe-
cients, an a decimation factor of M = 4, so we expect 4 polyphase components.
% Polyphase decomposition
c l e a r all , c l o s e a l l
% Input s i g n a l
n = 0 : 6 3 ;
h = zeros ( s i z e (n ) ) ; h (11:39) = 0 . 9 5 . ^ ( 1 : 2 9 ) ; % Generating the sequence 'h '
% Polyphase downsampling with the phase of f s e t
10

12. h0 = downsample (h , 4 ) ;
h1 = downsample (h , 4 , 1 ) ;
h2 = downsample (h , 4 , 2 ) ;
h3 = downsample (h , 4 , 3 ) ;
subplot (4 ,1 ,1) , stem ( 0 : length ( h0)1,h0 ) , ylabel ( ' h0 [m] ' )
subplot (4 ,1 ,2) , stem ( 0 : length ( h1)1,h1 ) , ylabel ( ' h1 [m] ' )
subplot (4 ,1 ,3) , stem ( 0 : length ( h2)1,h2 ) , ylabel ( ' h2 [m] ' )
subplot (4 ,1 ,4) , stem ( 0 : length ( h3)1,h3 ) , ylabel ( ' h3 [m] ' )
xlabel ( ' Time index m' )
The original signal is shown in gure 22 and the polyphase components are given in 23.
Figure 22: Filter Coecients.Figure 23: Polyphase components.
Upsampling
Upsampling increases the sampling rate by an integer factor L, by inserting L-1 equally
spaced zeros between each pair of samples.
y[m] =
(
x[m=L]; m= 0;L;2L; :::
0; Otherwise
(5)
Where L is called an interpolator factor. The symbol for this operation is a box with an
upward-pointing arrow, followed by the interpolator factor.[2]
Figure 24: Upsampler.
11

13. The following Matlab code illustrates an example of the upsampler with L = 2.
N = 21; % Length of the o r i g i n a l sequence
n=0:N1; % Time index
x=0.7sin (2pi0.0625n)+0.3sin (2pi0.2n ) ; % Original s i g n a l
L= 2; % Upsampling f ac t o r
y=zeros (1 ,Llength (x ) ) ;
y ( [ 1 : L: length (y)])= x ; % Up sampled s i g n a l
Ny = length (y)L+1; % Length of the up sampled s i g n a l
f i g u r e (2)
subplot (2 ,1 ,1)
stem (0:20 , x ( 1 : 2 1 ) )
xlabel ( ' Time index n ' ) , ylabel ( ' x [ n ] ' )
subplot (2 ,1 ,2)
stem (0:40 , y ( 1 : 4 1 ) )
xlabel ( ' Time index m' ) , ylabel ( ' y [m] ' )
We can see in gure 25, how the sampling rate is increased also we can see how L-1
equally space zeros are inserted between each pair of samples.
Figure 25: upsampling by L = 2.
12

14. Imaging
The interpolator consist of two stages as we can see in gure 26, the upsampler block
that was previously discussed and a low pass lter that is after the upsampler.
Figure 26: Interpolator
The process of upsampling introduces the replicas of the main spectra at every 2=L.
This is called imaging, since there are L-1 replicas (images) in 2. We illustrate this eect
next:
F = [ 0 , 0 . 2 , 0 . 9 , 1 ] ;
A = [ 0 , 1 , 0 , 0 ] ;
x = f i r 2 (128 ,F,A) ;
X = f f t (x , 1 0 2 4 ) ;
f = 0:1/1024:(5121)/1024;
y1 = upsample (x , 2 ) ;
y2 = upsample (x , 3 ) ;
y3 = upsample (x , 4 ) ;
Y1 = f f t (y1 , 1 0 2 4 ) ;
Y2 = f f t (y2 , 1 0 2 4 ) ;
Y3 = f f t (y3 , 1 0 2 4 ) ;
subplot (4 ,1 ,1) , plot (2f , abs (X( 1 : 5 1 2 ) ) ) , t i t l e ( ' o r i g i n a l signal ' )
subplot (4 ,1 ,2) , plot ( f2 , abs (Y1( 1 : 5 1 2 ) ) ) , t i t l e ( ' upsampled s i g n a l L = 2 ')
subplot (4 ,1 ,3) , plot (2f , abs (Y2( 1 : 5 1 2 ) ) ) , t i t l e ( ' upsampled s i g n a l L = 3 ')
subplot (4 ,1 ,4) , plot ( f2 , abs (Y3( 1 : 5 1 2 ) ) ) , t i t l e ( ' upsampled s i g n a l L = 4 ')
In gure 27 we can observe the dierent replicas for L = 2,3,4. We use the low pass lter
to remove these replicas, this lter is called an anti-imaging lter. In Matlab we can show
the use of this lter by using the command interp, like the command decimate, this already
includes the lter unlike the commands downsample and upsample which do not include it.
13

15. Figure 27: Interpolator
Next code show the use ofinterp, to remove the replicas for L = 4.
F = [ 0 , 0 . 2 , 0 . 9 , 1 ] ;
A = [ 0 , 1 , 0 , 0 ] ;
x = f i r 2 (128 ,F,A) ;
X = f f t (x , 1 0 2 4 ) ;
f = 0:1/1024:(5121)/1024;
L = 4;
xu = upsample (x ,L ) ;
y = interp (x ,L ) ;
Xu = f f t (xu , 1 0 2 4 ) ;
Y = f f t (y , 1 0 2 4 ) ;
subplot (3 ,1 ,1) , plot (2f , abs (X( 1 : 5 1 2 ) ) ) , t i t l e ( ' o r i g i n a l signal ' )
subplot (3 ,1 ,2) , plot ( f2 , abs (Xu( 1 : 5 1 2 ) ) ) , t i t l e ( ' upsampled s i g n a l L = 4 ')
subplot (3 ,1 ,3) , plot (2f , abs (Y( 1 : 5 1 2 ) ) ) , t i t l e ( ' interpolated signal ' )
f i g u r e
subplot (3 ,1 ,1) , stem (54:74 , x ( 5 5 : 7 5 ) ) ;
subplot (3 ,1 ,2) , stem(4551:4741,xu(4551:4741));
subplot (3 ,1 ,3) , stem(4551:4741,y(4551:4741));
14

16. In the following gure we observe the eect of imaging also is necessary to say that the
upsampling process consist of the increase in the sampling rate of the signal, this is known
as an expansion in time, by the principle of duality, in the frequency domain we expect a
compression as shown in gure 28.
Figure 28: Imaging and interpolation in frequency
In time something very interesting happens due to the lter, the zeros that are inserted
by the upsampler are interpolated, for this reason we call this lter an interpolation lter, we
can see it in gure 29, also the frequency spectrum is scaled by a factor of L this is presented
in gure 28.c.
Figure 29: Interpolation in time
15

17. Properties of upsampling
Now we show the linearity of the upsample operation, just like in downsampling we will
use the same code but now we use the command upsample instead of downsample as shown
next.
N = 16;
n = 0:N1;
a = 2;
b = 3;
x1 = 0.8sin (2pi0.0625n ) ;
x2 = 0.3sin (2pi0.2n ) ;
xs = a .x1 + b .x2 ;
xsd = upsample ( xs , 2 ) ; % Structure 1
xd1 = upsample (x1 , 2 ) ;
xd2 = upsample (x2 , 2 ) ;
x12d = a .xd1 + b .xd2 ; % Structure 2
We observe in gure 30 that xsd and xsd12 are equal so we can conclude that upsampling
is a linear operation.
Figure 30: Upsampling linear property
16

18. We also demonstrate this property using Simulink as follows:
Figure 31: Structure 1Figure 32: Structure 2
Now we will show the time varying property of the upsampling block. We suppose that
the input of the upsampler has a delay ofDsamples.
x(mD) (6)
The upsampled signal will be:
y((mD)L) =y(mLDL) =y(nDL)6=y(nD) (7)
Consequently, upsampling is a time-dependent operation.
N = 32;
n = 0:N1;
k = 9;
x = 0.8sin (2pi0.0625n ) ;
xs = 0.8sin (2pi0.0625(nk ) ) ;
xu = upsample (x , 2 ) ;
xsu = upsample ( xs , 2 ) ;
f i g u r e
subplot (2 ,1 ,1) , stem (x , ' LineWidth ' , 2 ) , grid , t i t l e ( ' x ' ) , grid ;
subplot (2 ,1 ,2) , stem ( xs , ' LineWidth ' , 2 ) , grid , t i t l e ( ' xs ' ) , grid ;
f i g u r e
subplot (2 ,1 ,1) , stem (xd , ' r ' , ' LineWidth ' , 2 ) , grid , t i t l e ( ' xu ' ) , grid ;
subplot (2 ,1 ,2) , stem ( xsd , ' r ' , ' LineWidth ' , 2 ) , grid , t i t l e ( ' xsu ' ) , grid ;
In gure 33 xs represent a shifted version of x by 9 samples, in gure 34 xu is the
upsampled signal by L = 2 of x and xsu is the upsampled signal of xs. We expect the same
amount of delay in xsu and xs, but xsu is shifted by 14 samples unlike xs which is shifted
by 9 samples, thus we prove that upsampling is a time-varying operation.
17

19. It is important to mention that unlike downsampling, the upsampled delayed signal and
the upsampled signal with no delay will always have the same shape, this is because the
delayDLis always an integer.
Figure 33: Delayed signal by 9 samplesFigure 34: Upsampled delayed signal
Upsampling identities
We have already seen three useful identities of the downsampled signals, and now we will
state the corresponding identities associated with upsampling.
Fourth identity
We have already demonstrated the fourth identity as it follows from the principle of
superposition.
Fifth identity
This identity states that a delay of one sample before upsampling is equivalent to the
delay of L samples after upsampling.
Figure 35: Structure 1Figure 36: Structure 2
18

20. Figure 37 and 38 are equal, so the equivalence of the two structures is conrmed.
Figure 37: Structure 1 plotFigure 38: Structure 2 plot
Sixth identity
States that the ltering followed by upsampling is equivalent to having upsampling rst
followed by expanded ltering.
Figure 39: Sixth identity structures
n = 0 : 1 5 ; % Time index
x = cos (0.2pin ) ; % Generating the o r i g i n a l s i g n a l
h = f i r 1 ( 1 0 , 0 . 5 ) ; % Designing the f i l t e r t r a n s f e r function H( z )
hu = upsample (h , 2 ) ; % Transfer function H( z^L)
y1 = f i l t e r (h ,1 , x ) ; % F i l t e r i n g
y = upsample (y1 , 2 ) ; % Up sampling
m = 0: length (y)1; % Time index
f i g u r e (1)
subplot (3 ,2 ,1) , stem (n , x ) , ylabel ( ' x [ n ] ' )
subplot (3 ,2 ,3) , stem (n , y1 ) , ylabel ( ' y1 [ n ] ' )
subplot (3 ,2 ,5) , stem (m, y ) , ylabel ( ' y [m] ' )
xlabel ( ' Time index ' )
axis ([0 ,30 ,1 ,1])
y2 = upsample (x , 2 ) ; % Up sampling
y = f i l t e r (hu ,1 , y2 ) ; % F i l t e r i n g
19

21. subplot (3 ,2 ,2) , stem (n , x , ' r ' ) , ylabel ( ' x [ n ] ' )
subplot (3 ,2 ,4) , stem (m, y2 , ' r ' ) , ylabel ( ' y2 [m] ' )
axis ([0 ,30 ,1 ,1])
subplot (3 ,2 ,6) , stem (m, y , ' r ' ) , ylabel ( ' y [m] ' )
xlabel ( ' Time index ' )
axis ([0 ,30 ,1 ,1])
In gure 40 the left-hand side shows the signals for structure 1, and the right-hand side
presents the signals for structure 2. This results shown in g. 40 demonstrate the equivalence
of the cascade connections dened by the sixth identity.
Figure 40: Illustration of the sixth identity
Polyphase interpolation
The convolution at the higher sampling rate can be replaced by independent convolutions
at the lower input sampling rate using polyphase decomposition.[1]
H(z) =
L1
X
k=0
z
k
Hk(z
L
) (8)
In gure 41 we can see the polyphase components of the signal.
20

22. Next we will show the implementation of the sixth identity in Matlab.
n = 0 : 6 3 ;
h = zeros ( s i z e (n ) ) ; h (11:39) = 0 . 9 5 . ^ ( 1 : 2 9 ) ; % F i l t e r c o e f f i c i e n t s
hd0 = downsample (h , 4 ) ;
hd1 = downsample (h , 4 , 1 ) ;
hd2 = downsample (h , 4 , 2 ) ;
hd3 = downsample (h , 4 , 3 ) ;
% Upsampling polyphase components with the phase o f f s e t
h0 = upsample (hd0 , 4 ) ;
h1 = upsample (hd1 , 4 , 1 ) ;
h2 = upsample (hd2 , 4 , 2 ) ;
h3 = upsample (hd3 , 4 , 3 ) ;
subplot (4 ,1 ,1) , stem ( 0 : length ( h0)1,h0 ) , ylabel ( ' h0 [m] ' )
subplot (4 ,1 ,2) , stem ( 0 : length ( h1)1,h1 ) , ylabel ( ' h1 [m] ' )
subplot (4 ,1 ,3) , stem ( 0 : length ( h2)1,h2 ) , ylabel ( ' h2 [m] ' )
subplot (4 ,1 ,4) , stem ( 0 : length ( h3)1,h3 ) , ylabel ( ' h3 [m] ' )
xlabel ( ' Time index m' )
Figure 41: Polyphase components
21

23. References
[1] Jovanovic-Dolecek, G. (2001).Multirate Systems: Design and Applications: Design and
Applications. IGI Global.
[2] Milic, L. (2009).Multirate Filtering for Digital Signal Processing: MATLAB Applications:
MATLAB Applications. IGI Global.
[3] Vaidyanathan, P. P. (1993).Multirate systems and lter banks. Pearson Education India.
22