Search tips

# Resampling filter performance

Language: Matlab

Processor: Not Relevant

Submitted by kadhiem Ayob on Jan 14 2012

# Resampling filter performance

A popular method for upsampling is zero insertion. A regular zero interleaving of a signal results in an extra copy in its frequency domain. If this copy is removed by a filter that keeps original signal then the filter in effect creates interpolation points. If the copy is not removed efficiently then this implies less efficient interpolation in time domain.

Decimation can be done on a signal by discarding samples provided aliasing is avoided by removing frequencies above half of final sampling frequency.

If both upsampling and downsampling are planned then a single filter after zero insertion and before decimation may be enough.

The following code can check filter performance for interpolation and decimation. Additionally it checks the effect of moving signal band or filter envelope.

%%%%%%%%%%%%%%%%%% inputs for model %%%%%%%%%%%%%%%%
clear all; close all;

%example bandlimited impulse input & parameters
x = filter(fir1(70,.1),1,[1 zeros(1,2^15-1)]);
Fs = 120;           %MHz original sample rate
h = fir1(30,.3);    %filter used for upsampling
up = 3;             %Interpolation factor
dn = 2;             %Decimation factor
Fc = 12;             %MHz band centre (-Fs/2 ~ +Fs/2)
Fch = 0;             %MHz filter centre (-Fs*up/2 ~ +Fs*up/2)

%move signal to its centre
x = x.*exp(j*2*pi*(0:length(x)-1)*Fc/Fs);

%shift filter
h = h.*exp(j*2*pi*(0:length(h)-1)*Fch/(Fs*up));

%%%%%%%%%%%%%%%%%%%%% model %%%%%%%%%%%%%%%%%%%%%%
%check signal in upsampled domain
x_up = zeros(1,length(x)*up);
x_up(1:up:end) = x;
[P, F] = pwelch(x_up, [], 0, 2^16, Fs*up,'twosided');
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100;  %avoid log of zero
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P));  %normalise

%check filter response in upsampled domain
H = fftshift(20*log10(abs(fft(h,2^16))));

subplot(2,1,1);
hold;grid;
plot(F, P_dB,'.-');
plot(F,H,'m--');
axis([min(F)-1 max(F)+1 -80 1]);
legend('upsampled signal','upsampling filter');

%check signal in downsampled domain
x_f = filter(h,1,x_up);
x_dn = x_f(1:dn:end);
[P, F] = pwelch(x_dn, [], 0, 2^16, Fs*up/dn,'twosided');
F = F - max(F)/2;
P = fftshift(P);
y(find(P == 0)) = -100;   %avoid log of zero
y(find(P ~= 0)) = 10*log10(P(find(P ~= 0)));
P_dB = y - 10*log10(max(P)); %normalise

subplot(2,1,2)
plot(F,P_dB,'r.-');
grid;
axis([min(F)-1 max(F)+1 -80 1]);
legend('downsampled signal');

Rate this code snippet:
0

Experienced FPGA Engineer, focussed on DSP functionality within FPGAs

# kaz wrote:

1/17/2012

Just to add a bit more. As you can see the interpolation does not like that tidy in the stop band. This is true for many other practical hardware methods. An NCO can be thought of as an ideal interpolator; if, for example, you double the sampling frequency and double the points you choose from its LUT then you get same output frequency interpolated ideally (precomputed interpolation).

# Rick Lyons wrote:

2/12/2012

I ran your code.  In the top panel of Figure 1 there's supposed to be a solid blue curve and a dashed magenta curve.  But only the magenta curve is displayed in the top panel.  In the bottom panel of Figure 1 there's supposed to be a solid red curve, but the bottom panel contains no curve at all.  Am I doing something wrong?

[-Rick-]

# kaz wrote:

2/12/2012

Hi Rick,

Seems like copy paste issue of those characters like 'r.-'
You may need to edit them or remove the dot.

Regards