function demo_doppler( 
	env_sel,			% environment selection
	channel_sel,			% channel selection	
	plot_time,			% time to plot
	add_fluorescent_effect,		% check to add fluorescent effect to model d and e
	isolate_fluorescent_effect,
	country_sel,			% country selection for PowerLineFrequencys
	alpha,
	filename			% output file for webdemo 
)

% create figure and set style 
fig = figure(1); 
set(fig, "visible", "on"); 
landscape = "-S930,350"; 
portrait  = "-S640,480"; 
output_format = portrait; 
set(0, "defaultlinelinewidth", 2); 
set(0, "defaultaxesfontsize", 8); 
set(0, "defaultaxesfontname", "Helvetica"); 
set(0, "defaulttextfontsize", 8); 
set(0, "defaulttextfontname", "Helvetica"); 
const_lw = 2; 
const_msize = 2; 

[d_bp, sigma_S, tap_delay, tap_power, cl_AoA_rAS_AoD_tAS, cl_power, cl_9, cl_1113, cl_15, cl_1012, cl_13, cl_14, rtp] = gep(env_sel);
[f_min, f_max, fc] = gcp(channel_sel);

% Bell Filter polynomial coefficients taken from Laurent Schuhmacher
%                   Filters designed by Stefano Valle and
%                   Angelo Poloni, STMicroelectronics Srl, Italy
%                   ( Bas Dijkstra and Laurent Schumacher )
NUM_B = [ 2.785150513156437e-004
    -1.289546865642764e-003
    2.616769929393532e-003
    -3.041340177530218e-003
    2.204942394725852e-003
    -9.996063557790929e-004
    2.558709319878001e-004
    -2.518824257145505e-005];
DEN_B = [   1.000000000000000e+000
    -5.945307133332568e+000
    1.481117656568614e+001
    -1.985278212976179e+001
    1.520727030904915e+001
    -6.437156952794267e+000
    1.279595585941577e+000
    -6.279622049460144e-002];

NUM_S= [     7.346653645411014e-002
    -2.692961866120396e-001-3.573876752780644e-002i
    5.135568284321728e-001+1.080031147197969e-001i
    -5.941696436704196e-001-1.808579163219880e-001i
    4.551285303843579e-001+1.740226068499579e-001i
    -2.280789945309740e-001-1.111510216362428e-001i
    7.305274345216785e-002+4.063385976397657e-002i
    -1.281829045132903e-002-9.874529878635624e-003i];

DEN_S=[ 1.000000000000000e+000
    -4.803186161814090e+000-6.564741337993197e-001i
    1.083699602445439e+001+2.655291948209460e+000i
    -1.461669273797422e+001-5.099798708358494e+000i
    1.254720870662511e+001+5.723905805435745e+000i
    -6.793962089799931e+000-3.887940726966912e+000i
    2.107762007599925e+000+1.503178356895452e+000i
    -2.777250532557400e-001-2.392692168502482e-001i];


% Parameters for time-varying impulse response
colormap(jet);
A=9;
NumberOfTaps = length(tap_delay);
v0kmh = 1.2;
v1kmh = 40;
v0ms = v0kmh/3.6;
v1ms = v1kmh/3.6;
lambda = 3*10^8/fc;
fD = v0ms/lambda;
fspike = v1ms/lambda;
f_D_norm     = 1/300;
tap_power = tap_power./sqrt(sum(tap_power));
FadingSamplingFrequency_Hz = fD/f_D_norm;
FadingSamplingTime_s       = 1/FadingSamplingFrequency_Hz;
FadingLength = round((plot_time*10^-3)/FadingSamplingTime_s);
if (FadingLength <= 100) % For faster computing
	f = linspace(f_min, f_max, FadingLength);
else
	f = linspace(f_min, f_max, 100);
endif
tax = linspace(0, plot_time, FadingLength);
Tc= (sqrt(A)/(2*pi*fD))*log(2);
%Tc = 1/fD;

