Untitled

by Ljiljana Milic

Supplemental material for Chapter IV

4. FIR Filters for Sampling Rate Conversion

Table of Contents

4. FIR Filters for Sampling Rate Conversion 4.1 Introduction 4.2 Structure Verification and Simulation Using MATLAB Example 4.1 Example 4.2 4.3 Memory Saving Structures for FIR Polyphase Decimators and Interpolators Example 4.3 4.4 Computational Efficiency of Single-Stage and Multistage Systems Example 4.4

4.1 Introduction

The transfer function of an FIR filter is a the polynomial in

, and is usually writen in the form

Here, we express

by means of M (or L) polyphase branches

where

Based on the polyphase implementation, operations in sampling-rate conversion systems can be executed at the lower sampling rate. The efficient commutative implementation structures for polyphase decimator| and interpolator are shown and verified in the next subsections.

4.2 Structure Verification and Simulation Using MATLAB

Commutative polyphase structure of a decimator.

Commutative polyphase structure of an interpolator.

Example 4.1

Decimation and interpolation by integer factors M = L = 5
Design of optimal linear-phase FIR filter
Creation of the input signal
Polyphase decomposition to 5 polyphase components
Polyphase decimation according to the Figure shown above
Polyphase interpolation according to the Figure shown above
Structure verification

clear all, close all

STEP 1: Lowpass FIR filter design to satisfy Case a tolerance scheme.

Filter specifications: Fp = 0.08, Fs = 1/M - Fp, Passband ripple ap = 0.01 dB, minimum stopband

attenuation as = 60 dB.Specifications are met with an optimal FIR filter of the length N = 60.

h = firpm(59,[0,2*400/10000,1000/5000,1],[1,1,0,0]); % Filter coefficients
[H,F1] = freqz(h,1,1024,10000);  % Computation of filter frequency response
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))), legend('Filter magnitude response'),...
    ylabel('Gain, dB'),axis([0,5000,-80,5])

STEP 2: Generating the input signal

The input signal is a causal sequence formed as a sum of three sinusoidal sequencies corrupted in noise.

where

is random signal with the normal distribution.

The sampling frequency at the input

. The sequence length N = 511 samples.

n = 0:510; % Time index
x = cos(2*pi*200*n/10000)+cos(2*pi*400*n/10000)+0.7*cos(2*pi*3800*n/10000)...
    +randn(size(n));
[X,F1] = freqz(x,1,1024,10000); % Computing the signal spectrum
subplot(4,1,2), plot(F1,abs(X)),... 
legend('Spectrum: input signal'), ylabel('|X(F)|'), axis([0,5000,0,300])

STEP 3: Polyphase decomposition

E = reshape(h,5,length(h)/5);

STEP 4: Polyphase downsampling and filtering

x0 = x(1:5:length(x)); x1 = [0,x(5:5:length(x))]; x2 = [0,x(4:5:length(x))];
x3 = [0,x(3:5:length(x))]; x4 = [0,x(2:5:length(x))];
y0 = filter(E(1,:),1,x0); y1 = filter(E(2,:),1,x1); y2 = filter(E(3,:),1,x2);...
    y3=filter(E(4,:),1,x3); y4=filter(E(5,:),1,x4);
ydec = y0 + y1 + y2 + y3 + y4; % Decimated signal
[Ydec,F2] = freqz(ydec,1,512,2000);
subplot(4,1,3), plot(F2,abs(Ydec)),legend('Spectrum: decimated signal'),...
     ylabel('|Y_d_e_c(F)|'), axis([0,1000,0,60])

STEP 5: Poliphase filtering and up-sampling

x = ydec; % Input to the interpolator

u0 = 5*filter(E(1,:),1,x); u1 = 5*filter(E(2,:),1,x); u2 = 5*filter(E(3,:),1,x);

u3 = 5*filter(E(4,:),1,x); u4 = 5*filter(E(5,:),1,x);

