What is the correct method to upsample?

Question

I have an array of samples at 75 Hz, and I want to store them at 128 Hz. If it was 64 Hz and 128 Hz it was very simple, I would just double all samples. But what is the correct way if the samplerates are not a fraction of eachother?

Answer 1

When you want to avoid Filtering then you can:

1. handle signal as set of joined interpolation cubics curves

   but this point is the same as if you use linear interpolation. Without knowing something more about your signal and purpose you can not construct valid coefficients (without damaging signal accuracy) for example of how to construct such cubic look here:

   • my interpolation cubic

   in bullet #3 inside that link are coefficients I use. I think there are sufficient even for your purpose so you can try them. If you want to do custom interpolation look here:

   • how to construct custom interpolation curve

2. create function that can return point in your signal given time from start of sampling

   so do something like

   double signal(double time);
   

   where time is time in [s] from start of sampling. Inside this function compute which 4 samples you need to access.

   ix = floor(time*75.0);
   

   gives you curve start point sample index. Cubic need 4 points one before curve and one after ... so for interpolation cubic points p0,p1,p2,p3 use samples ix-1,ix,ix+1,ix+2. Compute cubic coefficients a0,a1,a2,a3 and compute cubic curve parameter t (I use range <0,1>) so

   t=(time*75.0); t-=floor(t);
   

   

   • green - actual curve segment
   • aqua - actual curve segment control points = 75.0 Hz samples
   • red - curve parametric interpolation parameter t
   • gray - actual time

   sorry I forgot to draw the actual output signal point it should be the intersection of green and gray

3. simply do for loop through sampled data with time step 1/128 s

   something like this:

   double time,duration=samples*75.0,dt=1.0/128.0;
   double signal128[???];
   for (time=0.0,i=0;time<duration;i++,time+=dt)
    signal128[i]=signal(time);
   

   samples are the input signal array size in samples sampled by 75.0 Hz

[notes]

• for/duration can be done on integers ...
• change signal data type to what you need
• inside signal(time) you need to handle edge cases (start and end of signal)
• because you have no defined points in signal before first sample and after last sample. You can duplicate them or mirror next point (mirror is better).
• this whole thing can be changed to process continuously without buffers just need to remember 4 last points in signal so you can do this in RT. Of coarse you will be delayed by 2-3 75.0 Hz samples ... and when you put all this together you will see that this is a FIR filter anyway :)
• if you need to preserve more then first derivation add more points ...

Answer 2

You do not need to upsample and then downsample.

Instead, one can interpolate all the new sample points, at the desired spacing in time, using a wide enough low-pass interpolation kernel, such as a windowed Sinc function. This is often done by using a pre-calculated polyphase filter bank, either directly, or with an addition linear interpolation of the filter table. But if performance is not critical, then one can directly calculate each coefficient for each interpolated point.

Answer 3

The easiest way is to upsample to a sample rate which is the LCM of your two sample rates and then downsample - that way you get integer upsample/downsample ratios. In your case there are no common factors in the two sample rates so you would need to upsample by a factor of 128 to 9.6 kHz and then downsample by a factor of 75 to 128 Hz. For the upsampling you insert 127 0 samples in between each sample, then apply a suitable filter (37 Hz LPF, Fs = 9.6 kHz), and then downsample by taking every 75th sample. The filter design is the only tricky part, but there are online tools for taking the hard work out of this.

Alternatively look at third-party libraries which handle resampling, e.g. sox.

Answer 4

You need to upsample and downsample with an intermediate sampling frequency, as @Paul mentioned. In addition, it is needed to filter the signal after each transformation, which can be achieved by linear interpolation as:

% Parameters
F = 2;
Fs1 = 75;
Fs3 = 128;
Fs2 = lcm(Fs1,Fs3);

% Original signal
t1 = 0:1/Fs1:1;
y1 = sin(2*pi*F*t1);

% Up-sampled signal
t2 = 0:1/Fs2:1;
y2 = interp1(t1,y1,t2);

% Down-sampled signal
t3 = 0:1/Fs3:1;
y3 = interp1(t2,y2,t3);

figure;
subplot(3,1,1);
plot(t1,y1,'b*-');
title(['Signal with sampling frequency of ', num2str(Fs1), 'Hz']);
subplot(3,1,2);
plot(t2,y2,'b*-');
title(['Signal with sampling frequency of ', num2str(Fs2), 'Hz']);
subplot(3,1,3);
plot(t3,y3,'b*-');
title(['Signal with sampling frequency of ', num2str(Fs3), 'Hz']);