%% Method 1: Decimation of Autocorrelation Function
% DetectionWithDecimation calls OFDM to simulate the cyclostationary
% feature detector with decimation.
% A simulated OFDM signal is multiplied with a lagged version of itself,
% where lag is l bits. The resultant autocorrelation is decimated by a
% factor of m and the Nifft-point DFT of the result is obtained. The f_bin
% frequency bin of the fft is used as a means of comparison, R.
% For various SNR and M values, the variance of R is measured when no signal
% is present. The mean of R is assumed to be 0. var_trials are used to measure
% the variance.
% The variance is then used to obtain the chi-square test statistic, which is then
% compared with thresholds to form a signal presence hypothesis. detection_trials
% simulations of R, both when signal is and isn't present, are used to measure
% the probabilities of detection and false positives.
% The chi-square test statistic threshold is the quantile function of probability, 1-p,
% for 2 degrees of freedom. p is the chosen probability of false positives.
SNR = -20:1:5; % signal of interest to noise ratio (in dB)
fs_y = 500*10^6; % sampling frequency of generated OFDM signals
fs_ym = 500*10^6; % sampling frequency of detector
Nifft = 2048; % number of samples for IFFT
M = [1 2 4 8 16]; % decimation factor
var_trials = 150; % total number of trials for computing variance
detection_trials = 150; % total number of trials for detection
l = 500; % lag (in number of samples)
p = [.20 .1 .05 .01]; % probabilities of false alarm
f_bin = round(10*fs_ym/fs_y); % frequency bin of interest
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Measure Mean and Variance of Cyclic Autocorrelation Function %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
data_point = 0;
h = waitbar(0,'Measuring noise statistics...');
for snr = SNR
% Generate Noise
y = [];
t = [];
for k = 1:ceil((max(M)*Nifft*var_trials+l)/(1.25*10^-3*fs_ym))
[y_sample t_sample theta] = OFDM(0,snr,fs_ym);
y = [y;y_sample];
end
% Compute the autocorrelation of y for lag, l bits
y_lagged = [zeros(l,1);y(1:length(y)-l)];
ry = y.*conj(y_lagged);
for m = M
data_point = data_point + 1;
%Decimate autocorrelation
ycorr = decimate(ry,m);
for trial = 1:var_trials
% Compute cyclic autocorrelation
Ycorr = fft(ycorr(((trial-1)*Nifft+1):trial*Nifft))/Nifft;
% Record cyclic autocorrelations at frequency bin of interest
R(trial) = Ycorr(f_bin*m);
end
% Compute variances
var_real(data_point) = mean(real(R).^2);
var_im(data_point) = mean(imag(R).^2);
var_realim(data_point) = mean(imag(R).*real(R));
end
waitbar(find(snr==SNR)/length(SNR),h,'Measuring noise statistics');
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% compute thresholds for R %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
thresholds = chi2inv(1-p,2);
%%%%%%%%%%%%%%%%%%
% test detection %
%%%%%%%%%%%%%%%%%%
waitbar(0,h,'Testing detection...');
data_point = 0;
for snr = SNR
% Generate Noise and OFDM Signal with Noise
y_present = [];
y_npresent = [];
for k = 1:ceil((max(M)*Nifft*detection_trials+l)/(1.25*10^-3*fs_ym))
% Generate signal with data
[y_sample t_sample theta] = OFDM(1,snr,fs_ym);
y_present = [y_present;y_sample];
% Generate signal without data
[y_sample t_sample theta] = OFDM(0,snr,fs_ym);
y_npresent = [y_npresent;y_sample];
end
% Compute the autocorrelation of y for a lag of l samples
y_lagged_present = [zeros(l,1);y_present(1:length(y_present)-l)];
ry_present = y_present.*conj(y_lagged_present);
y_lagged_npresent = [zeros(l,1);y_npresent(1:length(y_npresent)-l)];
ry_npresent = y_npresent.*conj(y_lagged_npresent);
for m = M
data_point = data_point + 1;
% generate variance matrix
sigma = [var_real(data_point),var_realim(data_point);var_realim(data_point),var_im(data_point)];
% Decimate autocorrelations
ycorr_present = decimate(ry_present,m);
ycorr_npresent = decimate(ry_npresent,m);
for trial = 1:detection_trials
% Compute cyclic autocorrelations and record them at the frequency bin of interest
Ycorr_present = fft(ycorr_present(((trial-1)*Nifft+1):trial*Nifft))/Nifft;
Ycorr_npresent = fft(ycorr_npresent(((trial-1)*Nifft+1):trial*Nifft))/Nifft;
R_present = Ycorr_present(f_bin*m);
R_npresent = Ycorr_npresent(f_bin*m);
R_present = [real(R_present),imag(R_present)];
R_npresent = [real(R_npresent),imag(R_npresent)];
% Track Detection
detected(trial,:) = R_present*inverse(sigma)*R_present'>thresholds;
% Track False Positives
falsepositive(trial,:) = R_npresent*inverse(sigma)*R_npresent'>thresholds;
end
% Compute Probability of Detection
pdetection(data_point,:) = mean(detected);
% Compute Probability of False Positives (their averages are approximately the expected values)
pfalsepositive(data_point,:) = mean(falsepositive);
end
waitbar(find(snr==SNR)/length(SNR),h,'Testing detection...');
end
pdetection1 = reshape(pdetection(:,1),[length(M),length(SNR)]);
pdetection2 = reshape(pdetection(:,2),[length(M),length(SNR)]);
pdetection3 = reshape(pdetection(:,3),[length(M),length(SNR)]);
pdetection4 = reshape(pdetection(:,4),[length(M),length(SNR)]);
pfalsepositive1 = reshape(pfalsepositive(:,1),[length(M),length(SNR)]);
pfalsepositive2 = reshape(pfalsepositive(:,2),[length(M),length(SNR)]);
pfalsepositive3 = reshape(pfalsepositive(:,3),[length(M),length(SNR)]);
pfalsepositive4 = reshape(pfalsepositive(:,4),[length(M),length(SNR)]);
close(h);