U = [u0;u1;u2;u3;u4];

yint = U(:); % interpolated signal

[Yint,F3] = freqz(yint,1,1024,10000);

subplot(4,1,4), plot(F3,abs(Yint)), legend('Spectrum: interpolated signal'),...

xlabel('Frequency, Hz'), ylabel('|Y_i_n_t(F)|')

axis([0,5000,0,300])

disp('END OF EXAMPLE 4.1')

END OF EXAMPLE 4.1

Example 4.2

Decimation and interpolation by integer factors M = L = 5
Design of optimal linear-phase FIR filter
Creation of the input signal
Polyphase decomposition to 5 polyphase components
Simulation of real-time decimation
Simulation of real-time interpolation.

clear all, close all

STEP 1: Lowpass FIR filter design to satisfy Case a tolerance scheme.

Filter specifications from example 4.1.

h = firpm(59,[0,400/5000,1000/5000,1],[1,1,0,0]); % Filter coefficients
[H,F1] = freqz(h,1,1024,10000); % Computing the signal spectrum
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))),...
    legend('Filter magnitude response'), ylabel('Gain, dB'),axis([0,5000,-80,5])

STEP 2: Generating the input signal as given in Example 4.1

n = 0:510;  % Time index
x = cos(2*pi*200*n/10000)+cos(2*pi*400*n/10000)+0.7*cos(2*pi*3800*n/10000)...
    +randn(size(n));
[X,F1] = freqz(x,1,1024,10000); % Computing the spectrum of the input signal
subplot(4,1,2)
plot(F1,abs(X)), legend('Spectrum: input signal'), ylabel('|X(F)|')
axis([0,5000,0,300])
% STEP 3
% Polyphase decomposition

Polyphase decomposition

E = reshape(h,5,length(h)/5);

STEP 4: Simulation of decimator

vk_i = zeros(size(E)); %Initial states of polyphase components
rot = 0; % Initial rotator position
ydec = [];
for n=1:length(x)
xn = x(n);
vk = vk_i(rot+1,:); ek = E(rot+1,:);
vk = [xn,vk(1:length(vk)-1)]; yn = sum(vk.*ek); % filtering
vk_i(rot+1,:)=vk;
if n==1, yy=0; else, end
if rot == 0 ss = 4; yy = yy+yn; ydec = [ydec,yy]; else end
if rot == 1 ss = 0; yy = yy+yn; else end
if rot == 2 ss = 1; yy = yy+yn; else end
if rot == 3 ss = 2; yy = yy+yn; else end
if rot == 4 ss = 3; yy = yn; else end
rot = ss;
end
[Ydec,F2] = freqz(ydec,1,512,2000);
subplot(4,1,3), plot(F2,abs(Ydec)),...
    legend('Spectrum: decimated signal'), ylabel('|Y_d_e_c(F)|')
axis([0,1000,0,60])

STEP 5: Simulation of interpolator

yi = [0]; xk_i = zeros(size(E));

rot = 0;

k = 1;

for n=1:length(ydec)*5

xn=ydec(k);

xk = xk_i(rot+1,:); ek = E(rot+1,:);

xk=[xn,xk(1:length(xk)-1)]; yn=sum(xk.*ek); % filtering

xk_i(rot+1,:)=xk; yn=5*yn;

if rot == 0 ss = 1; yi = [yi,yn]; else end

if rot == 1 ss = 2; yi = [yi,yn]; else end

if rot == 2 ss = 3; yi = [yi,yn]; else end

if rot == 3 ss = 4; yi = [yi,yn]; else end

if rot == 4 ss = 0; k = k+1; yi = [yi,yn]; else end

rot = ss;

end

[Yint,F3] = freqz(yi,1,1024,10000);

subplot(4,1,4)

plot(F3,abs(Yint)); legend('Spectrum: interpolated signal')

xlabel('Frequency, Hz'),ylabel('|Y_i_n_t(F)|')

axis([0,5000,0,300])

