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

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.

| improve this answer | |
  • $\begingroup$ The problem is I have received noisy signal (-1 and 1 plus noise)and I would like to estimate its SNR in frequency domain, therefore I need to know what range of frequencies corresponds to signal . Yes you are right if I design pulse shaping and do it in continuous domain then everything makes sense but all I have is a noisy signal. $\endgroup$ – justin Aug 12 '16 at 20:14
  • $\begingroup$ @justin Let $r(t)=s(t)+n(t)$ be the received signal. In a receiver, the signal $r(t)$ goes through a matched filter. Let $u(t)$ be the matched-filtered signal. What you have is the sequence $u[kT_p]$, but what you need to do your frequency analysis is the signal $u(t)$ itself. $\endgroup$ – MBaz Aug 12 '16 at 22:44
  • 1
    $\begingroup$ You are better off estimating the SNR in the time domain in this situation because the noise (assuming it is white) can't be isolated in the spectrum. It overlaps with the signal. You can get a time domain estimate by taking the mean and variance of the absolute value of the time samples. $\endgroup$ – John Aug 13 '16 at 2:44

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


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.

| improve this answer | |
  • $\begingroup$ Thanks. Indeed this works well when SNR is around 5 or 6 dB or higher. But it does not work for low SNR values Do you know any technique when SNR is 0 to 2 dB. $\endgroup$ – justin Aug 15 '16 at 21:06
  • $\begingroup$ @justin: You can't do much if the SNR is bad. If you receive a sample value of $0.2$, the noise value can be either $-0.8$ (if $+1$ was sent), or it can be $+1.2$ (if $-1$ was sent). How do you want to tell which of the two cases is true? Of course, you could heuristically choose some threshold and discard all noisy data with a magnitude smaller than the threshold, and then apply formula $(1)$. This would get rid of data for which the sign function is prone to lead to wrong decisions. $\endgroup$ – Matt L. Aug 15 '16 at 21:19

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);
    power_spectrum = power_spectrum + abs(fft_bpsk.^2); % accumulate power
power_spectrum = power_spectrum / num_iter; % average out 

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 $$

| improve this answer | |
  • $\begingroup$ Thanks. But shouldn't we have a sum before abs(fft_psk.^2) so it should be sum(abs(fft_bpsk.^2)) and also isn't the power defined as sum of squares divided by length of the signal. So shouldn't the overall expression be sum(abs(fft_bpsk.^2))/length(fft_bpsk); $\endgroup$ – justin Aug 15 '16 at 20:18
  • $\begingroup$ 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 right scaling. So three second scale that you ask is being taken care of already. The first question about sum: no, we don't want to sum it up. We want the power of signal at each frequency position $\endgroup$ – Vignesh Aug 16 '16 at 4:51

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.