% clear all variables
clear;

sig_num = 300;
fd = 1;
fs = 7;
N = 64;
fm = 10;  %unit: 1Hz
Rs = 1;  %unit: 1MHz
Ts = 1/Rs;

% generate random QPSK signal
qpsk_D=randint(sig_num,1,4);
qpsk_C=dmodce(qpsk_D,fd,fs,'psk',4);
qpsk_flt = rcosflt(qpsk_C,fd,fs);

% draw eye diagram and scatter plot of QPSK signal
qpsk=qpsk_flt(1+3*fs:end-3*fs);
%eyediagram(qpsk,2*fs,2*fs);
%scatterplot(qpsk,fs);

% QPSK signal passes through AWGN channel
snr=8;
qpsk_AWGN_C=awgn(qpsk_C,snr);
qpsk_AWGN_flt=rcosflt(qpsk_AWGN_C,fd,fs);
qpsk_AWGN=qpsk_AWGN_flt(1+3*fs:end-3*fs);
%eyediagram(qpsk_AWGN,2*fs,2*fs,0,'r-');
%scatterplot(qpsk_AWGN,fs,0,'r.');

% generate Doppler power spectrum
deltaf = 2*fm/(N-1);
T = 1/deltaf;
sf0 = 1.5/(pi*fm);
for n = 1:(N-2)/2
    sf(n) = 1.5/(pi*fm*sqrt(1-(n*deltaf/fm)^2));
end
classicf = [fliplr(sf),sf0,sf];
%figure;
%plot(classicf);

for channels = 1:3
    % generate two normally distributed random variables
	gaussN_re1 = randn(1,(N-2)/2);
    gaussN_im1 = randn(1,(N-2)/2);
    gaussN_pos1 = gaussN_re1 + i*gaussN_im1;
	gaussN_neg1 = conj(gaussN_pos1);
	%% gaussN1 = [fliplr(gaussN_pos1),0,gaussN_pos1];
	gaussN1 = [fliplr(gaussN_neg1),0,gaussN_pos1];
	
	%% gaussN_pos2 = randn(1,(N-2)/2);
	%% gaussN2 = [fliplr(gaussN_pos2),0,gaussN_pos2];
	gaussN_re2 = randn(1,(N-2)/2);
    gaussN_im2 = randn(1,(N-2)/2);
    gaussN_pos2 = gaussN_re2 + i*gaussN_im2;
	gaussN_neg2 = conj(gaussN_pos2);
	gaussN2 = [fliplr(gaussN_neg2),0,gaussN_pos2];
	
	% generating flat Rayleigh fading channel
	x = ifft(sqrt(classicf).*gaussN1);
	y = ifft(sqrt(classicf).*gaussN2);
	%rayleigh_amp = sqrt(abs(x).^2+abs(y).^2);
	%rayleigh_db = 10*log10(rayleigh_amp);
	%figure;
	%plot(rayleigh_db);
	
	rayleighOrg = x+i*y;
	
	% frequency matching
	repNum = round(T*1e6*fs/Ts);
	dotNeed = floor(sig_num*fs/repNum);
	for k = 1:dotNeed
	    for j = 1:repNum
	        rayleighCh(j,k) = rayleighOrg(k);
	    end
	end
	remain = sig_num*fs-dotNeed*repNum;
	if remain > 0
	    for j=1:remain
	        rayleighCh(j,dotNeed+1) = rayleighOrg(dotNeed+1);
	    end
	end
	
    if channels == 1
	    rayleighCh1 = reshape(rayleighCh,1,[]);
    elseif channels == 2
	    rayleighCh2 = reshape(rayleighCh,1,[]);
    else
 	    rayleighCh3 = reshape(rayleighCh,1,[]);
    end

    %clear temporary variants for next use
    clear gaussN_re1 gaussN_im1 gaussN_re2 gaussN_im2 gaussN_neg1 gaussN_neg2;
    clear gaussN_pos1 gaussN_pos2 gaussN1 gaussN2;
    clear x y rayleighOrg rayleighCh;
end

% generate 3-path Rayleigh channel
a0 = 0.476;
a1 = 0.048;
a2 = 0.476;

% add delay time of Rayleigh channel
delay1 = 1;     %unit: 1us
delay2 = 2;
r1 = round(delay1/Ts*fs);
r2 = round(delay2/Ts*fs);

rayleigh0 = a0*rayleighCh1(1:fs*sig_num);
rayleigh1 = a1*rayleighCh2(1:fs*sig_num);
rayleigh2 = a2*rayleighCh3(1:fs*sig_num);

signal = reshape(qpsk_C,1,fs*sig_num);
sig_delay1 = [zeros(1,r1),signal(1:fs*sig_num-r1)];
sig_delay2 = [zeros(1,r2),signal(1:fs*sig_num-r2)];

sig_ray0 = signal.*rayleigh0;
sig_ray1 = sig_delay1.*rayleigh1;
sig_ray2 = sig_delay2.*rayleigh2;

% output of 3-path rayleigh channels
sig_ray = sig_ray0+sig_ray1+sig_ray2;

sig_flt= rcosflt(sig_ray,fd,fs);
out=sig_flt(1+3*fs:end-3*fs);

%% out1 = out(1:fs*sig_num);

eyediagram(out,2*fs,2*fs,0,'m-');
scatterplot(out,1,0,'m.');