disp('END OF EXAMPLE 4.2')

END OF EXAMPLE 4.2

4.3 Memory Saving Structures for FIR Polyphase Decimators and Interpolators

Memory saving structures based on polyphase implementation of decimator and interpolator are shown bellow for sampling-rate conversion factor M = L = 3.

Memory saving structure of the polyphase factor-of-3 decimator.

Memory saving structurefor the polyphase factor-of-3 interpolator.

Example 4.3

This exampl is developed to simulate the memory saving commutative decimator and interpolator for the sampling-rate conversion factor M = L = 3. The sampling-frequency of the test signal at the input amounts 10000 Hz.

Program is organized in 5 steps.

clear all, close all

STEP 1: Lowpass FIR filter design:

Optimal FIR filter of the length N = 12. Specifications: Case b tolerance scheme.

h = firpm(11,[0,800/5000,2*1666/5000-800/5000,1],[1,1,0,0]);
[H,F1] = freqz(h,1,1024,10000);
figure 
subplot(4,1,1), plot(F1,20*log10(abs(H))), legend('Filter magnitude response'), ylabel('Gain, dB')
axis([0,5000,-60,5])

STEP 2: Polyphase decomposition

E = reshape(h,3,length(h)/3); E = fliplr(E);

STEP 3: Input signal

n = 0:510; % Time index
x = cos(2*pi*350*n/10000) + cos(2*pi*750*n/10000) + ...
    0.7*cos(2*pi*3800*n/10000) + randn(size(n));
[X,F1] = freqz(x,1,1024,10000); % Computing the spectrum of the input signal
subplot(4,1,2)
plot(F1,abs(X)), legend('Spectrum: input signal'), ylabel('|X(F)|')
axis([0,5000,0,300])

STEP 4

Setting the initial states for decimator

xk = zeros(size(1:length(h)/3));

Decimation

rot = 0; % Initial rotator position
ydec = [];
for n=1:length(x)
xn = x(n); xk = xn*E((rot+1),:)+xk;
if rot == 0 ss = 2; yn = xk(length(xk)); ydec = [yn,ydec];...
        xk = [0,xk(1:length(xk)-1)]; else end
if rot == 1 ss = 0; else end
if rot == 2 ss = 1; else end
rot = ss; end
[Ydec,F2] = freqz(ydec,1,512,3332);
subplot(4,1,3), plot(F2,abs(Ydec)), legend('Spectrum: decimated signal'),...
    ylabel('|Y_{dec}(F)|')
axis([0,1662,0,100])

STEP 5: Interpolation

E = 3*fliplr(E);

yi = [];

xk = zeros(size(1:length(h)/3));;

rot = 0; % Initial rotator position

k = 1;

for n=1:length(ydec)*3

xn = ydec(k); ek = E(rot+1,:);

if k == 1, xk = [xn,xk(1:length(xk)-1)]; else end

yn = sum(xk.*ek); yi = [yi,yn];

if rot == 0 ss = 1; else end

if rot == 1 ss = 2; else end

if rot == 2 ss = 0; k = k+1; xk = [xn,xk(1:length(xk)-1)]; else end

rot = ss;

end

[Yint,F3] = freqz(yi,1,1024,10000);

subplot(4,1,4)

plot(F3,abs(Yint)); legend('Spectrum: interpolated signal')

xlabel('Frequency, Hz'),ylabel('|Y_{int}(F)|')

axis([0,5000,0,300])

disp('END OF EXAMPLE 4.3')

END OF EXAMPLE 4.3

4.4 Computational Efficiency of Single-Stage and Multistage Systems

Example 4.4

Implementation of 2-stage factor-of-15 decimator where

Sampling rate conversion from 30 kHz to 2 kHz.

Two-stage decimator: (a) Two atage implementation for

(b) Single-stage equivalent for

(developed ftom the 3rd identity).

close all, clear all

Filter