% Parameters for fluorescent lights effect
gm = 0.0203;						% gaussian mean
gsd = 0.0107;						% gaussian standard deviation
AC = [0 -15 -20];					% amplitude modulation coefficients
uipft=linspace(0, plot_time*10^-3, FadingLength);	% uninterpolatedfadingtime-vector
ipft=linspace(0, plot_time*10^-3, 10*FadingLength);	% interpolatedfadingtime-vector
dpos=[2 4 6];						% position of modulated taps in model d
epos=[3 5 7];						% position of modulated taps in model e											
phil = rand(1)*2*pi;					% random phase
X = randn(1)*gsd + gm;					% gaussian distributed random variable for ICrand 
ICrand = X^2;						% Random Interferer to Carrier Ratio
fm = country_sel;					% powerline frequency

% compute modulation function mf
mf_mat = zeros(3,length(ipft));					
for index = 1:3
	mf_mat(index,:) = (10^(AC(index)/20)).*exp(i.*(4.*pi.*(2*(index-1)+1.*fm.*ipft+phil))); % modulation function
endfor
mf = sum(mf_mat); 


% Simulate Doppler Spread by "Filtered Gaussain Noise" Method
% 1) Generate complex gaussian noise (two uncorrelated Gaussian processes with zero mean and common variance 1/2
cgn = sqrt(0.5)*(randn(NumberOfTaps, FadingLength) + i*randn(NumberOfTaps, FadingLength));

cgnt = sqrt(0.5)*(randn(NumberOfTaps, 1000) + i*randn(NumberOfTaps, 1000)); % cgn for transient time of Filter

% 2) Doppler Filtering "shaping the generated complex Gaussian Noise to get the desired Doppler Power Spectrum
% 2.1) Filter Noise with Bell-Shape spectrum
FilterStatesIn = zeros(max(length(NUM_B),length(DEN_B))-1,1);
fmt = zeros(NumberOfTaps, 1000);
for tap_index = 1:NumberOfTaps
	[fmt(tap_index,:), FilterStatesOut] = filter(NUM_B, DEN_B, cgnt(tap_index,:), FilterStatesIn);
endfor

fmt = zeros(NumberOfTaps, FadingLength);
for tap_index = 1:NumberOfTaps
	[fmt(tap_index,:), FilterStatesOut] = filter(NUM_B, DEN_B, cgn(tap_index,:), FilterStatesOut);
endfor

% 2.1) Spike Filter
if (env_sel == 6)
	FilterStatesIn = zeros(max(length(NUM_B),length(DEN_B))-1,1);
	fmts = zeros(1,1000);
	[fmts, FilterStatesOut] = filter(NUM_S, DEN_S, cgnt(3,:), FilterStatesIn);
	fmts = zeros(1,FadingLength);
	[fmts, FilterStatesOut] = filter(NUM_S, DEN_S, cgn(3,:), FilterStatesOut);
	fmt(3,:) = fmts;
endif

% Oversample
if ((env_sel == 3) || (env_sel == 4)) && (add_fluorescent_effect == 1)
	fmtos = zeros(NumberOfTaps, length(ipft));
	for tap_ind = 1:NumberOfTaps
		fmtos(tap_ind,:) = interp1(uipft, fmt(tap_ind,:), ipft, 'linear');
	endfor
% generate a matrix with modulated taps for energy summation (only the modulated taps)
	modulated_fmtos_matrix = zeros(3, length(ipft));
	if(env_sel == 3)
		for dpos_ind = 1:3
			modulated_fmtos_matrix(dpos_ind,:) = fmtos(dpos(dpos_ind),:).*mf;
		endfor
	else
		for epos_ind = 1:3
			modulated_fmtos_matrix(epos_ind,:) = fmtos(epos(epos_ind),:).*mf;
		endfor
	endif	
% Energy summation
	modulated_fmtos_matrix_energy = sum(sum(abs(modulated_fmtos_matrix).^2));
	fmtos_energy = sum(sum(abs(fmtos).^2));
