The Cyclic Autocorrelation for Rectangular-Pulse BPSK

The cyclic autocorrelation for rectangular-pulse BPSK can be derived as a relatively simple closed-form expression (see My Papers [6] for example or The Literature [R1]). It can be estimated in a variety of ways, which we will discuss in future posts. The non-conjugate cycle frequencies for the signal are harmonics of the bit rate, $k f_{bit}$, and the conjugate cycle frequencies are the non-conjugate cycle frequencies offset by the doubled carrier, or $2f_c + k f_{bit}$.

Recall that the simulated rectangular-pulse BPSK signal has $10$ samples per bit, or a bit rate of $0.1$, and a carrier offset of $0.05$, all in normalized units (meaning the sampling rate is unity). We’ve previously selected a sampling rate of $1.0$ MHz to provide a little physical realism. This means the bit rate is $100$ kHz and the carrier offset frequency is $50$ kHz. From these numbers, we see that the non-conjugate cycle frequencies are $k 100$ kHz, and that the conjugate cycle frequencies are $2(50) + k 100$ kHz, or $100 + k 100$ kHz.

The blindly estimated CAF for our noisy rectangular-pulse BPSK signal is shown here:

The non-conjugate CAF for $\alpha = 0$ is the conventional autocorrelation function. Here we can see it is triangular in shape, with an additional inflation of the peak for $\tau = 0$, due to the presence of noise, which has non-zero power. For comparison, here is the numerically evaluated theoretical formula for the CAF for rectangular-pulse BPSK:

The match is excellent.

The theoretical formulas for the cyclic autocorrelation and spectral correlation function for BPSK signals (and other digital QAM/PSK signals) can be found in several places. One of the first places it was published is in the book Statistical Spectral Analysis by W. A. Gardner (The Literature [R1], Chapter 12). The formulas for both the $n$th-order reduced-dimension cyclic cumulant and the $n$th-order cyclic polyspectrum for PAM/PSK/QAM can be found in My Papers [6] and my dissertation. The $n$-th order cyclic cumulant reduces to a lag-shifted version of the cyclic autocorrelation for $n=2$.

Interpreting the CAF

The non-conjugate CAF can be interpreted as the correlation between the lag product $x(t+\tau/2)x^*(t-\tau/2)$ and the complex-valued sine wave $e^{i 2 \pi \alpha t}$. From that point of view, the CAF $R_x^\alpha(\tau)$ is the complex amplitude of the additive sine wave that is present in the lag product. That is, that sine wave has amplitude $| R_x^\alpha(\tau)|$ and phase $\angle R_x^\alpha(\tau)$. Note that the amplitude and phase depend on the delay $\tau$. When the CAF is zero, there is no finite-strength additive sine-wave component in the lag product for the value of $\tau$ chosen; otherwise there is.

The conjugate CAF can be interpreted as the correlation between the lag product $x(t+\tau/2)x(t-\tau/2)$ and the complex-valued sine wave $e^{i 2 \pi \alpha t}$.

Some signals may have only non-zero non-conjugate CAF values (such as large-alphabet digital QAM), some only non-zero conjugate CAF values (such as analog amplitude modulation, discounting the ever-present non-conjugate cycle frequency of zero), and some have both (such as our favorite signal, the rectangular-pulse BPSK signal).

Comparison with a Bandwidth-Efficient BPSK Signal

The rectangular-pulse signal has infinite bandwidth, and therefore it possesses an infinite number of cycle frequencies. However, the strength of the cyclic autocorrelation for most of those cycle frequencies is very small, so that in practical terms, the rectangular-pulse signal possesses ten or so significant features. That’s actually quite a lot compared to practical real-world signals like BPSK with square-root raised-cosine (SRRC) pulses. Here is the estimated CAF for a BPSK with SRRC pulses and a pulse roll-off factor of 0.3 (the excess bandwidth is 30%):

