function [] = equalizerFun(N_in,SNR, eq_type,delay_eq, F_in, imp_resp_nb, show_imp_resp_of_equalizer, longer_sim_4_snr, filename)

% change from filter coefficients to filter order
F = F_in + 1;
% number of BPSK symbols
N = N_in;
% length impulse response
L = 3;
% variance exponential function
V = 1;
% noise variance
n_var = 10.^(-SNR./10);
% equalization selection
%eq_sel = 'ZF-trunc'; % {'ZF-ideal', 'ZF-trunc', 'ZF-trunc-opt', 'MMSE', 'DF'}
switch eq_type
case 0
eq_sel = 'ZF-ideal';
case 1
eq_sel = 'ZF-trunc';
case 2
eq_sel = 'ZF-trunc-opt';
case 3
eq_sel = 'MMSE';
case 4
eq_sel = 'DF';
otherwise
eq_sel = 'ZF-ideal';
endswitch
%% max FIR filter length
%F = 20;
% delay of equalizer
%d_eq = 2;
d_eq = delay_eq;
% IIR equalizer
act_iir_eq = false; % {'true','false'}
% show impulse response of equalizer
if show_imp_resp_of_equalizer==1
	show_imp_resp_eq = true; % {'true','false'}
else
	show_imp_resp_eq = false; % {'true','false'}
end

if longer_sim_4_snr==1
	N_snr = 100000;
else
	N_snr = N+d_eq;
end

%n = 0:N-1+d_eq;
n = 0:N_snr-1+d_eq;

% generate random data
x = randi(2,N_snr,1); % uniformly distributed values 1 and 2
% convert to BPSK
x = 2.*(x-1.5); % uniformly distributed values -1 and 1
x = [x; zeros(d_eq,1)];

figure(1)
stem(n,x)
title('transmit signal (E_s=1)')
% generate impulse response
l = 0:L-1;

switch imp_resp_nb
case 0
	h = [1 -3/4 -5/8];
case 1
	h = [1 1/2 -1/4];
otherwise
	h = [1 -3/4 -5/8];
