1
$\begingroup$

Since the process can be applied in either domain to increase the sampling rate in the other domain, I am trying to apply zero-padding in frequency space to recover a 'cleaner' interpolated signal in temporal space. To do so, I insert zero-valued frequencies in the spectrum at the location of higher frequencies, which is a common practice.

However I don't seem to recover the original signal very well (in black below) after zero-padding (in red).

import numpy as np
import matplotlib.pyplot as plt

# odd dimension for simplicity
n    = 19
npad = 99

x    = np.linspace(0.,4.*np.pi,n)
xpad = np.linspace(0.,4.*np.pi,npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)//2
f_fwd_pad[0:h+1]   = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[h+1:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

fig, ax = plt.subplots(1,2)

ax[0].plot(x,np.real(f),linestyle=None,marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),linestyle=None,marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

enter image description here

Is that result expected? Is there some fundamental understanding that I am missing?

Improve this question
$\endgroup$
5
  • $\begingroup$ I don't think your bin re-arrangement subscripts are correct. Maybe use h = (n-1)/2 as a substitute to simplify. For starters, since you are using odd (smart) the h point needs to be part of the lower range. $\endgroup$
    – Cedron Dawg
    Jul 17, 2020 at 19:54
  • $\begingroup$ For context, see: dsp.stackexchange.com/questions/59740/… "These formulas are equivalent to taking a DFT, zero padding at the Nyquist or chopping it, then taking the inverse DFT." $\endgroup$
    – Cedron Dawg
    Jul 18, 2020 at 15:18
  • $\begingroup$ You're welcome. MG will have the same issue if he copied the code. The bug is in dropping the bin, not in misinterpreting. Considering the upper half to be positive frequencies is not incorrect, it is just not as "band-limited". See my answer and comments here dsp.stackexchange.com/questions/59316/… then follow the fluffy cloud link for a visual demonstration in the complex case. Note that all the points get hit. The nub of the conversation concerns the impact of splitting the Nyquist as the band limit value. $\endgroup$
    – Cedron Dawg
    Jul 20, 2020 at 11:11
  • $\begingroup$ @CedronDawg mmm unless I misunderstood your comment I can confirm that the problem was considering a time domain with an overlapping discrete point in the period. Running the code of the OP with the corrected bin slices still yields the original problem. $\endgroup$
    – circuitbreaker
    Jul 20, 2020 at 11:21
  • $\begingroup$ You may have a look at dsp.stackexchange.com/questions/79356 for analytic derivation. $\endgroup$
    – Eric Johnson
    Jan 22, 2022 at 13:34

3 Answers 3

Reset to default
4
$\begingroup$

The answer above is correct. Just to clarify a bit further, using x = np.linspace(0,10,5) will produce 5 numbers from 0 to 10 inclusively

np.linspace(0,10,5)
array([ 0. ,  2.5,  5. ,  7.5, 10. ])

You don't want the last number because in your example the last number is the first number of the next period. A correct implementation would be:

periods = 4.*np.pi
x    = np.arange(0., periods, periods/n)
xpad = np.arange(0., periods, periods/npad)

Result with appropriate period

Improve this answer
$\endgroup$
6
  • $\begingroup$ Please double check this: $$ \text{h=(n-1) / 2} $$ $$ \text{f_fwd_pad[0:h+1] = f_fwd[0:h+1]} $$ $$ \text{f_fwd_pad[npad-h:] = f_fwd[h+1:]} $$ Sort of the same issue. I'll be deleting this comment shortly. $\endgroup$
    – Cedron Dawg
    Jul 18, 2020 at 2:43
  • $\begingroup$ Ok, so the period needs to be exactly one period (or a multiple) without overlapping. And yes Cedron there was a typo in my indexes I forgot to change some of them for the 0-based indexing in python $\endgroup$
    – circuitbreaker
    Jul 19, 2020 at 9:43
  • $\begingroup$ @CedronDawg I corrected the typos in the OP thanks for spotting that $\endgroup$
    – circuitbreaker
    Jul 19, 2020 at 10:31
  • $\begingroup$ @circuitbreaker I'm replying here since your responded here and MG will see it too. You still missed a "+1", I edited it in for you. Python ranges don't include the last point. You want your bin set to include bin "h". Then the upper range starts at "h+1" as the syntax shows. $\endgroup$
    – Cedron Dawg
    Jul 19, 2020 at 13:39
  • $\begingroup$ Thanks you @CedronDawg that was another mistake indeed. So first range is DC component + h positive frequencies so 0:h+1 in python And second range is h negative frequencies so h+1: $\endgroup$
    – circuitbreaker
    Jul 20, 2020 at 9:34
1
$\begingroup$

Your problem is here:

x    = np.linspace(0.,4.*np.pi,n)

The period should be N+1 samples from 0 to 2pi, for n=0..N. Then take x(k) for k=0..N-1

Currently your FFT is not a pure single tone, because the sinusoid does not have a perfect period within the FFT period. And so padding with zeros would not be the correct padding. The fix above will a perfect period within the FFT period, to make padding with zero the correct padding.

Improve this answer
$\endgroup$
0
$\begingroup$

Here is your code with MG's fix applied. A few other tweaks, and a fluffy cloud. The pics look fine to me. Command line Python 2.7

import numpy as np
import matplotlib.pyplot as plt

# odd dimension for simplicity
n    = 19
npad = 299

duration = 0.85*2.0*np.pi
x    = np.arange(0., duration, duration/n)
xpad = np.arange(0., duration, duration/npad)

f = np.cos(x) + 1j*np.sin(x)

f_fwd = np.fft.fft(f)

f_fwd_pad = np.zeros(npad,dtype=complex)

h = (n-1)/2
h = n - 3
f_fwd_pad[0:h+1]   = f_fwd[0:h+1]
f_fwd_pad[npad-h:] = f_fwd[n-h:]

f_interpolated = np.fft.ifft(f_fwd_pad)*npad/n

fig, ax = plt.subplots(1,2)

ax[0].plot(x,np.real(f),marker='x',color='k')
ax[0].plot(xpad,np.real(f_interpolated),color='r')

ax[1].plot(x,np.imag(f),marker='x',color='k')
ax[1].plot(xpad,np.imag(f_interpolated),color='r')

plt.show()

plt.plot( np.real(f_interpolated), np.imag(f_interpolated) )
plt.show()

enter image description here

Improve this answer
$\endgroup$

Your Answer

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

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