for i=1+d1:length(recv)
	copy1(i)=a1*recv(i-d1);
end
copy2=zeros(size(recv));
for i=1+d2:length(recv)
	copy2(i)=a2*recv(i-d2);
end
recv=recv+copy1+copy2;

% 3.2.3、%ch_noise                     %施加信道噪聲
% ch_noise (operate on recv)
% random noise defined by noise_level amplitude
if already_made_noise == 0	% only generate once and use for both QAM and OFDM
	noise = (rand(1,length(recv))-0.5)*2*noise_level;
	already_made_noise = 1;
end
recv = recv + noise;

% 3.3 %rx.m                        %接收程序
% rx
disp('Receiving')
rx_chunk

% perform fft to recover original data from time domain sets
recv_spaced_chunks = zeros(num_chunks,fft_size);
for i = 1:num_chunks
	recv_spaced_chunks(i,1:fft_size) = fft(recv_td_sets(i,1:fft_size));
	% Note: 'round()' gets rid of small numerical error in Matlab but a threshold will be needed for a practical system
	% 2001-4-17 -- Got rid of 'round()' to do decoding more intelligently
end
rx_dechunk
output = pol2bin(output);	% Converts polar to binary
write

% 3.3.1 %rx_chunk.m
% rx_chunk
% break received signal into parellel sets for demodulation
recv_td_sets = zeros(num_chunks,fft_size);
for i = 1:num_chunks
	for k = 1:fft_size
		recv_td_sets(i,k) = recv(k + (i-1)*fft_size);
	end
end

% 3.3.2 % rx_dechunk             %并串轉換
% rx_dechunk
% take out zeros_between from recv_spaced_chunks --> recv_padded_chunks
recv_padded_chunks = zeros(num_chunks, num_carriers+num_zeros);
i = 1;
for k = zeros_between +1:zeros_between +1:fft_size/2
	recv_padded_chunks(1:num_chunks,i) = recv_spaced_chunks(1:num_chunks,k);
	i = i+1;
end

% take out num_zeros from padded chunks --> recv_chunks
recv_chunks = zeros(num_chunks, num_carriers);
recv_chunks = recv_padded_chunks(1:num_chunks, num_zeros+1:num_carriers+num_zeros);

% Recover bit stream by placing reconstructed frequency domain data in series
recv_dechunked = zeros(1, num_chunks*num_carriers);
for i = 1:num_chunks
	for k = 1:num_carriers
		recv_dechunked(k + (i-1)*num_carriers*2) = real(recv_chunks(i,k));
		recv_dechunked(k + (i-1)*num_carriers*2 + num_carriers) = imag(recv_chunks(i,k));
	end
end

% take out trailing zeros from output --> output
output_analog = recv_dechunked(1:data_length);
output = sign(output_analog);

% 3.3.3 %write               %save received data
% write
% ******************TEST OUTPUT*********************************
if input_type == 1
	if test_input_type == 1
		%already binary - do nothing
	end
	
	if (test_input_type == 2) | (test_input_type == 3)
		%random input     OR     sine wave samples
		output_samples = zeros(1,floor(length(output)/8));  %extra zeros are not original data
		for i = 1:length(output_samples)
			output_samples(i) = bin2eight(output(1 + (i-1)*8:(i-1)*8 + 8));
		end
		if do_QAM == 1
			QAM_output_samples = zeros(1,floor(length(QAM_data_out)/8));
			for i = 1:length(QAM_output_samples)
				QAM_output_samples(i) = bin2eight(QAM_data_out(1 + (i-1)*8:(i-1)*8 + 8));
			end
		end
	end
end

% ******************FILE OUTPUT*********************************
if input_type == 2

	if file_input_type == 1
		%binary file output - not implemented
	end

	if file_input_type == 2
		%text file output
		output_samples = zeros(1,floor(length(output)/8));  %extra zeros are not original data
		for i = 1:length(output_samples)
			output_samples(i) = bin2eight(output(1 + (i-1)*8:(i-1)*8 + 8));
		end
		file = fopen('OFDM_text_out.txt','wt+');
		fwrite(file,output_samples,'char');
		fclose(file);

		if do_QAM == 1
			%extra zeros are not original data
			QAM_output_samples = zeros(1,floor(length(QAM_data_out)/8));  
			