h1 = firpm(18,[0,0.5/15,5/15,7/15,11/15,13/15,1,1],[1,1,0,0,0,0,0,0],[1,5,5,5]);
[H1,F1] = freqz(h1,1,1024,30); % Frequency response  of H1(z)

Filter

h2 = firpm(35,[0,5*0.5/15,5*1/15,1],[1,1,0,0],[1,5]); % Filter H2(z)
[H2,F2] = freqz(h2,1,1024,6); % Frequency response of H2(z)

Periodic filter

h2u = upsample(h2,5); % Periodic filter H_2(z^5)
[H2u,F1] = freqz(h2u,1,1024,30); % Spectrum of the periodic filter H_2(z^5)

The overall frequency response

H = H1.*H2u; % The overall frequency response

figure

subplot(4,1,1), plot(F2,20*log10(abs(H2))), ylabel('|H_2(F)|'), axis([0,3,-80,2])

title('Magnitude responses, dB')

subplot(4,1,2), plot(F1,20*log10(abs(H2u))), ylabel('|H_2_u(F)|'),axis([0,15,-80,2])

subplot(4,1,3), plot(F1,20*log10(abs(H1))), ylabel('|H_1(F)|'), axis([0,15,-80,2])

subplot(4,1,4), plot(F1,20*log10(abs(H))), ylabel('|H(F)|'), axis([0,15,-80,2])

xlabel('Frequency, kHz')

Creating the input signal:

ff = [0,1/15,1.5/15,0.6,0.8,1]; a = [0,0,0.8,0.7,0.9,0];
x = fir2(1000,ff,a);
x = 5*x+0.01*sin(2*pi*(1:1001)*300/30000) + 0.01*sin(2*pi*(1:1001)*400/30000);
[X,F1] = freqz(x,1,2048,30); % Spectrum of the input signal
figure 
subplot(3,1,1), plot(F1,abs(X)), ylabel('|X(F)|'), legend('Original')
title('Signal spectrum: original, 1st stage, 2nd stage')

Polyphase decomposition

h1 = [h1,0];
E1 = reshape(h1,5,length(h1)/5);

Polyphase downsampling and filtering,

v0 = x(1:5:length(x)); v1 = [0,x(5:5:length(x))]; v2 = [0,x(4:5:length(x))]; 
v3 = [0,x(3:5:length(x))]; v4 = [0,x(2:5:length(x))];
y0 = filter(E1(1,:),1,v0); y1 = filter(E1(2,:),1,v1); y2 = filter(E1(3,:),1,v2);
y3 = filter(E1(4,:),1,v3); y4 = filter(E1(5,:),1,v4);
yd1 = y0 + y1 + y2 + y3 + y4; % Decimated signal
[Yd1,F2]=freqz(yd1,1,2048,6);
subplot(3,1,2), plot(F2,abs(Yd1)),
legend('Decimated, 1st stage'), ylabel('|Y_d_1(F)|')

Polyphase decomposition

E2 = reshape(h2,3,length(h2)/3);

Polyphase downsampling and filtering,

x = [yd1,0];

v0 = x(1:3:length(x)); v1 = [0,x(3:3:length(x))]; v2 = [0,x(2:3:length(x))];

y0 = filter(E2(1,:),1,v0); y1 = filter(E2(2,:),1,v1);

y2 = filter(E2(3,:),1,v2);

y = y0 + y1 + y2; % decimated signal

[Y,F3] = freqz(y,1,2048,2);

subplot(3,1,3), plot(F3,abs(Y)), legend('Decimated, 2nd stage'), ylabel('|Y(F)|')

xlabel('Frequency, kHz')

Comment: Efficiency of implementation

Single-stage factor-of-15 FIR decimator needs 180 multiplication constanta
Two-stage 5x3 FIR decimator needs 55 multiplication constants as shown above.

disp('END OF EXAMPLE 4.4')
END OF EXAMPLE 4.4
disp('                    END OF CHAPTER IV')
                    END OF CHAPTER IV