3
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ There should be a MATLAB SE somewhere. Maybe SO. I've used resample before, but this forces me to review it. $\endgroup$ Commented Feb 23, 2023 at 3:57

2 Answers 2

4
$\begingroup$

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);
$\endgroup$
8
  • $\begingroup$ (Nice answer - I would think it would hurt to multiply by a hamming window since least squares will give you the least square solution directly?) $\endgroup$ Commented Feb 23, 2023 at 4:47
  • $\begingroup$ While I agree it hurts the transition band, it does flatten out the pass-band and lowers the stop-band. A different window probably gives a better trade-off than Hamming though... I'll add that in. $\endgroup$
    – Jdip
    Commented Feb 23, 2023 at 5:02
  • $\begingroup$ I don't think so-- if you relax your transition band within the least squares target to match the same transition band you degraded to, you should find that the least squares in that case wins (since it is the optimized solution in the least squares sense for that given transition band). Kaiser comes close to the least squares result but doesn't beat it. So you see the improvement which (I think) is fooling you---I believe you'll get that improvement and more if you relax your targets to that degraded transition. $\endgroup$ Commented Feb 23, 2023 at 5:14
  • 1
    $\begingroup$ @DanBoschen Matlab's resample applies a Kaiser window on a least square design with adjacent band edges. I think that's their way to flatten the pass band and lower the stop band without having to specify a transition bandwidth. I see your point and I agree, the latter method would, in a least-square sense, give better performance. I've edited! $\endgroup$
    – Jdip
    Commented Feb 23, 2023 at 16:23
  • 1
    $\begingroup$ @DanBoschen here is the relevant snippet: h = firls( L-1, [0 2*fc 2*fc 1], [1 1 0 0]).*kaiser(L,bta)' ; $\endgroup$
    – Jdip
    Commented Feb 23, 2023 at 17:25
0
$\begingroup$

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

enter image description here

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])

enter image description here

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

enter image description here

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

enter image description here

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'

enter image description here

enter image description here

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

enter image description here

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

% 007

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

$\endgroup$
1
  • $\begingroup$ This is very helpful, thank you $\endgroup$
    – Levi
    Commented Mar 3, 2023 at 4:03

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.