for i = 1:length(QAM_output_samples)
				QAM_output_samples(i) = bin2eight(QAM_data_out(1 + (i-1)*8:(i-1)*8 + 8));
			end
			file = fopen('QAM_text_out.txt','wt+');
			fwrite(file,QAM_output_samples,'char');
			fclose(file);
		end
	end

	if file_input_type == 3
		output_samples_big = zeros(1,floor(length(output)/8));  %extra zeros are not original data
		for i = 1:length(output_samples_big)
			output_samples_big(i) = bin2eight(output(1 + (i-1)*8:(i-1)*8 + 8));
		end
		%convert dynamic range from 0:255 to -1:1
		output_samples = (output_samples_big-127)/128;
		%sound file output
		wavwrite(output_samples, 11025, 8, 'OFDM_out.wav')
		if do_QAM == 1
			QAM_data_out_big = zeros(1,floor(length(QAM_data_out)/8));
			for i = 1:length(QAM_data_out_big)
				QAM_data_out_big(i) = bin2eight(QAM_data_out(1 + (i-1)*8:(i-1)*8 + 8));
			end
			%convert dynamic range from 0:255 to -1:1
			QAM_output_samples = (QAM_data_out_big-127)/128;
			%sound file output
			wavwrite(QAM_output_samples, 11025, 8, 'QAM_out.wav')
		end
	end
	
	if file_input_type == 4
		%image file output - not implemented		
	end
end

% 4、%Analysis.m
% Analysis
disp(' '), disp('------------------------------------------------------------')
disp('Preparing Analysis')
figure(1), clf
if (input_type == 1) & (test_input_type == 1)
	subplot(221), stem(data_in), title('OFDM Binary Input Data');
	subplot(223), stem(output), title('OFDM Recovered Binary Data')
else
	subplot(221), plot(data_samples), title('OFDM Symbol Input Data');
	subplot(223), plot(output_samples), title('OFDM Recovered Symbols');
end
subplot(222), plot(xmit), title('Transmitted OFDM');
subplot(224), plot(recv), title('Received OFDM');


% dig_x_axis = (1:length(QAM_tx_data))/length(QAM_tx_data);
% 	figure(4), clf, subplot(212)
% 	freq_data = abs(fft(QAM_rx_data));
% 	L = length(freq_data)/2;
	
dig_x_axis = (1:length(xmit))/length(xmit);
figure(2), clf

if channel_on ==1
	num = [1, zeros(1, d1-1), a1, zeros(1, d2-d1-1), a2];
	den = [1];
	[H, W] = freqz(num, den, 512);
	mag = 20*log10(abs(H));
	phase = angle(H) * 180/pi;
	
	subplot(313)
	freq_data = abs(fft(recv));
	L = length(freq_data)/2;
	plot(dig_x_axis(1:L), freq_data(1:L))
	xlabel('FFT of Received OFDM')
	axis_temp = axis;
	
	subplot(311),
	freq_data = abs(fft(xmit));
	plot(dig_x_axis(1:L), freq_data(1:L)), axis(axis_temp)
	title('FFT of Transmitted OFDM')
	
	subplot(312)
	plot(W/(2*pi),mag),
	ylabel('Channel Magnitude Response')
else
	subplot(212)	
	freq_data = abs(fft(recv));
	L = length(freq_data)/2;
	plot(dig_x_axis(1:L), freq_data(1:L))
	xlabel('FFT of Received OFDM')
	axis_temp = axis;
	
	subplot(211),
	freq_data = abs(fft(xmit));
	plot(dig_x_axis(1:L), freq_data(1:L)), axis(axis_temp)
	title('FFT of Transmitted OFDM')
end

% if file_input_type == 4
% 	figure(5)
% 	subplot(211)
% 	image(data_in);
% 	colormap(map);
% 	subplot(212)
% 	image(output);
% 	colormap(map);
% end