end
% normalize h
h = h./sqrt(h*h');

figure(2)
stem(l,h)
title('Impulse response')

% compute channel response
r_wo_noise = filter (h, 1,x);
figure(3)
stem(n,r_wo_noise)

% add noise
r = r_wo_noise + normrnd(0,sqrt(n_var),N_snr+d_eq,1);
hold on
stem(n,r,'r');
hold off
legend('received signal','received signal with noise')
title('signal before equalizer');


% compute zeros of transfer function
zeros_tf = roots(h);
if (act_iir_eq)
	stable_roots = zeros_tf(abs(zeros_tf)<1);
	unstable_roots = zeros_tf(abs(zeros_tf)>=1);
else
	stable_roots = [];
	unstable_roots = zeros_tf;
end

stable_equalizer = true;
dirac_impulse = zeros(N_snr,1);
dirac_impulse(1) = 1;

dirac_response = filter(h,1,dirac_impulse);

% for comparison compute always ZF-ideal impulse response of equalizer
ideal_dirac = [zeros(1,d_eq),1]

% compute impulse response (eq) for ZF-ideal
imp_resp_ref = filter(ideal_dirac, h, dirac_impulse);

% do equalization 
switch eq_sel
case 'ZF-ideal' % ideal equalizer (might be unstable)
	
	% perform equalization
	eq_out = filter(ideal_dirac,h,r);
	% compute impulse response (ch+eq)
	imp_resp = filter(ideal_dirac,h,dirac_response);
	% compute impulse response (eq)
	imp_resp_eq = filter(ideal_dirac, h, dirac_impulse);
	
	if ~isempty(unstable_roots)
		stable_equalizer = false;
	end
	
case 'MMSE'
	% IIR part of equalizer (i.e. the stable roots are processed here)
	if ~isempty(stable_roots)
		A = poly(stable_roots);
	else
		A = 1;
	end
	
	% FIR part of equalizer (i.e. the unstable roots are processed here)
	if ~isempty(unstable_roots)
		a = poly(unstable_roots);
		a=a./norm(a);
		% generate linear equation system for determination of equalizer filter coefficients
		M = convmtx(a(end:-1:1),F+length(unstable_roots));
		M = M(:,length(unstable_roots)+1:length(unstable_roots)+F);
		
		% target impulse response
		b = zeros(length(unstable_roots)+F,1);
		b(d_eq+1) = 1;
		
		% get filter coefficients (MMSE solution)
		B = inv(M'*M + n_var.*eye(F))*M'*b
	else
		B = zeros(F,1);
		B(d_eq+1) = 1;
	end	
	
	% perform equalization
	eq_out = filter(B,A,r);
	% compute impulse response (ch+eq)
	imp_resp = filter(B,A,dirac_response);
	% compute impulse response (eq)
	imp_resp_eq = filter(B,A, dirac_impulse);


case 'ZF-trunc-opt' % minimum mean squared error
	% IIR part of equalizer (i.e. the stable roots are processed here)
	if ~isempty(stable_roots)
		A = poly(stable_roots);
	else
		A = 1;
	end
	
	% FIR part of equalizer (i.e. the unstable roots are processed here)
	if ~isempty(unstable_roots)
		a = poly(unstable_roots);
		a=a./norm(a);
		% generate linear equation system for determination of equalizer filter coefficients
		M = convmtx(a(end:-1:1),F+length(unstable_roots));
		M = M(:,length(unstable_roots)+1:length(unstable_roots)+F);
		
		% target impulse response
		b = zeros(length(unstable_roots)+F,1);
		b(d_eq+1) = 1;
	
		% get filter coefficients
		B = M\b % this is the pseudo inverse
	else
		B = zeros(F,1);
		B(d_eq+1) = 1;
	end	
	
	% perform equalization
	eq_out = filter(B,A,r);
	% compute impulse response (ch+eq)
	imp_resp = filter(B,A,dirac_response);
	% compute impulse response (eq)
	imp_resp_eq = filter(B,A, dirac_impulse);
	
case 'ZF-trunc'
	% IIR part of equalizer (i.e. the stable roots are processed here)
	if ~isempty(stable_roots)
		A = poly(stable_roots);
	else
		A = 1;
	end
	
	% FIR part of equalizer (i.e. the unstable roots are processed here)
	if ~isempty(unstable_roots)
		a = poly(unstable_roots);
		a=a./norm(a);
		% generate linear equation system for determination of equalizer filter coefficients
		M = convmtx(a(end:-1:1),F+length(unstable_roots));
		M = M(:,length(unstable_roots)+1:length(unstable_roots)+F);
		
		% target impulse response
		b = zeros(length(unstable_roots)+F,1);
		b(d_eq+1) = 1;
		
		% compute number of samples to force to zeros
		degreeOfFreedom = F-1; % one coefficient is already used for forcing the dirac impulse
		
		numZerosPred = 0;
		numZerosSucc = 0;
		succ = true;
		while (degreeOfFreedom>0)
			if succ % last zero was placed as a successor to the dirac impulse
				if d_eq+1-numZerosPred-1>0
					numZerosPred = numZerosPred + 1;
					degreeOfFreedom = degreeOfFreedom-1;
				end
				succ = false;
			else
				if d_eq+1+numZerosSucc+1<= length(b)-1
					numZerosSucc = numZerosSucc + 1;
					degreeOfFreedom = degreeOfFreedom-1;
				end
				succ = true;
			end
		end
		
		numZerosPred
		numZerosSucc
		d_eq
		size(b)
		
		b = b(d_eq+1-numZerosPred:d_eq+1+numZerosSucc)
		M = M(d_eq+1-numZerosPred:d_eq+1+numZerosSucc , :)
	
		% get filter coefficients
		B = M\b
	else
		B = zeros(F,1);
		B(d_eq+1) = 1;
	end	

	% perform equalization
	eq_out = filter(B,A,r);
	
	% compute impulse response (ch+eq)
	imp_resp = filter(B,A,dirac_response);
	
	% compute impulse response (eq)
	imp_resp_eq = filter(B,A, dirac_impulse);
	
case 'DF'
	delay_line = zeros(F,1);
	eq_coeffs = zeros(F-1,1)
	eq_output = 0;
	eq_out = zeros(length(r),1);
	imp_resp = zeros(length(dirac_response),1)
	
	% load equalizer coefficients
	for d = 1:(min(F,L-d_eq)-1)
		eq_coeffs(d) = h(1+d+d_eq);
	end
	
	% compute output for rx signal
	for c = 1:length(r)
		cur_sample  = r(c)-eq_output;
		eq_out(c) = cur_sample;
		% make decision on current sample
		if cur_sample<0
			dec_cur_sample = -1;
		else
			dec_cur_sample = 1;
		end
		% put to equalizer delay line 
		delay_line(2:end) = delay_line(1:end-1);
		delay_line(1) = dec_cur_sample;
		
		% compute DF output for next samples
		eq_output = 0;
		for d = 1:(min(F,L)-1)
			eq_output = eq_output + delay_line(d)*eq_coeffs(d);
		end
	end
	
otherwise
	error('Equalization method unknown');
endswitch

if ~strcmp(eq_sel,'DF')
	figure(4)
	stem(0:length(imp_resp)-1,imp_resp)
	grid on
	title('Impulse response of channel + equalizer')	
end

figure(5)
stem(n+d_eq,x,'b');
hold on
stem(n,eq_out,'r')
hold off
legend('transmit signal','equalized signal')

% compute signal-to-noise ratio after equalization
delta = eq_out(d_eq+1:end)-x(1:end-d_eq);
	
snr = -10.*log10(mean(delta.^2)); % var(delta)

% compute BER after hard decision

bit_diff = abs(sign(eq_out(d_eq+1:end)) - sign(x(1:end-d_eq)))./2; 
ber = sum(bit_diff)./length(bit_diff);

if stable_equalizer
	title(sprintf('Equalizer is stable, estimated SINR = %g dB',snr))
else
	title(sprintf('Equalizer is NOT stable, estimated SINR = %g dB',snr))
end 


fh = figure(6);

set(0, "defaultaxesfontsize", 12);
set(0, "defaultaxesfontname", "Helvetica");
set(0, "defaulttextfontsize", 12);
set(0, "defaulttextfontname", "Helvetica");
set(0, "defaultaxesfontweight", "bold");
set(0, "defaultaxeslinewidth", 3);
const_lw = 5;
const_ms = 7;


sfh = subplot(2,2,1);
set(sfh,'outerposition',[0 0.5 0.5 0.525]);
stem(n(1:N+d_eq),x(1:N+d_eq),'Color','b','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','b')
title('Transmit Signal (E_s=1)')

xlabel('t / T')
ylabel('Amplitude')
grid on

sfh = subplot(2,2,2);
set(sfh,'outerposition',[0.5 0.5 0.5 0.525])
stem(n(1:N+d_eq),r_wo_noise(1:N+d_eq),'Color',[0 0.65 0],'LineWidth',const_lw,'Marker','d','MarkerSize',const_ms,'markerfacecolor',[0 0.65 0])
hold on
stem(n(1:N+d_eq),r(1:N+d_eq),'Color','b','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','b');
%stem(n,r,'rx');
hold off
legend('w/o noise','w/ noise')
title('Received Signal before Equalizer');
xlabel('t / T')
ylabel('Amplitude')
grid on

sfh = subplot(2,2,3);
set(sfh,'outerposition',[0 0 0.5 0.525])
stem(n(1:N+d_eq),eq_out(1:N+d_eq),'Color','b','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','b')
if (d_eq>0)
	hold on
	stem(n(1:d_eq),eq_out(1:d_eq),'Color','m','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','m')
	hold off
end

hold on
	if ~isempty(find(bit_diff(1:N)==1))
		stem(n(find(bit_diff(1:N)==1))+d_eq,eq_out(find(bit_diff(1:N)==1)+d_eq),'Color','r','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','r')
	end
hold off

title(sprintf('Equalized Received Signal - SINR = %g dB  BER = %g',snr,ber))
xlabel('t / T')
ylabel('Amplitude')
grid on

sfh = subplot(2,2,4);
set(sfh,'outerposition',[0.5 0 0.5 0.525])
if ~strcmp(eq_sel,'DF')
	if ~show_imp_resp_eq
		stem(0:length(imp_resp(1:N+d_eq))-1,imp_resp(1:N+d_eq),'Color','b','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','b');
	else
		if (~strcmp(eq_sel,'ZF-ideal'))
			stem(0:length(imp_resp_ref(1:N+d_eq))-1,imp_resp_ref(1:N+d_eq),'Color',[0.87 0.49 0],'LineWidth',const_lw,'Marker','d','MarkerSize',const_ms,'markerfacecolor',[0.87 0.49 0]);
		end
		hold on
		stem(0:length(imp_resp(1:N+d_eq))-1,imp_resp_eq(1:N+d_eq),'Color','b','LineWidth',const_lw,'Marker','o','MarkerSize',const_ms,'markerfacecolor','b');
		hold off	
	end	
else
	text(0.025, 0.5,'Not available for this equalization method!', 'fontsize', 14, 'color', 'red'); 
end
grid on
xlabel('t / T')
ylabel('Amplitude')
if (show_imp_resp_eq)&&(~strcmp(eq_sel,'ZF-ideal'))
	legend('ZF-ideal',eq_sel)
end

if (show_imp_resp_eq)&&(imp_resp_nb==0)&& (~strcmp(eq_sel,'ZF-ideal'))
	minVal = min(imp_resp_eq(1:N+d_eq));
	maxVal = max(imp_resp_eq(1:N+d_eq));
	minVal = floor(minVal.*2)./2;
	maxVal = ceil(maxVal.*2)./2;
	ylim([minVal maxVal]);
end

if ~show_imp_resp_eq
	title('Impulse Response of Channel + Equalizer')	
else
	title('Impulse Response of Equalizer')
end


%print(filename, '-dpng', '-S930,350', '-tight');
print(filename, '-dpng', '-S1860,700', '-tight');

