I am simulating a first order delta sigma modulator (DSM) used as an ADC. The input to my DSM is a full scale sine wave. The power spectral density I expect from the DSM output is as follows : The sinusoid peak as well as some noise that is shaped to the higher frequencies.
In order to avoid any spectral leakage, I define my sinusoid from the FFT bin it will occupy as well as the number of samples for the FFT (this is also suggested for simulations in the book "Understanding delta sigma data converters" by R. Schreier and G. C. Temes).
bin = 15
Nfft = 8192
t = np.arange(0, Nfft)
signal = np.sin(2*np.pi*bin/Nfft*t)
Then I pass it through the DSM, for which the code is as follows :
def mod1(signal):
Nfft = signal.size
y = np.zeros(Nfft)
x = np.zeros(Nfft)
v = np.zeros(Nfft)
for i in range(1, Nfft):
v[i] = np.sign(y[i-1])
x[i] = signal[i] - v[i]
y[i] = x[i-1] + y[i-1]
return v
Here is my issue :
When I choose Nfft
as a multiple of the frequency bin, I will see harmonics of sine frequency appearing. These do not appear when Nfft
is not multiple of the frequency bin.
If these harmonics should really appear and are not some artifact, as they are placed at k*bin
, they should also fall exactly into a bin and avoid spectral leakage regardless of the number of fft samples.
Here is a plot illustrating this issue. On the top left, a sine wave and its DSM modulation. On the bottom left, their PSD when the number of samples is not a multiple of the frequency bin. On this side, the noise shaping can clearly be seen.
The same plots can be found on the right when the number of samples is a multiple of the frequency bin. On this side, the harmonics can be seen.
So which PSD is more correct and why do these harmonic appear only in the particular case of the number of FFT samples being a multiple of the fundamental frequency bin ?
Here is the full code to perform the plots :
import numpy as np
import matplotlib.pyplot as plt
import scipy.signal as sg
def generate_sine_from_bin(bin, Nfft):
t = np.arange(0, Nfft)
signal = np.sin(2*np.pi*bin/Nfft*t)
return signal
bin1 = 15
Nfft1 = 8192
bin2 = 15
Nfft2 = 9000
signal1 = generate_sine_from_bin(bin1, Nfft1)
signal2 = generate_sine_from_bin(bin2, Nfft2)
def mod1(signal):
Nfft = signal.size
y = np.zeros(Nfft)
x = np.zeros(Nfft)
v = np.zeros(Nfft)
for i in range(1, Nfft):
v[i] = np.sign(y[i-1])
x[i] = signal[i] - v[i]
y[i] = x[i-1] + y[i-1]
return v
v1 = mod1(signal1)
v2 = mod1(signal2)
ws1 = sg.windows.hann(signal1.size, sym = True)
ws2 = sg.windows.hann(signal2.size, sym = True)
def compute_psd(signal, w):
w1 = np.linalg.norm(w, 1)
V =np.abs(np.fft.fft(signal*w)/(w1/2))**2
return V
Ps1 = compute_psd(signal1, ws1)
Ps2 = compute_psd(signal2, ws2)
P1 = compute_psd(v1, ws1)
P2 = compute_psd(v2, ws2)
fs1 = np.linspace(0, 1, Nfft1)
f1 = np.linspace(0, 1, Nfft1)
fs2 = np.linspace(0, 1, Nfft2)
f2 = np.linspace(0, 1, Nfft2)
fig, ax = plt.subplots(2, 2, figsize = (10, 5))
ax[0, 0].set_title(f"sine wave bin={bin1} Nfft={Nfft1}")
ax[0, 0].plot(signal1)
ax[0, 0].plot(v1)
ax[0, 0].legend(["signal", "dsm"])
ax[0, 0].set_xlabel("Samples")
ax[1, 0].set_title("DSM PSD (Hann window)")
ax[1, 0].semilogx(fs1, 10*np.log10(Ps1))
ax[1, 0].semilogx(f1, 10*np.log10(P1))
ax[1, 0].legend(["signal", "dsm"])
ax[1, 0].set_xlabel("Normalized freq")
ax[1, 0].set_ylim([-100, 50])
ax[0, 1].set_title(f"sine wave bin={bin2} Nfft={Nfft2}")
ax[0, 1].plot(signal2)
ax[0, 1].plot(v2)
ax[0, 1].legend(["signal", "dsm"])
ax[0, 1].set_xlabel("Samples")
ax[1, 1].set_title("DSM PSD (Hann window)")
ax[1, 1].semilogx(fs2, 10*np.log10(Ps2))
ax[1, 1].semilogx(f2, 10*np.log10(P2))
ax[1, 1].legend(["signal", "dsm"])
ax[1, 1].set_xlabel("Normalized freq")
ax[1, 1].set_ylim([-100, 50])
plt.tight_layout()
plt.show()
Thank you for your help !