if do_QAM == 1 	% analyze if QAM was done
	
	figure(3), clf
	if (input_type == 1) & (test_input_type == 1)
		subplot(221), stem(data_in), title('QAM Binary Input Data');
		subplot(223), stem(QAM_data_out), title('QAM Recovered Binary Data')
	else
		subplot(221), plot(data_samples), title('QAM Symbol Input Data');
		subplot(223), plot(QAM_output_samples), title('QAM Recovered Symbols');
	end
	subplot(222), plot(QAM_tx_data), title('Transmitted QAM');
	subplot(224), plot(QAM_rx_data), title('Received QAM');
	
	dig_x_axis = (1:length(QAM_tx_data))/length(QAM_tx_data);
	figure(4), clf
	
	if channel_on ==1
		subplot(313)
		freq_data = abs(fft(QAM_rx_data));
		L = length(freq_data)/2;
		plot(dig_x_axis(1:L), freq_data(1:L))
		xlabel('FFT of Received QAM')
		axis_temp = axis;
		
		subplot(311),
		freq_data = abs(fft(QAM_tx_data));
	 	plot(dig_x_axis(1:L),freq_data(1:L)), axis(axis_temp)
		title('FFT of Transmitted QAM')
		
		subplot(312)
		plot(W/(2*pi),mag)
		ylabel('Channel Magnitude Response')
	else
		subplot(212)
		freq_data = abs(fft(QAM_rx_data));
		L = length(freq_data)/2;
		plot(dig_x_axis(1:L), freq_data(1:L))
		title('FFT of Received QAM')
		axis_temp = axis;
		
		subplot(211),
		freq_data = abs(fft(QAM_tx_data));
	 	plot(dig_x_axis(1:L),freq_data(1:L)), axis(axis_temp)
		title('FFT of Transmitted QAM')
	end
		
	% Plots the QAM Received Signal Constellation
	figure(5), clf, plot(xxx,yyy,'ro'), grid on, axis([-2.5 2.5 -2.5 2.5]), hold on

% 	% Overlay plot of transmitted constellation
% 	x_const = [-1.5 -0.5 0.5 1.5 -1.5 -0.5 0.5 1.5 -1.5 -0.5 0.5 1.5 -1.5 -0.5 0.5 1.5];
% 	y_const = [-1.5 -1.5 -1.5 -1.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5 1.5 1.5 1.5 1.5];
% 	plot(x_const, y_const, 'b*')

	% Overlay of constellation boundarys
	x1 = [-2 -2]; x2 = [-1 -1]; x3 = [0 0]; x4 = [1 1]; x5 = [2 2]; x6 = [-2 2];
	y1 = [-2 -2]; y2 = [-1 -1]; y3 = [0 0]; y4 = [1 1]; y5 = [2 2]; y6 = [-2 2];
	plot(x1,y6), plot(x2,y6), plot(x3,y6), plot(x4,y6), plot(x5,y6)
	plot(x6,y1), plot(x6,y2), plot(x6,y3), plot(x6,y4), plot(x6,y5)	

	hold off
	title('16-QAM Received Signal Constellation and Decision Boundarys')
	
	binary_err_bits_QAM = 0;
	for i = 1:length(data_in)
		err = abs(data_in(i)-QAM_data_out(i));
		if err > 0
			binary_err_bits_QAM = binary_err_bits_QAM + 1;
		end
	end
	BER_QAM = 100 * binary_err_bits_QAM/data_length;
end

figure(6), clf
if channel_on == 1
	subplot(211), plot(W/(2*pi),mag),title('Channel Magnitude Response')
	xlabel('Digital Frequency'),ylabel('Magnitude in dB')
	subplot(212), plot(W/(2*pi),phase),title('Channel Phase Response')
	xlabel('Digital Frequency'),ylabel('Phase in Degrees')
else
	title('Channel is turned off - No frequency response to plot')
end

% Compare output to input and count errors
binary_err_bits_OFDM = 0;
for i = 1:length(data_in)
	err = abs(data_in(i)-output(i));
	if err > 0
		binary_err_bits_OFDM = binary_err_bits_OFDM +1;
	end
end
BER_OFDM = 100 * binary_err_bits_OFDM/data_length;
disp(strcat('OFDM: BER=', num2str(BER_OFDM,3), ' %'))
disp(strcat('      Number of error bits=', num2str(binary_err_bits_OFDM)))

if (do_QAM == 1)
	disp(strcat('QAM:  BER=', num2str(BER_QAM,3), ' %'))
	disp(strcat('      Number of error bits=', num2str(binary_err_bits_QAM)))
end

% Display text file before and after modulation
if (input_type == 2) & (file_input_type == 2)
	original_text_file = char(data_samples')
	if do_QAM ==1
		edit QAM_text_out.txt
	end
	edit OFDM_text_out.txt
end

% Listen to sounds
if (input_type == 2) & (file_input_type == 3)
	do_again = '1';
	while ( ~(isempty(do_again)) )
		disp(' ')
		disp('Press any key to hear the original sound'), pause
		sound(data_samples,11025)
		disp('Press any key to hear the sound after OFDM transmission'), pause
		sound(output_samples,11025)
		if do_QAM == 1
			disp('Press any key to hear the sound after QAM transmission'), pause
			sound(QAM_output_samples,11025)
		end
		do_again = '';
		do_again = input('Enter "1" to hear the sounds again or press "Return" to end  ', 's');
	end
end