%a_run_design.m文件
% setup				%系統(tǒng)設(shè)置
% QAM				%QAM調(diào)制方式
% OFDM				%OFDM調(diào)制方式
% Analysis				%QAM和OFDM兩種方式對比

% 1setup            %初始化部分
% setup
disp(' '), disp('------------------------------------------------------------')
disp('Simulation Setup')

% OFDM Setup -----------------------------------------------------------
fft_size = 128		% should be a power of 2 for fast computation
				% more points = more time domain samples (smoother & more cycles)
num_carriers = 32	% should be <= fft_size/4
				% number of carriers used for each data chunk
% new var - denotes even spacing or variations of carriers among fft points
input_type = 2;
% 1 = test input
	test_input_type = 1;
	% 1 = bit specified (binary)
		binary_data = [0 1 0 1 0 1 0 1];
	% 2 = random data stream (samples in the range of 0-255)
		num_symbols = 9;
	% 3 = sinusoidal
		frequency = 2;
		num_samples = 50;
% 2 = external file input
	file_name = 'shortest.wav';		% Name of input file
	file_input_type = 3;
		% 1 = binary (not implemented)
		% 2 = text					% Demo file:  'text.txt'
		% 3 = sound				% Demo files: 'shortest.wav' & 'shorter.wav'
		% 4 = image (not implemented)

% QAM Setup ------------------------------------------------------------
do_QAM = 1;		% (1=on, 0=off)
QAM_periods = 10;	% defines the number of periods per QAM Symbos (1=2*pi)

% Channel Simulation Parameters --------------------------------------------
channel_on = 1;		% 1=on, 0=off
clip_level = 1.0;		% 0.0 - 1.0 (0-100%)
	% Max magnitude of the signal is 'clip_level' times the full magnitude of the signal
noise_level = 0.0;	% 0.0 - 1.0 (0-100%)
	% Magnitude of noise is 'noise_level' times the magnitude of the signal
% Multipath Channel Simulation
	% Good defaults when fft_size = 128 and num_carriers = 32:
	%	d1=6; a1=0.30; d2=10; a2=0.25
	d1 = 6;		% delay in units
	a1 = 0.30; 	% attenuation factor - multipath signal is x% of size or original signal
	d2 = 10;		% delay for second multipath signal
	a2 = 0.25;	% attenuation factor for second multipath signal


% ****************** TEST INPUT SETUP - DO NOT MODIFY **************************
if input_type == 1
	if test_input_type == 1
		%specify BINARY input bit-by-bit
		data_in = binary_data;
	end
	if test_input_type == 2
		%random input defined by parameters
		num_levels = 255;		%number of possible levels of a symbol
						%must be integer between 1-255
		data_samples = round(rand(1,num_symbols)*(num_levels-1));
		data_in = zeros(1,8*length(data_samples));
		for i = 1:length(data_samples)
			data_in(1 + (i-1)*8:(i-1)*8 + 8) = eight2bin(data_samples(i));
		end
	end
	if test_input_type == 3
		%data stream represents sine wave samples
		t = linspace(0,1,num_symbols);	%evenly space number of samples
		%take 8-bit samples of sine wave
		data_samples = round(127.5*sin(frequency*2*pi*t) +127.5);
		data_in = zeros(1,8*length(data_samples));
		for i = 1:length(data_samples)
			data_in(1 + (i-1)*8:(i-1)*8 + 8) = eight2bin(data_samples(i));
		end
	end
end

