0
$\begingroup$

Start with a 2 Hz signal. The signal is sampled at a 4.167 Hz sample rate.

enter image description here

The intent is to reconstruct a sampled signal (top right side) to be identical to the original signal (top left side).

Insert 200 zeros between the max frequency in the middle of the frequency spectrum (middle right side image). Then use Inverse DFT to attempt to reconstruct the original signal.

enter image description here

The reconstructed signal (bottom right side) is similar, but not identical to the original signal (top left side).

Increasing the number of sine wave cycles being sampled, and changing the sample delta X from 0.01 sec to 0.001 sec, does not eliminate the problem.

enter image description here

It should be possible to reconstruct the original, smooth sine wave (according to Nyquist).

Why does this not work perfectly?

Perhaps there are minor errors in the DFT code equations.

DFT:

public Complex[] acplxConvertSignalToDft(double[] adSignalYValues)
{
  int iNumOfValues = adSignalYValues.Length;
  Complex[] acplxFrequencies = new Complex[iNumOfValues];
  for (int ii = 0; ii < iNumOfValues; ii++) {
    acplxFrequencies[ii] = 0;
    for (int iii = 0; iii < iNumOfValues; iii++) {
      acplxFrequencies[ii] += adSignalYValues[iii] *
      Complex.Exp(-Complex.ImaginaryOne *
      2 * Math.PI * (ii * iii) /
      Convert.ToDouble(iNumOfValues));
    }
    acplxFrequencies[ii] = acplxFrequencies[ii] / iNumOfValues;
  }
  return acplxFrequencies;
}//acplxConvertSignalToDft

Inverse DFT:

public double[] adConvertDftToSignal(Complex[] acplxFrequencies)
{
  //---Number of spectrum elements
  int iNumOfValues = acplxFrequencies.Length; 
  double[] adInverseDftSignalYValues = new double[iNumOfValues];
  for (int ii = 0; ii < iNumOfValues; ii++) {
    Complex cplxSum = 0;
    for (int iii = 0; iii < iNumOfValues; iii++) {
      cplxSum += acplxFrequencies[iii] *
                 Complex.Exp(Complex.ImaginaryOne *
                 2 * Math.PI * (iii * ii) /
                 Convert.ToDouble(iNumOfValues));
    }
    //---As a result we expect only real values (if our calculations 
    //    are correct, imaginary values should be equal or close to zero).
    adInverseDftSignalYValues[ii] = cplxSum.Real;
  }
  return adInverseDftSignalYValues;
}//adConvertDftToSignal

This DFT code comes from Jakub Szymanowski => https://www.codeproject.com/Articles/1077529/Fourier-Transform-in-Digital-Signal-Processing.

One key comment in the Inverse DFT code: "As a result we expect only real values (if our calculations are correct, imaginary values should be equal or close to zero)." Maybe the imaginary values cannot be ignored.

In addition, the frequency spectrum of the sine wave (middle left image) shows side bands next to the 2 Hz frequency point (easier to see in the first image at the top of this post). Perhaps there is a minor error in the DFT code equations, as only the 2 Hz point should appear in a perfect case.

$\endgroup$
4
  • $\begingroup$ iiuc you are interpolating the time-domain signal by zero-padding in the FFT domain. Have a look at this dspguru.com/dsp/howtos/… $\endgroup$
    – ZR Han
    Nov 12 '21 at 7:29
  • 1
    $\begingroup$ The point that you've added is the DC (zero frequency point). It definitely should NOT be replicated at the top of the frequency range. The original is perfectly symmetrical about the DC point. $\endgroup$
    – Peter K.
    Nov 12 '21 at 13:26
  • $\begingroup$ @ttom: I've re-worded your question in a manner that I think gets at the heart of the question you're asking. Please feel free to revert it if I've missed the mark. $\endgroup$
    – Peter K.
    Nov 13 '21 at 0:29
  • $\begingroup$ Thanks Peter K. for clarifying my original post! $\endgroup$
    – ttom
    Nov 14 '21 at 13:23
0
$\begingroup$

So, I've had a bit of a think about this. I believe the problem is that your original 201 point signal:

import matplotlib.pyplot as plt
import numpy as np

T = 201
fs = 100
t = np.arange(T)/fs

x = np.sin(2*np.pi*2*t)
x_fft = np.fft.fft(x)

is not bandlimited.

So when you decimate it by a factor of 24

decimation_factor = 24
x_decimated = x[1:T:decimation_factor]
t_decimated = t[1:T:decimation_factor]

you're going to get aliasing.

Original and subsampled signal

Resampled signal (and original and subsampled)

The FFT of the original 201 point signal on a log scale is below.

Log scale FFT

As you can see, the values in between the peaks are non-zero.

To look at this even more closely:

Zoom in on peak of FFT

the red line shows the new Nyquist frequency (fs/2). As you can see, the frequency content of the original signal is non-zero above this frequency. So subsampling at this new sampling frequency will result in aliasing.

I'm still working on this to find a signal that can be perfectly reconstructed after subsampling to this level. Watch this space for another update.

$\endgroup$
3
  • $\begingroup$ Is that because the original sine wave is not perfect? A delta X of 0.001 is better than a delta X of 0.01? For 0.001, there are 2001 points. $\endgroup$
    – ttom
    Nov 18 '21 at 20:01
  • 1
    $\begingroup$ After thinking more about this, the sampled x and y values are precise, not an approximation, although the accuracy of the computer will never result in complete precision. Per Nyquist, it should be possible to identify the frequencies => create the original sine wave, regardless of the delta X used. Is there a slight error in the DFT equations? $\endgroup$
    – ttom
    Nov 18 '21 at 21:12
  • $\begingroup$ @ttom Here's a minor update. Still doesn't have an example of something that could be reconstructed perfectly, but explains why your example can't be. $\endgroup$
    – Peter K.
    Nov 24 '21 at 17:39

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.