Matlab "resample" function

Question

Say I have a signal x(n) and I want to resample from rate 30Hz to 20Hz. In Matlab, I can make a filter with:

fs1 = 30;  
fs2 = 20;

f = [0 6 7 fs2/2]/(fs2/2);  
a = [1 1 0 0];  
n = 31;  
b = firpm(n, f, a);

This way my cutoff frequency is 7Hz. Now I want to use this as my interpolation filter. I can run

x2 = resample(x, fs2, fs1, b);

The filter b from firpm I understand is normalized to rate fs2. But, is this what the resample function expects, or does it want a filter normalized to fs1? In other words, should I have run

f = [0 6 7 fs2/2]/(fs2/2);  

or

f = [0 6 7 fs1/2]/(fs1/2);

The documentation does not make this clear, nor could I find anything elsewhere online.

Answer 1

First, a couple definitions I use:

• Interpolation is the process by which a signal is up-sampled (by inserting 0s, or "zero-stuffing"), then low-pass filtered.
• Decimation is the process by which a signal is low-pass filtered, then down-sampled (throwing away samples).

---

There are multiple ways the resampling can be done since the ratio is non-integer. A straight-forward path is:

1. interpolate-by-2 to fs3 = 60Hz
2. decimate-by-3 to fs2 = 20Hz

1. As you are aware, when interpolating, you want to filter after the zero-insert operation, since it creates images at multiples of the original sampling frequency. In our case, the filter should have a cut-off at $f_{s_1}/2 = 15 \,\texttt{Hz}$
2. When decimating, we need to filter before the down-sampling operation, since it introduces aliasing of the components above the resulting sampling rate Nyquist frequency. In our case, this filter should have a cut-off at $f_{s_2}/2 = 10 \texttt{Hz}$

Notice, both filters run on the data at $f_{s_3}= 60 \texttt{Hz}$

