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