already_made_noise = 0;	% initialization (don't change)

% 2、QAM.m
% QAM.m compares OFDM (multicarrier) to multi-level QAM (single carrier)
% when they transmit the same # of bits in a given time period

read		% read data for QAM - does not affect OFDM
data_in_pol = bin2pol(data_in);		% Converts binary data to polar data

% check to see if num_carriers is a power of 2
is_pow_2 = num_carriers;
temp_do_QAM = 0;
if is_pow_2 ~= 2
	while temp_do_QAM == 0
		temp_do_QAM = rem(is_pow_2,2);
		is_pow_2 = is_pow_2/2;
		if is_pow_2 == 2
			temp_do_QAM = -99;		% it is a power of 2 -> can do QAM
		end
	end
else
	temp_do_QAM = -99;	% 2 is a power of 2
end
if temp_do_QAM ~= -99
	do_QAM = 0;	% don't do it if it's not possible
	disp(' '),disp('ERROR: Cannot run QAM because num_carriers is not valid.')
	disp('       Please see "setup.m" for details.')
end

%------------------------------------------------------------------------
%                           單獨的QAM調(diào)制通過信道及解調(diào)
%------------------------------------------------------------------------
if do_QAM == 1
	tic	% Start stopwatch to calculate how long QAM simulation takes
	
	disp(' '), disp('------------------------------------------------------------')
	disp('QAM simulation'), disp('Transmitting')
	
	% Pad with zeros so data can be divided evenly
	data_length = length(data_in_pol);
	r = rem(data_length,num_carriers);
	if r ~= 0
		for i = 1:num_carriers-r
			data_in_pol(data_length+i) = 0;	%pad input with zeros to complete last data set
		end							%speed improve possible
	end
	data_length = length(data_in_pol);		%update after padding
	
	num_OFDM_symbols = ceil(data_length / (2*num_carriers));
	% num QAM symbols that represent equal amount of data to one OFDM symbol
	num_QAM_symbols = num_carriers / 2;
	% num samples per QAM symbol
	num_symbol_samples = fft_size / num_QAM_symbols;
	
	% convert polar data [-1, 1] to 4 level data [-3, -1, 1, 3]
	data_in_4 = zeros(1,data_length/2);
	for i = 1:2:data_length
		data_in_4(i - (i-1)/2) = data_in_pol(i)*2 + data_in_pol(i+1);
	end
	
	% define sample points between 0 and 2*pi
	ts = linspace(0, 2*pi*QAM_periods, num_symbol_samples+1);
	
	% Generate 16-QAM data
	% total length of 16-QAM transmission
	tx_length = num_OFDM_symbols * num_QAM_symbols * num_symbol_samples;
	QAM_tx_data = zeros(1,tx_length);
	for i = 1:2:data_length/2
		for k = 1:num_symbol_samples
			QAM_tx_data(k+((i-1)/2)*num_symbol_samples) = data_in_4(i)*cos(ts(k)) + data_in_4(i+1)*sin(ts(k));
		end
	end
	
	% Do channel simulation on QAM data
	xmit = QAM_tx_data;		% ch uses 'xmit' data and returns 'recv'
	ch
	QAM_rx_data = recv;		% save QAM data after channel
	clear recv				% remove 'recv' so it won't interfere with OFDM
	clear xmit				% remove 'xmit' so it won't interfere with OFDM
	
	disp('Receiving')			% Recover Binary data (Decode QAM)
	cos_temp = zeros(1,num_symbol_samples);		%
	sin_temp = cos_temp;					%
	xxx = zeros(1,data_length/4);				% Initialize to zeros for speed
	yyy = xxx;								%
	QAM_data_out_4 = zeros(1,data_length/2);		%
	
	for i = 1:2:data_length/2	% "cheating"
		for k = 1:num_symbol_samples
			% multiply by carriers to produce high frequency term and original data
			cos_temp(k) = QAM_rx_data(k+((i-1)/2)*num_symbol_samples) * cos(ts(k));
			sin_temp(k) = QAM_rx_data(k+((i-1)/2)*num_symbol_samples) * sin(ts(k));
		end
		% LPF and decide - we will do very simple LPF by averaging
		xxx(1+(i-1)/2) = mean(cos_temp);
		yyy(1+(i-1)/2) = mean(sin_temp);
		% Reconstruct data in serial form
		QAM_data_out_4(i) = xxx(1+(i-1)/2);
		QAM_data_out_4(i+1) = yyy(1+(i-1)/2);
	end
	
	% Make decision between [-3, -1, 1, 3]
	for i = 1:data_length/2
		if QAM_data_out_4(i) >= 1, QAM_data_out_4(i) = 3;
		elseif QAM_data_out_4(i) >= 0, QAM_data_out_4(i) = 1;
		elseif QAM_data_out_4(i) >= -1, QAM_data_out_4(i) = -1;
		else QAM_data_out_4(i) = -3;
		end
	end
	
	% Convert 4 level data [-3, -1, 1, 3] back to polar data [-1, 1]
	QAM_data_out_pol = zeros(1,data_length);	% "cheating"
	for i = 1:2:data_length
		switch QAM_data_out_4(1 + (i-1)/2)
			case -3
				QAM_data_out_pol(i) = -1;
				QAM_data_out_pol(i+1) = -1;
			case -1
				QAM_data_out_pol(i) = -1;
				QAM_data_out_pol(i+1) = 1;
			case 1
				QAM_data_out_pol(i) = 1;
				QAM_data_out_pol(i+1) = -1;
			case 3
				QAM_data_out_pol(i) = 1;
				QAM_data_out_pol(i+1) = 1;
			otherwise
				disp('Error detected in switch statment - This should not be happening.');
		end
	end
	QAM_data_out = pol2bin(QAM_data_out_pol);	% convert back to binary
	
	% Stop stopwatch to calculate how long QAM simulation takes
	QAM_simulation_time = toc;
	if QAM_simulation_time > 60
		disp(strcat('Time for QAM simulation=', num2str(QAM_simulation_time/60), ' minutes.'));
	else
		disp(strcat('Time for QAM simulation=', num2str(QAM_simulation_time), ' seconds.'));
	end
end
%-------------------------------------------------------------------------
%                             OFDM仿真
%-------------------------------------------------------------------------
% 3 % Run OFDM simulation
tic	% Start stopwatch to calculate how long QAM simulation takes
disp(' '),disp('------------------------------------------------------------')
disp('OFDM Simulation')
tx         
ch
rx
% Stop stopwatch to calculate how long QAM simulation takes
OFDM_simulation_time = toc;
if OFDM_simulation_time > 60
	disp(strcat('Time for OFDM simulation=', num2str(OFDM_simulation_time/60), ' minutes.'));
else
	disp(strcat('Time for OFDM simulation=', num2str(OFDM_simulation_time), ' seconds.'));
end

% 3.1發(fā)送
% tx
disp('Transmitting')
read                            %read original data
data_in_pol = bin2pol(data_in);	   % Converts binary data to polar data
tx_chunk                        %convert polar data into chunks.      

% perform ifft to create time domain waveform representing data
td_sets = zeros(num_chunks,fft_size);
for i = 1:num_chunks
	td_sets(i,1:fft_size) = real(ifft(spaced_chunks(i,1:fft_size)));
end

tx_dechunk      % Construct signal to transmit by placing time domain sets in series

% 3.1.1 % tx_chunk    %雙極性數(shù)組轉(zhuǎn)化成OFDM字符串
data_length = length(data_in_pol)			    	%number of symbols in original input
num_carriers
num_chunks = ceil(data_length/(2*num_carriers))	%2 data on each carrier (real and imaginary)
r = rem(data_length,2*num_carriers)
if r ~= 0
	for i = 1:num_carriers*2-r
		data_in_pol(data_length+i) = 0;	%pad input with zeros to complete last data set
	end							%speed improve possible
end

% break data into chunks
chunks = zeros(num_chunks,num_carriers);	% for speed
for i = 1:num_chunks
	% *********************chunk done
	for k = 1:num_carriers
		chunks(i,k) = data_in_pol(2*num_carriers*(i-1)+k) + data_in_pol(2*num_carriers*(i-1)+k+num_carriers)*j;
	end
end

chunks
% Padding chunks with zeros so num_carriers and fft_size are compatible
% Once compatible, further spacing is simplified
num_desired_carriers = num_carriers;
num_zeros = 0;
thinking = 1;
while thinking == 1 % Continue if num_carriers and fft_size are not compatible
	if rem(fft_size/2,num_desired_carriers) == 0
		thinking = 0;
	else
		num_desired_carriers = num_desired_carriers + 1;
		num_zeros = num_zeros + 1;
	end
end

padded_chunks = zeros(num_chunks,num_carriers + num_zeros); % for speed
padded_chunks(1:num_chunks,num_zeros + 1:num_carriers + num_zeros) = chunks;


%compute zeros_between
zeros_between = ((fft_size/2) - (num_carriers + num_zeros))/(num_carriers + num_zeros);

spaced_chunks = zeros(num_chunks,fft_size); % for speed - extra room for folding later
%add zeros_between
i = 1;
for k = zeros_between +1:zeros_between +1:fft_size/2
	spaced_chunks(1:num_chunks,k) = padded_chunks(1:num_chunks,i);
	i = i+1;
end

% folding data to produce an odd function for ifft input
for i = 1:num_chunks
	% Note: index = 1 is actually DC freq for ifft -> it does not get copied over y-axis
	spaced_chunks(i,fft_size:-1:fft_size/2+2) = conj(spaced_chunks(i,2:fft_size/2));
end

% 3.1.2 tx_dechunk
% tx_dechunk

% Construct signal to transmit by placing time domain sets in series
xmit = zeros(1,num_chunks*fft_size);
for i = 1:num_chunks
	for k = 1:fft_size
		xmit(k + (i-1)*fft_size) = td_sets(i,k);
	end
end

% 3.2 %ch.m
% ch
recv = xmit;	% channel is applied to recv, don't modify transmitted data
if channel_on == 1
	disp('Simulating Channel')
	norm_factor = max(abs(recv));	% Normalize all data before applying
	recv = (1/norm_factor) * recv;	% channel for a fair comparison
	ch_clipping                     %clipp data
	ch_multipath                    %
	ch_noise                        %
	recv = norm_factor * recv;		% Restore data magnitude for proper decoding
end

% 3.2.1 %ch_clipping.m
% ch_clipping
for i = 1:length(recv)
	if recv(i) > clip_level
		recv(i) = clip_level;
	end
	if recv(i) < -clip_level
		recv(i) = -clip_level;
	end
end

% 3.2.2 % ch_multipath                  % 產(chǎn)生多經(jīng)的方法
copy1=zeros(size(recv));