Search code examples
cmathsignal-processingsamplingresampling

What is the correct method to upsample?


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?


Solution

  • 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:

      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:

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

      enter image description here

      • 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 ...