The non-conjugate CAF possesses exactly three cycle frequencies, $\{-f_{bit}, 0, f_{bit}\}$ (only two of those are shown in the plot), and the conjugate CAF also possesses three cycle frequencies $\{2f_c-f_{bit}, 2f_c, 2f_c+f_{bit}\}$. Note also that the width of the CAF is much larger than for the rectangular-pulse signal (the $\tau$ axes are different between the plots for the two signals). This is a direct consequence of the fact that the rectangular-pulse signal is strictly time-limited whereas the SRRC-pulse signal is strictly bandlimited.

I'm a signal processing researcher specializing in cyclostationary signal processing (CSP) for communication signals. I hope to use this blog to help others with their cyclo-projects and to learn more about how CSP is being used and extended worldwide.

33 thoughts on “The Cyclic Autocorrelation for Rectangular-Pulse BPSK”

1. wahiba says:

hello, could you please explain why cyclic frequencies in your example are equal to k* (0.1) MHz?

1. please I need an example code for cyclostastionary features detection for bpsk modulation

1. Hi Ben. The CSP Blog is a self-help blog. I don’t give out much code. Let me know if you get stuck, though, trying to write your own. That’s the only way you’ll really understand CSP.

1. Ben mohamed says:

this is my own code but i don’t know if it’s good to calculate the cyclic autoucorrelation function or not;

function [caf2, taus, alphas] = caf2(s)
N1=length(s);
N=floor(N1/3);
caf2 = zeros(2*N+1,N);
prods = zeros(2*N+1,N);
taus=-N:N;
alphas=-pi+2*pi/(N):2*pi/(N):pi;
alphas=alphas/pi;
for tau = taus
for alpha=alphas
stat = (s(N+1:2*N)).*exp(-j*pi*alpha*t(N+1:2*N));
shifted = conj(s(N+1-tau:2*N-tau)).*exp(j*pi*alpha*t(N+1-tau:2*N-tau));
prod=stat.*shifted.*exp(-2*j*pi*alpha*t(N+1:2*N));
prods(tau+N+1,:)=prod;
end
caf2(tau+N+1,:)=fftshift(fft(prod));

end
end

2. Your code does not run to completion:

>> caf2(s)
Undefined function or variable ‘t’.

Error in caf2 (line 11)
stat = (s(N+1:2*N)).*exp(-j*pi*alpha*t(N+1:2*N));

So you need to fix that first.

I see that you are trying to compute a statistic for a large number of cycle frequencies and delays. I would like to suggest the following: Divide and conquer. Do the simplest things first, then build up to the most complicated. For example, try to estimate the cyclic autocorrelation for just a single cycle frequency (say, zero), and a few delays (tau) near zero. Then try to extend to all delays. Then extend so that you can specify any cycle frequency. Finally, when that is all correct, put it in a loop and call it for as many cycle frequencies as you want.

2. ASMAA says:

hey
i would like to know about the theorical formula for bpsk cyclic autocorrelation ?
is there any listed reference?

1. I’ve updated the post with references that should help you find the answer to your question. Let me know if you still have problems or questions.

3. Ben mohamed says:

thanks ser
in witch frequency we have peak

Fs = 5000;
t = 0:1/Fs:1-1/Fs;
s = (cos(2*pi*600*t)+randn(size(t)));
N1=length(s);
N=floor(N1/3);
caf2 = zeros(2*N+1,N);
prods = zeros(2*N+1,N);
taus=-0.5:0.5;
alpha=0.2;
stat = (s(N+1:2*N)).*exp(-j*pi*alpha*t(N+1:2*N));
shifted = conj(s(N+1-taus:2*N-taus)).*exp(j*pi*alpha*t(N+1-taus:2*N-taus));
prod=stat.*shifted.*exp(-2*j*pi*alpha*t(N+1:2*N));

caf2=(fft(prod));

1. Ben, that code also does not run correctly:
>> caf2_2
Subscript indices must either be real positive integers or logicals.

Error in caf2_2 (line 11)
shifted = conj(s(N+1-taus:2*N-taus)).*exp(j*pi*alpha*t(N+1-taus:2*N-taus));

Is it running without errors on your system? I think you should at least make sure the code runs before sending it to someone else.

4. Ben Mohamed says:

thanks mister but my code run without errors

caf2 not caf_2

