0
$\begingroup$

I am new here and also am not very knowledgeable about DSP so this might be dumb and easy. I am aware of the fact that in order to reconstruct a signal, I need to sample it at a frequency that is more than twice its highest frquency.

In the example below, I use N=50 samples and expect the reconstruction to fail if I use freq = 25 or so (perhaps due to numerical round-off). But I get that the reconstruction breaks around freq = 12.5. I feel I am missing something fundamental.

Here are plots and the code to generate them:

Success

Failure

from numpy.fft import rfft as rfft
from numpy.fft import irfft as irfft
import numpy as np
import matplotlib.pyplot as plt
import math

def interpolant( f_hat, pts ):
    '''
    This is the trigonometric polynomial interpolating
    the signal that is encoded in the values of f_hat.
    it is evaluated at the points pts

    The commented part is the true calculation to be done,
    the uncommented just does it more accurately.
    '''
    N = len( f_hat )
    f = f_hat[0]
    for k in range( 1, N/2 ):
        f = f + 2 * f_hat[k].real * np.cos( k * pts )
        f = f - 2 * f_hat[k].imag * np.sin( k * pts )
        #f = f + f_hat[k]             * np.exp( 1j * k * pts )
        #f = f + f_hat[k].conjugate() * np.exp(-1j * k * pts )
    f = f + f_hat[N/2] * np.cos( N/2 * pts )
    return f / N / 2

def f( pts, freq ):
    '''
    this is "the conitinuous" signal, evaluated (sampled) 
    at pts
    '''
    return  np.sin( freq * pts ) 

# Number of sampled points
N = 50

# Where we sample the signal
pts = np.linspace(0, 2 * np.pi, num = N, endpoint = False )

# A much finer grid, used solely to display results
oversampled = np.linspace( 0, 2 * np.pi, 20 * N, endpoint = False )

# The grid we use for the plotting
grid = oversampled

# The frequncy of the sine wave below
freq = 7

# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )

plt.plot( grid, interpolant( f_hat, grid )  , color = "g" )
plt.plot( grid, f( grid, freq )      , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green is interpolant. Reconstruction succeeds"
plt.title( title1 + title2 )
plt.show()

# The frequncy of the sine wave below
freq = 12.9

# Do the FFT on the sampled signal
f_hat = rfft( f( pts, freq ) )

plt.plot( grid, interpolant( f_hat, grid )  , color = "g" )
plt.plot( grid, f( grid, freq )      , color = "r" )
title1 = str(N) + " samples. Signal frequency is " + str(freq)+" \n"
title2 = "Red is true, green interpolant. Reconstruction fails"
plt.title( title1 + title2 )
plt.show()
Improve this question
$\endgroup$

2 Answers 2

Reset to default
2
$\begingroup$

While hotpaw2's answer is correct in general, in your case it is just a bug in your code. The N in your interpolant function is just 25 as far as I can tell, and then you divide it by 2 again.

Improve this answer
$\endgroup$
5
  • $\begingroup$ Can you be more specific? I tested @hotpaw2's suggestion and it actually works. Do you mean the division by 2 in the interpolant funtion? That seems to be necessary, though I've yet figured out why. $\endgroup$
    – Yair Daon
    Sep 1, 2015 at 19:10
  • $\begingroup$ I think it should be N-1 instead of N/2. Otherwise, you are omitting half of the array. $\endgroup$
    – Sebastian Reichelt
    Sep 1, 2015 at 19:31
  • $\begingroup$ You mean the iteration? I'm pretty confident that's OK - I am using a real DFT and for real DFT the coefficients satisfy $Y_k = \bar{Y_{N-k}}$ the numpy implementation outputs onlly $N/2$ coefficients. See the two commented lines in the loop for the actual calculation that is being done - the first corresponds to the $k$ summand whereas the second corresponds to the $N-k$ summand. My comment above referred to the final division by $2$ that seems to be necessary. $\endgroup$
    – Yair Daon
    Sep 1, 2015 at 20:41
  • $\begingroup$ Yes, but you have two variables named "N", one is 50 and the other is just 25. Then you divide by 2 and get 12. I just checked it; N-1 is correct and produces the graph you would expect. $\endgroup$
    – Sebastian Reichelt
    Sep 1, 2015 at 21:05
  • $\begingroup$ I see, you are right! $\endgroup$
    – Yair Daon
    Sep 2, 2015 at 13:43
2
$\begingroup$

The closer the sampling frequency gets to half the (highest spectral) frequency of the signal, the longer one has to sample to get a sufficient bounded fit (for some bound) of a reconstruction using those finite number of samples. This required sampling time extends to infinity as one approaches Fs/2.

So, you simply didn't sample for long enough in your experiment. Try 5 thousand and 5 million samples for sampling a 12.9 Hz sinusoid at 26 sps, and you will see a closer fitting looking plot over most of an X axis of 5 thousand or better million reconstructed sample points. Then try a sinusoid even closer to Fs/2 and you will find you need to sample even longer to get the same visual quality of fit.

Improve this answer
$\endgroup$
2
  • $\begingroup$ thanks, I feel enlightened! Is this a theorem with known bounds or just a heuristic practitioners use? If the former is right, I'd love to see a reference. $\endgroup$
    – Yair Daon
    Sep 1, 2015 at 17:16
  • 1
    $\begingroup$ Note that the sampling theorem applies to band-limited signals, and perfectly band-limited signals (transform with finite support) have to be infinite in length. Thus any finite length sample vectors will not be band-limited, but will have "leakage" from the windowing that can cause aliasing artifacts. The shorter the time domain window, the more leakage there is in the frequency domain, especially good for creating aliasing artifacts near Fs/2. $\endgroup$
    – hotpaw2
    Sep 1, 2015 at 18:25

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.