Code

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