Fs = 5000;
t = 0:1/Fs:1-1/Fs;
s = (cos(2*pi*600*t)+randn(size(t)));
N1=length(s);
N=floor(N1/3);
caf2 = zeros(2*N+1,N);
prods = zeros(2*N+1,N);
taus=-0.5:0.5;
alpha=0.2;
stat = (s(N+1:2*N)).*exp(-j*pi*alpha*t(N+1:2*N));
shifted = conj(s(N+1-taus:2*N-taus)).*exp(j*pi*alpha*t(N+1-taus:2*N-taus));
prod=stat.*shifted.*exp(-2*j*pi*alpha*t(N+1:2*N));
caf2=fftshift(fft(prod));
plot(abs(caf2))

1. The code you posted on April 4 produces the following output using MATLAB R2015b:

>> bm_code_040417
Subscript indices must either be real positive integers or logicals.

Error in bm_code_040417 (line 11)
shifted = conj(s(N+1-taus:2*N-taus)).*exp(j*pi*alpha*t(N+1-taus:2*N-taus));

The code you posted on April 7 produces the following output:

>> bm_code_040717
Subscript indices must either be real positive integers or logicals.

Error in bm_code_040717 (line 11)
shifted = conj(s(N+1-taus:2*N-taus)).*exp(j*pi*alpha*t(N+1-taus:2*N-taus));

I’m pretty sure the problem is with

taus=-0.5:0.5;

which yields

>> taus

taus =

-0.500000000000000 0.500000000000000

So you are mixing fractions with integers in your indexing.

This isn’t the kind of debugging I normally help with as it doesn’t have anything to do with understanding CSP.

5. Stephen Sun says:

I use rectangle pulse for BPSK. But the CAF shows followings: 1)there are only 3 lines (represent alpha=-2fc, 0, 2fc)when delay tau=0, and it do nothing with symbol rate; 2)CAF quickly decreases and tends to constant as delay increases when alpha is 2fc. Are those right? Why can’t find cycle frequencies related with symbol rate? Thanks in advance.

1. Did you use your own rectangular-pulse BPSK signal? What happens if you use the one I provided at the end of the post on rectangular-pulse BPSK? See the end of that post for details.

1. Stephen Sun says:

Thanks for your reply. After using your BPSK signal code, the CAF shows more peaks. But it’s still a big difference from your diagram. I have two questions here:
1) What is the CA formula of BPSK signal? I only find the cyclic spctrum formula in the book Statistical Spectral Analysis; 2) Should s[t] multiply by e^-j*pi*alpha*t before calculating delay product, like code written by Ben Mohamed aboved?

1. You can look at the references I cite for the formula for the cyclic autocorrelation for PSK/QAM (one of the references
is my doctoral dissertation, which can be downloaded from the CSP Blog here). Or, look at formula (4) in the post on digital QAM and PSK. Specialize that formula to the case of $n=2$. That formula is slightly more general than the formula for the “symmetrical” cyclic autocorrelation. The symmetrical autocorrelation uses the lags $-\tau/2$ and $\tau/2$, whereas the formula (4) uses generic lags $\tau_1, \cdots, \tau_n$.

You don’t have to multiply by $e^{-i \pi \alpha t}$ before averaging the lag product. You can either multiply the lag product by $e^{-i 2 \pi \alpha t}$ or multiply $s(t+\tau/2)$ by $e^{-i \pi \alpha t}$ and $s^*(t-\tau/2)$ by $e^{-i\pi\alpha t}$ and then multiply them together.

I suggest that you first try to get the cyclic autocorrelation correct for a single cycle frequency before you try to write code that computes the entire cyclic autocorrelation surface. I have seen many attempts at the latter that fail due to not understanding how to do the former. By the way, the surface I show in the post includes cyclic autocorrelation estimates for only those cycle frequencies that are known to exist for the rectangular-pulse BPSK signal I used.

6. msaserg says:

