DSSS spectrum not flat

Question

Why is the PSD of the simulated BPSK/DSSS in AWGN not flat? Here are the simulation parameters:

• BPSK symbols are randomly generated and are {+1,-1}
• Spreading factor = 128
• Spreading code is a randomly generated sequence {+1,-1} that is the length of the spreading factor.
• Each BPSK symbol is spread by the spreading code
• Signal is oversampled by 4 using a raised cosine filter
• AWGN is added
• ~1M samples used in simulation

The PSD of the signal is:

Answer 1

The OP may be seeing more variability than expected due to the quality of the random number generator in creating a uniform spectrum. Below shows two spectrums I created with similar parameters as the OP but using two different spreading sequences.

"randint": length 128 +1/-1 sequence generated from Python's numpy random module.

"max-len": length 127 maximum length LFSR based pseudo-random sequence.

Below shows the spectrums after spreading a random BPSK data sequence (also generated similar to randint), and interpolating by 4 with a raised cosine pulse shaping filter. I did not add noise; the result would trivially be the superposition of the noise and waveform spectrums, and I didn't want to mask the characteristics of the waveform itself as generated.

The length of the randint and max-len spread waveforms were 1,024,000 samples and 1,016,000 samples respectively.

Further insight is gained by reviewing the FFT of the spreading sequences repeated with no data modulation (basically continuously sending "1" for each BPSK symbol). This reveals the spectrum of the unmodulated spreading sequence, and for that we would expect to see discrete tones spaced by the rate at which the entire sequence repeats. It is exactly this spectrum that would convolve with the much narrower spectrum of the data modulation (multiplication in time is convolution in frequency). The spectrum of the data modulation, if it were perfectly random, would approach a Sinc magnitude response, with the width from the peak to first null equal to the spacing between the tones of the spectrum for the repeated spreading sequence.

We see for the case on the left where randint was used as the spreading sequence, a variability reaching 20 dB from tone to tone in the spectrum. This indicates a poor quality of "white-ness" for the random sequence used. When the randint sequence is modulated with the data, the convolution of the data spectrum with the randint spectrum will follow this variability.

In comparison on the right is the very flat (white) spectrum for the max-len case. The dominant contributor to the remaining variability in the initial spectrum that I presented with using the max-len sequence is due using randint for the BPSK data pattern.

Ultimately a longer time duration will cause the spectral contributions for the data itself to become increasingly flatter assuming random data, but that will not help with the spreading sequence as demonstrated with the FFT plots below (the same plot will result for the unmodulated sequence regardless of long we repeat it for- so the quality is in the sequence of the 127 or 128 digits used).

Regardless, the spread of the estimate of the signal's PSD will always be larger than the spread of the estimate for the noise PSD alone, given in the OP's case we have 1 million independent samples of noise. In contrast, we only have about 1950 independent data samples (given the interpolation by 4 and then spreading by 128, a total factor of 512). When averaging white noise processes, the standard deviation of the average goes down as the square of the number of samples in the average and therefore we would expect to have about 27 dB less variability in the PSD of the noise versus the PSD of the data (to get the data spectrum as clean as the noise spectrum that we see, we would need 512 million total samples!).

Python code used:

import numpy as np
import numpy.random as rand
import scipy.signal as sig
import matplotlib.pyplot as plt

# spreading_codes
sf = 128
rng = np.random.default_rng()
randint = 2*rng.integers(low=0, high=2, size = 128) - 1
max_len = sig.max_len_seq(7)[0] * 2 - 1

# create random BPSK data
nsamps = 200000
data = 2 * rng.integers(low=0, high=2, size = nsamps) - 1

spread_r = (data[:, None] * randint).flatten()  # spread with randint
spread_m = (data[:, None] * max_len).flatten()  # spread with max-len

pulse_shaped_r = pshape(spread_r)    # interpolate by 4 and pulseshape
pulse_shaped_m = pshape(spread_m)    # interpolate by 4 and pulseshape


# PSD plot
plt.figure()
f, psd_r= sig.welch(pulse_shaped_r,return_onesided=False, nperseg=1024)
f2, psd_m= sig.welch(pulse_shaped_m,return_onesided=False, nperseg=1024)
plt.plot(fft.fftshift(f), 10*np.log10(fft.fftshift(psd_r)), label= 'randint')
plt.plot(fft.fftshift(f2), 10*np.log10(fft.fftshift(psd_m)), label='max-len')

## FFTs with no data modulation to evaluate whiteness of random sequence
data2 = np.ones(100)
repeat_r = (data2[:, None] * randint).flatten()
repeat_m = (data2[:, None] * max_len).flatten()

pulse_repeat_r = pshape(repeat_r)
pulse_repeat_m = pshape(repeat_m)

fout_r = fft.fft(pulse_repeat_r)/len(pulse_repeat_r)
freq_r = fft.fftfreq(len(pulse_repeat_r))
fout_m = fft.fft(pulse_repeat_m)/len(pulse_repeat_m)
freq_m = fft.fftfreq(len(pulse_repeat_m))

Answer 2

1. "Signal is oversampled by 4 using a raised cosine filter" remember what that means for the spectrum!

2. what you show is not a PSD, it's a PSD estimate, the difference is that a PSD would be a property of the infinitely extending stochastic signal, and you're showing us something that looks like an FFT of a short observation.