FFT of BPSK signal interpretation

Question

I am plotting spectrum of BPSK signal just series of -1 and 1, and I can see random peaks (looks to me to be random) at various locations but I am really confused if I do not multiply my signal by carrier shouldn't the peak of bpsk modulated signal be at DC? If not then how to identify the exact location that peaks should be ? Here is what I have done:

bits=randi([0 1],1,1000);
bpsk_mod=2*bits-1;
fft_bpsk=fft(bpsk_mod);
len=length(bpsk_mod);
f=[-len/2:len/2-1];
plot(f,abs(fftshift(fft(bpsk_mod)))

Answer 1

A random sequence of +1 and -1 is not really a BPSK signal; it's just a sequence of bits. As such, its FFT has not much meaning, at least not in the sense you're thinking of.

To create an actual BPSK signal, first denote your sequence of $N-1$ bits by $a_k$, where $k=0,\ldots,N-1$ is an integer. Second, choose a bit rate $R_p=1/T_p$. Then, choose a pulse shape, $p(t)$. The pulse shape is typically rectangular, a sinc, a raised cosine, or similar. The important thing is that the pulse is orthogonal to $T_p$-shifted versions of itself; that is, $$\int_{-\infty}^\infty p(t-mT_p)p(t-nT_p)=0$$ for all integers $m\neq n$. Now you're ready to create your BPSK signal, which is given by $$s(t)=\sum_{k=0}^{N-1} a_kp(t-kT_p).$$

Now you can find the FT of $s(t)$ (or the DFT of a properly sampled version of $s(t)$), and you can find the bandwidth and other spectral properties of your BPSK signal.

Answer 2

What you see is just a bad estimate of the spectral characteristic of your random number generator. Ideally, if the random samples were independent, the corresponding power spectrum should be flat (indicating white data).

Note that this is generally not the spectrum of a BPSK signal, because a BPSK signal is a continuous-time signal, and its spectrum also depends on the pulse shape, as pointed out in MBaz's answer.

Without any further knowledge, if you just have discrete-time samples, and if you can be sure that your desired signal ($\pm 1$) has not been scaled in any way, then the easiest method for estimating the SNR is to estimate the noise variance by first deciding on the actual data using a slicer, and then computing the estimated noise variance as

$$\sigma^2_n=\frac{1}{N-1}\sum_{k=0}^{N-1}(a[k]-\text{sign}(a[k]))^2\tag{1}$$

where $N$ is the number of noisy samples $a[k]$. The $\text{sign}$ function is $1$ if its argument is positive and $-1$ if it is negative. This is just a slicer to estimate the original clean data. Note that I've assumed that the noise is zero mean (which is not necessary, just easier). The SNR is then simply given by $1/\sigma_n^2$.

As with all decision-directed methods, this one will only work for high to moderate SNRs. For very low SNRs, the slicer makes wrong decisions, and, consequently, the noise estimate $a[k]-\text{sign}(a[k])$ used in $(1)$ is wrong.

Answer 3

You are quite right. BPSK signal, sampled at Nyquist rate, is indeed just a sequence of independently generated -1 and 1 - so you got that right. For convenience, let us denote this discrete time sequence by $x_n$ ($n$ denotes time domain).

Your second observation was that the frequency transform was taking up the entire spectrum, and had peaks at random frequency locations. That is expected as well. Below, I explain why this happens and how to estimate the power spectrum.

Consider the time sequence $x_n$. The power spectrum gives you how much energy it occupies in different frequencies. To calculate the power spectrum, we start with the time domain autocorrelation and take its Fourier transform.

Since the sequence $x_n$ is independent, its autocorrelation $R(\tau)$ has energy only at the zeroth term i.e. $$R(\tau) = E[x_n x_{n+\tau}^*] = \left\{ \begin{array}{cc} 1 &if \tau=0 \\ 0 & if \tau \neq 0 \end{array} \right.$$

Denoting the power spectrum of $x_n$ by $S(f)$, we have the following $$ S(f) = \sum_{k=-\infty}^{+\infty} R(k) e^{-j2\pi f k T_s} \\ = R(0) e^{-j2\pi f \times 0 \times T_s} \\ = 1 $$

Here, $T_s$ is the sampling period associated with the discrete time sequence $x_n$, and the frequency range is $f \in [-1/{2T_s}, 1/{2T_s}]$. The second step in the equation above makes use of $R(k)=0$ if $k\neq 0$.

So, we expect $x_n$ to have a flat frequency response; i.e. all frequencies are expected to have equal energy. We do not expect energy to peak at DC for this input sequence.

However, in your Matlab code, you observed randomly located peaks in the frequency transform. This happens because you took the FFT of $\{x_1,\cdots,x_{1000}\}$; i.e. this is because you are looking at the frequency response of a specific instance of a random signal.

The power spectrum does not really have these peaks. To get the power spectrum, you need to do some additional steps. What you need is the average energy of the frequency response. I suggest adding an outer loop to your code, as below, and averaging the energies you measure in each frequency.

num_iter = 10000;
sample_len_per_iter = 1000;
power_spectrum = zeros(1,sample_len_per_iter);
for iter_index = 1:num_iter
    bits=randi([0 1],1,sample_len_per_iter);
    bpsk_mod=2*bits-1;
    fft_bpsk=fft(bpsk_mod);
    power_spectrum = power_spectrum + abs(fft_bpsk.^2); % accumulate power
end
power_spectrum = power_spectrum / num_iter; % average out 
f=[-sample_len_per_iter/2:sample_len_per_iter/2-1];
plot(f,power_spectrum)

The above should give you a much more flat response.

In case you are curious, the fft output of the specific instance of BPSK sequence is (or converges to) an i.i.d. complex Gaussian sequence. This is a consequence of the central limit theorem. This is why you observed randomly located peaks.

In a practical system, this signal is up-sampled and filtered by a baseband pulse shape so that the bandwidth actually occupied meets some pre-specified requirement of RF emissions.

Edit about the power scaling:

The fft() function in Matlab includes a scaling of $1/\sqrt{L}$ where L is the fft length i.e. Length of fft_bpsk. This, after taking magnitude square, gives the desired power scaling of $1/L$. $$X_k = \frac 1 {\sqrt L} \sum_{n=1}^{L} e^{-j2\pi \frac{(n-1)(k-1)}{L}} x_n\\ E[|X_k|^2] = \frac 1 L \sum_{n=1}^{L} \sum_{m=1}^{L} e^{-j2\pi \frac{(n-m)(k-1)}{L}} E[x_n x_m^*] \\ = \frac 1 L L R(0) = 1 $$