Subsampling an array of numbers

Question

I have a series of 100 integer values which I need to reduce/subsample to 77 values for the purpose of fitting into a predefined space on screen. This gives a fraction of 77/100 values-per-pixel - not very neat.

Assuming the 77 is fixed and cannot be changed, what are some typical techniques for subsampling 100 numbers down to 77. I get a sense that it will be a jagged mapping, by which I mean the first new value is the average of [0, 1] then the next value is [3], then average [4, 5] etc. But how do I approach getting the pattern for this mapping?

I am working in C++, although I'm more interested in the technique than implementation.

Thanks in advance.

Answer 1

Either if you downsample or you oversample, you are trying to reconstruct a signal over nonsampled points in time... so you have to make some assumptions.

The sampling theorem tells you that if you sample a signal knowing that it has no frequency components over half the sampling frequency, you can continously and completely recover the signal over the whole timing period. There's a way to reconstruct the signal using sinc() functions (this is sin(x)/x)

sinc() (indeed sin(M_PI/Sampling_period*x)/M_PI/x) is a function that has the following properties:

1. Its value is 1 for x == 0.0 and 0 for x == k*Sampling_period with k == 0, +-1, +-2, ...
2. It has no frequency component over half of the sampling_frequency derived from Sampling_period.

So if you consider the sum of the functions F_x(x) = Y[k]*sinc(x/Sampling_period - k) to be the sinc function that equals the sampling value at position k and 0 at other sampling value and sum over all k in your sample, you'll get the best continous function that has the properties of not having components on frequencies over half the sampling frequency and have the same values as your samples set.

Said this, you can resample this function at whatever position you like, getting the best way to resample your data.

This is by far, a complicated way of resampling data, (it has also the problem of not being causal, so it cannot be implemented in real time) and you have several methods used in the past to simplify the interpolation. you have to constructo all the sinc functions for each sample point and add them together. Then you have to resample the resultant function to the new sampling points and give that as a result.

Next is an example of the interpolation method just described. It accepts some input data (in_sz samples) and output interpolated data with the method described before (I supposed the extremums coincide, which makes N+1 samples equal N+1 samples, and this makes the somewhat intrincate calculations of (in_sz - 1)/(out_sz - 1) in the code (change to in_sz/out_sz if you want to make plain N samples -> M samples conversion:

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

/* normalized sinc function */
double sinc(double x)
{
    x *= M_PI;
    if (x == 0.0) return 1.0;
    return sin(x)/x;
} /* sinc */

/* interpolate a function made of in samples at point x */
double sinc_approx(double in[], size_t in_sz, double x)
{
    int i;
    double res = 0.0;
    for (i = 0; i < in_sz; i++)
            res += in[i] * sinc(x - i);
    return res;
} /* sinc_approx */

/* do the actual resampling.  Change (in_sz - 1)/(out_sz - 1) if you
 * don't want the initial and final samples coincide, as is done here.
 */
void resample_sinc(
    double in[],
    size_t in_sz,
    double out[],
    size_t out_sz)
{
    int i;
    double dx = (double) (in_sz-1) / (out_sz-1);
    for (i = 0; i < out_sz; i++)
            out[i] = sinc_approx(in, in_sz, i*dx);
}

/* test case */
int main()
{
    double in[] = {
            0.0, 1.0, 0.5, 0.2, 0.1, 0.0,
    };

    const size_t in_sz = sizeof in / sizeof in[0];
    const size_t out_sz = 5;
    double out[out_sz];
    int i;

    for (i = 0; i < in_sz; i++)
            printf("in[%d] = %.6f\n", i, in[i]);
    resample_sinc(in, in_sz, out, out_sz);
    for (i = 0; i < out_sz; i++)
            printf("out[%.6f] = %.6f\n", (double) i * (in_sz-1)/(out_sz-1), out[i]);

    return EXIT_SUCCESS;
} /* main */

Answer 2

Everything depends on what you wish to do with the data - how do you want to visualize it.

A very simple approach would be to render to a 100-wide image, and then smooth scale the image down to a narrower size. Whatever graphics/development framework you're using will surely support such an operation.

Say, though, that your goal might be to retain certain qualities of the data, such as minima and maxima. In such a case, for each bin, you're drawing a line of darker color up to the minimum value, and then continue with a lighter color up to the maximum. Or, you could, instead of just putting a pixel at the average value, you draw a line from the minimum to the maximum.

Finally, you might wish to render as if you had 77 values only - then the goal is to somehow transform the 100 values down to 77. This will imply some kind of an interpolation. Linear or quadratic interpolation is easy, but adds distortions to the signal. Ideally, you'd probably want to throw a sinc interpolator at the problem. A good list of them can be found here. For theoretical background, look here.

Answer 3

There are different ways of interpolation (see wikipedia)

The linear one would be something like:

std::array<int, 77> sampling(const std::array<int, 100>& a)
{
     std::array<int, 77> res;

     for (int i = 0; i != 76; ++i) {
         int index = i * 99 / 76;
         int p = i * 99 % 76;

         res[i] = ((p * a[index + 1]) + ((76 - p) * a[index])) / 76;
    }
    res[76] = a[99]; // done outside of loop to avoid out of bound access (0 * a[100])
    return res;
}

Live example

Answer 4

Create 77 new pixels based on the weighted average of their positions.

As a toy example, think about the 3 pixel case which you want to subsample to 2.

Original (denote as multidimensional array original with RGB as [0, 1, 2]):

|----|----|----|

Subsample (denote as multidimensional array subsample with RGB as [0, 1, 2]):

|------|------|

Here, it is intuitive to see that the first subsample seems like 2/3 of the first original pixel and 1/3 of the next.

For the first subsample pixel, subsample[0], you make it the RGB average of the m original pixels that overlap, in this case original[0] and original[1]. But we do so in weighted fashion.

subsample[0][0] = original[0][0] * 2/3 + original[1][0] * 1/3  # for red
subsample[0][1] = original[0][1] * 2/3 + original[1][1] * 1/3  # for green
subsample[0][2] = original[0][2] * 2/3 + original[1][2] * 1/3  # for blue

In this example original[1][2] is the green component of the second original pixel.

Keep in mind for different subsampling you'll have to determine the set of original cells that contribute to the subsample, and then normalize to find the relative weights of each.

There are much more complex graphics techniques, but this one is simple and works.