0
$\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?

$\endgroup$
  • $\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 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 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 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 at 11:21
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$
  • $\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 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 at 9:43
  • $\begingroup$ @CedronDawg I corrected the typos in the OP thanks for spotting that $\endgroup$ – circuitbreaker Jul 19 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 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 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, privacy policy and cookie policy

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