You wrote in your blog https://cyclostationary.blog/2015/09/28/the-cyclic-autocorrelation-for-rectangular-pulse-bpsk/
“The non-conjugate CAF can be interpreted as the correlation between the lag product and the complex-valued sine wave ”.
Suppose I have a signal x1…x10.
I want to calculate the CAF for tau=2 and alf=0.1.
Then the calculations should be like this?
1. Lag product of the signal X for tau=2
x1 x2 x3 x4 x5 x6 x7 x8 .* conj(x3 x4 x5 x6 x7 x8 x9 x10)=
x1*conj(x3) x2* conj(x4) x3* conj(x5) x4* conj(x6) x5* conj(x7) x6* conj(x8) x7* conj(x9) x8* conj(x10)
2. Complex-valued sine wave
exp(2*pi*alf*t) t=0…7
e0 e1 e2 e3 e4 e5 e6 e7
3. Correlation between the lag product and the complex-valued sine wave
R(alf,tau)=R(0.1,2)=
= x1* conj(x3)*e0+x2* conj(x4)*e1+x3* conj(x5)*e2+x4* conj(x6)*e3+x5* conj(x7)*e4+x6* conj(x8)*e5+x7* conj(x9)*e6+ x8* conj(x10)*e7

1. Almost correct; you definitely have the right idea.

In Step 2, you’ll want a complex-valued exponential: $e^{-i2\pi\alpha t}$. (Don’t forget the minus sign and $i$.)

In Step 3, you’ll want to divide the sum by the number of elements that you have included. In this case $8$.

I think you chose a signal having only ten samples to make the comment formatting easy, but just to be sure I need to say that you’ll need to average over longer durations for any signal of practical interest.

1. msaserg says:

How averaging should be carried out?
Should one just take X long enough, for example, 1000
or, at the first stage, calculate the lag product for x1-x10, x11-x20, …, x991-x1000 and then average, getting the one vector of the lag product of 8 elements long?

1. The former is correct. Make sure your signal vector is long enough (pop-quiz: how?) and then use the method you wrote in your earlier comment with my slight modifications.

Each value of the lag product must be multiplied by a corresponding element of the complex sine wave to achieve the correlation between the lag product and the sine wave. If you take your long X vector and create a lag product and then average that, you’ll be destroying any sine-wave components (except for $\alpha = 0$) that might exist in the lag product.

1. msaserg says:

Well, but I do not understand then what does your phrase mean in this blog:
“You don’t have to multiply by exp() BEFORE AVERAGING the lag product”.
And in your blog https://cyclostationary.blog/2015/09/28/the-cyclic-autocorrelation/ is equations 4 and 7. In equation 7, multiply by exponent is actually performed after averaging (E[…]).
Then how is it that the former is correct (just take X long)?

pop-quiz: how? – Determined by the required resolution Frequency and Cycle Frequency.

2. Where did you find the quote you made:

“You don’t have to multiply by exp() BEFORE AVERAGING the lag product”.

I mean, where did I write that?

And in your blog https://cyclostationary.blog/2015/09/28/the-cyclic-autocorrelation/ is equations 4 and 7. In equation 7, multiply by exponent is actually performed after averaging (E[…]).

Yes, if you could average over the ensemble using $E[\cdot]$, then you could extract the cyclic autocorrelation from the resulting time-varying function through the Fourier-series operation.

But your original question was about data $x_1, x_2, \ldots$, not an ensemble. Don’t forget Equation (8) in that post! That’s what you use on data.

Does that help?

1. Thanks Serg! I see the problem now. The full quote is:

You don’t have to multiply by $e^{-i \pi \alpha t}$ before averaging the lag product. You can either multiply the lag product by $e^{-i 2 \pi \alpha t}$ or multiply $s(t+\tau/2)$ by $e^{-i \pi \alpha t}$ and $s^*(t-\tau/2)$ by $e^{-i\pi\alpha t}$ and then multiply them together.

And the context was some code written by readers in which the input signal was multiplied by $e^{-i\pi\alpha t}$ prior to forming the lag product. So what I said was a rather trivial observation that you can multiply the signal by an exponential related to $\alpha$ (and also multiply the conjugate by the appropriate exponential!), then form a lag product, then average. Or you can form the lag product, then multiply it by $e^{-i2\pi\alpha t}$.

7. msaserg says:

