Resampling filter performance

kadhiem Ayob January 14, 20123 comments Coded in Matlab

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.

The code uses test parameters but you can choose your own input, your filter, your sampling rate, your required interpolation/decimation rates and your input signal centre frequency.

%%%%%%%%%%%%%%%%%% 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');

Comments:

kaz
Said:
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).
4 years ago
0
Reply
Sorry, you need javascript enabled to post any comments.
Rick Lyons
Said:
Hello Kadhiem Ayob, 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-]
4 years ago
0
Reply
Sorry, you need javascript enabled to post any comments.
kaz
Said:
Hi Rick, Seems like copy paste issue of those characters like 'r.-' You may need to edit them or remove the dot. Regards kadhiem
4 years ago
0
Reply
Sorry, you need javascript enabled to post any comments.
Sorry, you need javascript enabled to post any comments.