1
$\begingroup$

Upsampling a spike [100 0 0 ...] with rfft - irfft gives quite different results for input lengths even or odd. Here's a plot for 100 and 101: 30oct-spike-upsample.png

Why the difference ?
An intuitive explanation which could predict this on paper, without running it, would be nice.


TL;DR

The curve for odd $\to $ 2*odd goes $\sim 0, 1/n, 0, -1/n, 0 ...$ in the middle; that for even $\to$ 2*even approaches 0 in the middle (quadratically ?)

See scipy.signal.resample . The code in resample .py appears to do:

m2 = n // 2 + 1  # number of relevant frequency bins of a one-sided FFT
X = rfft(x) [:m2]
if n is even:  # Account for unpaired bin at m//2:
    X[n//2] *= 0.5
0-pad  # inside irfft
irfft

#!/usr/bin/env python3
"""Upsample a spike [100 0 0 ...] with periods even / odd
rfft, 0-pad, irfft (scipy.signal.resample)
"""
# https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.resample.html
# https://github.com/scipy/scipy/blob/main/scipy/signal/_signaltools.py

import numpy as np
from scipy.signal import resample
import matplotlib.pyplot as pl
import sys

print( 80 * "▄" )
np.set_printoptions( edgeitems=10, threshold=10, linewidth=120, formatter = dict(
        float = lambda x: "0" if abs(x) < 1e-10  else "%.2g" % x ))  # 2 digits, tiny -> 0

def spike( n ):
    x = np.zeros( n )
    x[0] = 100
    return x

#...............................................................................
periods = [100, 101]
zoom = 2
plot = 0
    # to change these params, run this.py  a=1  b=2,None  'c = expr' ... in sh or ipython --
for arg in sys.argv[1:]:
    exec( arg )
print( f"{periods = }  {zoom = } " )

title = __doc__
print( title, "\n" )
if plot:
    fig, axes = pl.subplots( nrows=len(periods), figsize=[8, 5] )
    fig.suptitle( title, multialignment="left" )
else:
    axes = len(periods) * [None]

for period, ax in zip( periods, axes ):
    x = spike( period )
    nout = zoom * period
#...............................................................................
    y = resample( x, nout )
    jmid = nout // 2
    ymid = y[jmid-5 : jmid+5+1]
    print( "period %d  y middle * 100: %s " % ( period, ymid * 100 ))
    if plot:
        title = r"$ %d \to %d $" % ( period, nout )
        ax.set_title( title, loc="center", y=0.9, fontsize=12 )
        ax.set( ylim=[-2, 2] )
        ax.plot( y, "s-", lw=0.36, ms=1, color="crimson" )
"""
period 100  y middle * 100: [-7.9  0 4.7   0 -1.6   0  -1.6   0 4.7  0 -7.9]
period 101  y middle * 100:    [0 99   0 -99    0  99     0 -99   0 99    0]
"""
$\endgroup$
4
  • $\begingroup$ You can please post the complete code required to reproduce your graphs? Resampling is complicated and even small details matter a lot. Chances are the main difference is due how the Nyquist frequency is handled $\endgroup$ Commented Oct 31 at 14:31
  • $\begingroup$ @Hilmar, done, but it just calls scipy resample. (python3 this.py periods=500,501 runs with other periods.) The lines after "The code in resample .py appears to do ..." in the question are my take on some of that code. Can anyone reproduce this in matlab ? $\endgroup$ Commented Oct 31 at 16:41
  • $\begingroup$ @Hilmar, do you think this is a bug in scipy resample, or just a silly corner case ? $\endgroup$ Commented Nov 5 at 16:16
  • 1
    $\begingroup$ No. It's likely a fundamental property. Whenever you non-zero energy at the Nyquist frequency, you have aliasing that aliasing is affected by the resampling in different ways $\endgroup$ Commented Nov 6 at 18:22

1 Answer 1

0
$\begingroup$

The ripples you get from this resampling (equivalent to a sinc interpolation kernel) taper off quite slowly. The interpolation ripples

FFTs are circular operations (or: they assume their inputs/outputs are periodic). So what you're seeing here is the ripples from one direction overlayed onto the other - and since when doing a 2x oversample, those ripples alternate ±, it will make quite a difference whether those two are phase-aligned or not.

If you use a much larger size, you'll probably see how slow the ripples taper off, and get an intuitive look at how much they overlap at shorter lengths.

$\endgroup$

Your Answer

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

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.