I realized simple model of the implementation of CAF and SCF.
The result is generally similar to that in your blog.
But there are a few problems.
1.For a single input (y_of_t = x_of_t1 + 0.0*x_of_t2+ n_of_t;), non-conjugate CAF and conjugate CAF match your results.
But if there are two input signals (y_of_t = x_of_t1 + 0.99*x_of_t2+ n_of_t;), then two signals are visible on the conjugate CAF, and only one on the conjugate CAF, that is, one signal is lost. Why?
2. Your non-conjugate SCF has a look that is different from my result (lagdot=x1.*conj(x2);).
And on the graph of the conjugate SCF (lagdot=x1.*(x2);), the cyclic frequency axis is digitized in the other direction (if without fliplr, fliplr(axis_alf)).

% CAF and CSP
clc

%% Initialize key variables.

T_bit = 10; % 1/T_bit is the bit rate
num_bits = 4000; % Desired number of bits in generated signal
fc1 = 0.05; % Desired carrier frequency(normalized units)
fc2=0.055;
N0_dB = 5.0; % Noise spectral density(average noise power)

%% Create the baseband signal.

% Create bit sequence.

bit_seq1 = randi([0 1], [1 num_bits]);
bit_seq2 = randi([0 1], [1 num_bits]);

% Create symbol sequence from bit sequence.

sym_seq1 = 2*bit_seq1 – 1;
zero_mat1 = zeros((T_bit – 1), num_bits);
sym_seq1 = [sym_seq1 ; zero_mat1];
sym_seq1 = reshape(sym_seq1, 1, T_bit*num_bits);

sym_seq2 = 2*bit_seq2 – 1;
zero_mat2 = zeros((T_bit – 1), num_bits);
sym_seq2 = [sym_seq2 ; zero_mat2];
sym_seq2 = reshape(sym_seq2, 1, T_bit*num_bits);

% Create pulse function.

p_of_t = ones(1, T_bit);

% Convolve bit sequence with pulse function.

s_of_t1 = filter(p_of_t, [1], sym_seq1);
s_of_t2 = filter(p_of_t, [1], sym_seq2);

%% Frequency-shift the baseband signal and add noise.

% Apply the carrier frequency.

e_vec1 = exp(sqrt(-1)*2*pi*fc1*[1:length(s_of_t1)]);
e_vec2 = exp(sqrt(-1)*2*pi*fc2*[1:length(s_of_t2)]);
x_of_t1 = s_of_t1 .* e_vec1;
x_of_t2 = s_of_t2 .* e_vec2;

n_of_t = randn(size(x_of_t1)) + sqrt(-1)*randn(size(x_of_t1));
noise_power = var(n_of_t);
N0_linear = 10^(N0_dB/10);
pow_factor = sqrt(N0_linear / noise_power);
n_of_t = n_of_t * pow_factor;
y_of_t = x_of_t1 + 0.99*x_of_t2+ n_of_t;%!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x=y_of_t;

Fs=1000000;
tau_max_t=15*10^-6;
tau_max_s=round(tau_max_t*Fs);

CAF=[];
for alf=0:1000:600000
alf

alfn=alf/Fs;

CAF1=[];

t=(1:length(x)-2*tau_max_s);
for tau=-tau_max_s:tau_max_s
x1=x(t+tau_max_s);
x2=x(t+tau+tau_max_s);

lagdot=x1.*conj(x2); % non-conjugate CAF
% lagdot=x1.*(x2); % conjugate CAF

ex=exp(-sqrt(-1)*2*pi*alfn*(0:length(lagdot)-1));

alftau=abs(sum(lagdot.*(ex))/length(lagdot));

CAF1=[CAF1 alftau];

end;
CAF=[CAF; CAF1]; % Cyclic Autocorrelation
end;

SCF=abs(fft(CAF’));
SCF=(fftshift(SCF)); % Spectral Correlation Function, Cyclic Wiener Relationship

axis_alf=[0:1000:600000];
axis_tau=[-tau_max_s:tau_max_s]/Fs;
bin=(Fs/length([-tau_max_s:tau_max_s]));
axis_freq=bin*[-tau_max_s:tau_max_s];

