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:

enter image description here


2 Answers 2


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).

FFT results

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
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))
  • $\begingroup$ Did you use an open source library to generate the maximum length LRS? I am not able to reproduce your results using scipy.signal.max_len_sequence(7,taps=[0,1,7])[0]*2 -1.0. I have tried longer and shorter sequences and tried different primitive polynomials from here: ece.unb.ca/tervo/ece4253/polyprime.shtml $\endgroup$ Oct 16, 2023 at 20:20
  • $\begingroup$ Hi BigBrown, I used sig.max_len with the default polynomial. $\endgroup$ Oct 16, 2023 at 20:31
  • 1
    $\begingroup$ I’ll add my code later $\endgroup$ Oct 16, 2023 at 20:46
  • $\begingroup$ where is the pshape function from? $\endgroup$ Oct 17, 2023 at 10:34
  • 1
    $\begingroup$ it is a custom function that does the interpolate by 4 with raised cosine pulse shaping, I didn't think you needed to see that so I could keep it short here. The pulse shaping filter I used has a length of 24 symbols or 96 taps with $\alpha=0.3$. $\endgroup$ Oct 17, 2023 at 12:02
  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.

  • $\begingroup$ It is a PSD in the traditional DSP sense, time averaged FFTs. More specifically, it's the output of matplotlib's psd function $\endgroup$ Oct 9, 2023 at 14:19
  • $\begingroup$ in the traditional PSD sense, a PSD is a property of a stochastic process :D Some programming languages are a bit sloppy in their documentation, but this should really be listed as PSD estimate. $\endgroup$ Oct 9, 2023 at 14:41
  • $\begingroup$ as the broken language of the documentation of that function says : by Welch's average periodogram method; this is an estimator! no function applied to data could ever truly know the PSD :) So, the amount of data you offered for this calculation was simply not enough for a very stable estimate, it seems. $\endgroup$ Oct 9, 2023 at 14:43
  • $\begingroup$ I'm not really following what you are saying. 1 million samples to produce a PSD is WAY more than enough, and this is not due to a lack of samples. I can use far less samples with other simulated waveforms and produce a flat spectrum. The difference is the underlying data is not random and the spreading code is not random, among some other implementation specifics. Point being, this is not related to the number of points used for the PSD estimate. $\endgroup$ Oct 9, 2023 at 14:46
  • $\begingroup$ Why are 1 Million samples "WAY more" than sufficient? Look at the "noise-only" out-of band: clearly, your estimate is still quite noisy. Your PSD estimate seems fairly high-resolved, and as 1. hints at, your PSD is the PSD of the pulse shaping filter multiplied with the PSD of the data, and as I said, the data isn't sufficient – in it's amount or its whiteness. So make sure you use a parametrization of the psd function that fits your needs well, and that your data is sufficiently large and uncorrelated. $\endgroup$ Oct 9, 2023 at 14:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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