% Calculation of alpha
	ICreal = modulated_fmtos_matrix_energy / fmtos_energy;
	ratio = ICrand/ICreal;
	if (alpha == 0)
		alpha = sqrt(ratio);
	endif
% final fmtos
	fmtos_final = fmtos;
	if (env_sel == 3)
		for index = 1:3
			fmtos_final(dpos(index),:) = fmtos_final(dpos(index),:)+(alpha*modulated_fmtos_matrix(index,:));
		endfor
	else
		for index = 1:3
			fmtos_final(epos(index),:) = fmtos_final(epos(index),:)+(alpha*modulated_fmtos_matrix(index,:));
		endfor
	endif
% h(tau,t) for fmtos_final
	for tap_ind = 1:NumberOfTaps
		hos(tap_ind,:) = fmtos_final(tap_ind,:) .* tap_power(tap_ind);
	endfor
% h(f, t) for fmtos_final
	for fading_unit = 1:FadingLength*10
		H_fos = 0;
		for tap_ind = 1:NumberOfTaps
			H_fos = H_fos + hos(tap_ind, fading_unit)*exp(-j*2*pi*f*tap_delay(tap_ind));
		endfor	
		H_f_dbos(fading_unit,:) = 20*log10(abs(H_fos));
	endfor
% Downsampling of H_f_dbos
	H_f_db_final = zeros(FadingLength,length(f));
	for fad_ind = 1:FadingLength
		H_f_db_final(fad_ind,:) = H_f_dbos(fad_ind*10,:);
	endfor
endif



% h(tau, t)
for tap_ind = 1:NumberOfTaps
	h(tap_ind,:) = fmt(tap_ind,:) .* tap_power(tap_ind);
endfor

% h(f, t)
for fading_unit = 1:FadingLength
	H_f = 0;
	for tap_ind = 1:NumberOfTaps
		H_f = H_f + h(tap_ind, fading_unit)*exp(-j*2*pi*f*tap_delay(tap_ind));
	endfor	
	H_f_db(fading_unit,:) = 20*log10(abs(H_f));
endfor

% Plot
if (add_fluorescent_effect == 0)
	meshz(f, tax, H_f_db);
	zlabel("|H(f,t)|^2 in dB");
	title(['Coherence Time Tc=' num2str(Tc) 's']);

else
	meshz(f,tax,H_f_db_final)
	zlabel("|H(f,t)|^2 in dB");
endif	
view(135,45);
set(gca,'XDir','reverse');
set(gca,'YDir','reverse');
set(gca,'XTick',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]);
set(gca,'XTickLabel',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]/10^6);
set(gca,'YTick',[0 round(max(tax)/5) 2*round(max(tax)/5) 3*round(max(tax)/5) 4*round(max(tax)/5)]);
xlabel(char("                           frequency", "                        in MHz"));
ylabel(char("time", "in ms"));
if (isolate_fluorescent_effect == 1) && (add_fluorescent_effect == 1)
	set(0, "defaultaxesfontsize", 7); 
	dif = (H_f_db_final - H_f_db);
	imagesc(f,tax, dif);
	colormap(jet)
	caxis([-0.15 0.15]);
	h = colorbar('SouthOutside');
	xlabel(h, 'modulation effect in dB')
	title("isolated fluorescent effect");
	set(gca,'XDir','normal');
	set(gca,'YDir','normal');
	set(gca,'XTick',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]);
	set(gca,'XTickLabel',[f_min f_min + (f_max-f_min)/4 f_min + 2*(f_max-f_min)/4 f_min + 3*(f_max-f_min)/4 f_max]/10^6);
	set(gca,'YTick',[0 round(max(tax)/5) 2*round(max(tax)/5) 3*round(max(tax)/5) 4*round(max(tax)/5)]);
	xlabel('frequency in MHz');
	ylabel('time in ms');
endif



% Create an image file 
print(filename, "-dpng",output_format); 
endfunction