figure(1)
mesh(axis_tau, axis_alf, CAF)

figure(2)
mesh(fliplr(axis_alf), axis_freq, SCF)

1. Regarding Question 1, which is “Why are there two signals visible in the conjugate CAF but only a single signal in the non-conjugate CAF,” the reason is that you have generated two signals with the same bit rate but different carrier frequencies. Recall that for these BPSK signals, the non-conjugate cycle frequencies are harmonics of the bit rate. So here that means harmonics of $1/10$, or $\alpha = k/10$. For example, for $k = 1$, each signal contributes its CAF for $\alpha = 0.1$ to the CAF for the data. For any two statistically independent signals $s_1(t)$ and $s_2(t)$ the sum signal $x(t) = s_1(t) + s_2(t)$ has a CAF or SCF that is the sum of the two individual CAFs or SCFs:

$\displaystyle S_x^\alpha(f) = S_{s_1}^\alpha(f) + S_{s_2}^\alpha(f).$

When only one signal possesses a particular cycle frequency, say $\alpha_1$ is possessed by $s_1(t)$, then

$\displaystyle S_x^{\alpha_1}(f) = S_{s_1}^{\alpha_1}(f).$

But if they both possess a cycle frequency, say $\alpha_0$, then the SCF of the sum signal is just the sum of the two SCFs (which may add constructively or destructively because they are complex-valued and depend on things like symbol-clock phase):

$\displaystyle S_x^{\alpha_0}(f) = S_{s_1}^{\alpha_0}(f) + S_{s_2}^{\alpha_0}(f).$

Using the code you posted, I blindly estimated the non-conjugate and conjugate CFs for the sum signal using the SSCA:

Then I modified the code slightly so that the second BPSK signal has a different bit rate than the first. The first one still has a bit rate of $1/10$, but the second one has rate $1/11$:

(I used $65536$ samples, so I had to increase the number of bits from $4000$.)

Regarding Question 2, I’m not sure what you mean by the reverse digitization of the cycle-frequency axis, but in most of my surface plots I tend to view the surface so that cycle frequency increases toward the viewer and spectral frequency increases to the left of the viewer. This is just for convenience and is a consequence of the fact that in the non-conjugate surface, the PSD is the largest slice and tends to block the view of the $\alpha \neq 0$ slices.

1. msaserg says:

Thank you!
I apologize for the not entirely specific second question, but apart from digitizing the axis, there is a difference in the form of a SCF. You have (https://cyclostationary.files.wordpress.com/2015/02/ww_1.jpg) non-conjugate SCF is different from my non-conjugate SCF (alpha = 0.3 f = 0.1, and mine f = 0)

1. As you compute your estimates of the cyclic autocorrelation, you have decided to take their magnitude and store in a matrix:
 alftau=abs(sum(lagdot.*(ex))/length(lagdot)); CAF1=[CAF1 alftau]; 
Later you take the Fourier transform of the obtained cyclic autocorrelation estimates, but remember that you are taking the Fourier transform of the magnitude (abs) of the estimates here. The phase information in the cyclic autocorrelation estimates is what leads to the frequency shift that you’re missing.

You’ll also need to be careful about transpositions once you restore the complex values:
 SCF=abs(fft(CAF')); 
And you might have a sign error in how you implement the delay:
 for tau=-tau_max_s:tau_max_s x1=x(t+tau_max_s); x2=x(t+tau+tau_max_s); 

8. Claudio Dias says:

Dear Prof. Spooner,

There is a typo in the begining of the post. There is “$100 + k 100$” text instead of latex rendering.

9. John Snoap says:

Hi Dr. Spooner,

Been working on implementing a cyclic-autocorrelation function in MATLAB based on what I’ve read so far on your blog/papers!

While I believe the implementation is close because when alpha=0 it looks very similar to the plots you have in this article, when alpha is not 0 my plots look the same as when alpha=0. I’m assuming I’m either not implementing something correctly or I’m misunderstanding something about cyclic-autocorrelation. Can I email you the details of what I have so far to get help in learning what I’m misunderstanding?

Thank you,
John Snoap