We can hence combine these two filters into one, specifically into the one with lower cut-off, i.e. our decimation filter with cut-off $f_{s_2}/2 = 10\texttt{Hz}$. The filtering will be done on the up-sampled data (not the interpolated data, since we haven't low-pass filtered it yet) with rate $f_{s_3}= 60\texttt{Hz}$, so if you really insist on specifying the filter yourself, you would want to design it with $f_{s_3}$ in mind.
Also, a least-square design is probably better suited for this task:

fs1 = 30;
fs2 = 20;
fs3 = 60;

f = [0 fs2/2 fs2/2 fs3/2] / (fs3/2);
a = [1 1 0 0];
n = 31; 
upsampFac = fs3/fs1; 
b = firls(n-1, f, a) .* hamming(n)'; % Matlab's resample uses a Kaiser window.

% % Dan Boschen's suggestion: specify the bandwidth yourself!
% B = 6; % transition bandwidth, in Hz. This transition bandwidth is very 
% %close to the one obtaine through the above windowing method, but gives better performance.
% f = [0 fs2/2-B/2 fs2/2+B/2 fs3/2] / (fs3/2); 
% b = firls(n-1, f, a); % No need for windowing!

b = upsampFac * b / sum(b);

Full disclosure, this is in essence what Matlab's resample function does if you call it without the extra argument:

x2 = resample(x, fs2, fs1);

Answer 2

1.- Signal Integrity :

2 samples/cycle on band edge is not enough

Let's see what happens with 2 samples/cycle on band edge :

close all;clear all;clc

f01 = 30;  
f02 = 20;

t = [0:1/f01:1-1/f01];
f=100./(1:5);

x = sum(cos(2*pi*t.*f'),1);

figure(1)
plot(t,x)
grid on
title(['fs = 2 x ' num2str(f01)])
xlabel('t');ylabel('x');

I use f01 f02 for 30kHz 20kHz instead of fs1 fs2 because 30kHz 20kHz are not the sampling frequencies.

fs=10*f01;
t = [0:1/fs1:1-1/fs1];

x = sum(cos(2*pi*t.*f'),1);

figure(2)
plot(t,x)
grid on
title(['fs = 10 x ' num2str(f01)])
xlabel('t');ylabel('x');
axis([0 .2 -5 5])

not enough yet

fs1=25*f01;
t = [0:1/fs1:1-1/fs1];

x = sum(cos(2*pi*t.*f'),1);

figure(3)
plot(t,x)
grid on
title(['fs = 25 x ' num2str(f01)])
xlabel('t');ylabel('x');
axis([0 .2 -5 5]);

Now there's enough signal integrity .

2.- Filter cut-off setup :

2.1.- From the question; with f = [0 6 7 fs2/2]/(fs2/2);a=[1 1 0 0]; this means f*fNyq=[0 6 7 10] and the cut-off frequency is 6Hz, not the 7Hz mentioned in the question.

2.2.- In Jdip's answer fs2=20;fs3=60;f=[0 fs2/2 fs2/2 fs3/2]/(fs3/2)];

This means f=[0 1/3 1/3 1] and for a=[1 1 0 0] there's no need to repeat the same frequency, and now f=[0 10 10 30] and the cut-off takes place at 10Hz, not the sought 7Hz.

2.3.- The suggestion closing the question

f = [0 6 7 fs2/2]/(fs2/2) this is f=[0 .6 .7 1] (with a=[1 1 0 0]) cut-off still at 6Hz.

f = [0 6 7 fs1/2]/(fs1/2) this means f=[0 6 7 15] cut-off still at 6Hz .

3.- kHz, not Hz, and more points in fn

Let's increase the amount frequencies to build the filter with, the amount of points in fn because 7Hz cut-off is a really low cut-off, given 30k and 20kHz.

There's need for higher frequency resolution indeed, moreover aiming at achieving this with a 31 coefficient FIR.

Setting fs1=30e3;fs2=20e3; instead of 30Hz 20Hz used in the question code.

The input signal is

f01=30e3;  % [Hz]
f02=20e3;
fs1=25*f01;
t = [0:1/fs1:30-1/fs1];
x = sum(cos(2*pi*t.*f'),1);

figure
plot(t,x)
grid on; xlabel('t');ylabel('x'); 

And now something curious happens

fc0=7
N1=10   % points on pass band
N2=20000  % points on stop band
f=[linspace(0,fc0,N1) linspace(fc0+100, fs1/2,N2)]/(fs1/2);
a=[ones(1,N1) zeros(1,N2)];

n = 31;  % filter order

b = firpm(n, f, a);

x2 = resample(x, f02, f01, b);
t2=[0:1:numel(x2)-1]*1/f02;
%[1:1/fs2:10-1/fs1];

figure
plot(t2,x2)
grid on; xlabel('t');ylabel('x2'); 

[h,w] = freqz(b,1,512);
figure
plot(f,a,w/pi,abs(h))
grid on
title(['N2 = ' num2str(N2)])
legend('Ideal','firpm Design')
xlabel 'Radian Frequency (\omega/\pi)', ylabel 'Magnitude'

Warning: Failed to converge after 3 iterations.  This is
likely due to machine rounding error.  If the number of
iterations exceeds 3, the design may be correct, but should
be verified with an FFT.

This warning is kind of : .. reaching the calculations limit of command resample.

And not even meeting the expected zeros on fn values along stop band.

4.- Stop band level is too high

I have repeated the previous cell with N2=[200 2000 20000] obtaining same graph

for all N2 values.

This frequency response is not satisfactory towards obtaining a resampled signal.

5.- More efficient alternative way

• Either we have reached the operational limits of command resample or too much time is needed to configure it correctly.

A working riffle has to be easy and quick to load, both at the same time, for a user with an average IQ.

• The resampling factor is p/q=f02/f01=2/3 meaning that going down in frequency the current sampling frequency fs is good enough for x2 as well.

• Good signal integrity : f_Nyquist >> max(fin)/2 .

Therefore, slowing down by 2/3 means increasing the amount of samples per cycle by 3/2.

f01=30e3;  % [Hz]  q
f02=20e3;  % [Hz]  p
fs1=25*f01;  % [Hz] sampling frequency

t = [0:1/fs1:30-1/fs1];    % [s]
x = sum(cos(2*pi*t.*fn'),1);   % input signal

T01=1/f01
T02=1/f02

T01*fs1  % samples per cycle
T02*fs1

figure
plot(t,x,'LineWidth',2)
hold on;grid on; xlabel('t');
plot(t,x2,'LineWidth',2)
axis([0 .2 -5 5])
legend({'xin' 'xout'})

Objective reached without using resample or any demanding filter that it is difficult or impossible to implement with firpm .