#not a function file:
1;

global f_sig = 1000;	#signal frequency
global f_sample = 2200;	#sample frequency

global T_sig = 1/f_sig;
global T_sample = 1/f_sample;

#define signal
global sig_form = "rectangle";	#possible values: sine, rectangle

function retval = signal (t)
global sig_form
if (strcmp(sig_form, "sine"))
global f_sig
retval = sin(2*pi()*f_sig*t);
elseif (strcmp(sig_form, "rectangle"))
global T_sig
retval = floor ((t + T_sig/2)/T_sig) - floor (t/T_sig);
else
retval = t;
endif
endfunction

#calculate needed sample_time to get a 'sample period'
function retval = sample_time ()
global T_sig
global T_sample
n = 1;
m = 1;
while (n*T_sample != m*T_sig)
	if (n*T_sample < m*T_sig)
		n++;
	else
		m++;
	endif
endwhile
retval = n*T_sample;
endfunction

#caluclate a_n fourier-coefficients
function retval = a (n, w, t, signal_sampled)
retval = 0;
for i = 1:n
coeffs(i) = 2/sample_time() * trapz(t, signal_sampled .* cos (i * w * t));
endfor
for i = 1:n
retval += coeffs(i) * cos (i * w * t);
endfor
endfunction

#calculate b_n fourier-coefficients
function retval = b (n, w, t, signal_sampled)
retval = 0;
for i = 1:n
coeffs(i) = 2/sample_time() * trapz(t, signal_sampled .* sin (i * w * t));
endfor
for i = 1:n
retval += coeffs(i) * sin (i * w * t);
endfor
endfunction

#calcualte sample-function with ideal samples at the beginning of the sample period
function retval = ideal_sample_beg (t)
global T_sample
sample_pulses = diff(floor(t/T_sample))./diff(t);
retval = (signal(t))(2:end).*sample_pulses;
endfunction

#calcualte sample-function with ideal samples at the middle of the sample period
function retval = ideal_sample_mid (t)
global T_sample
sample_pulses = diff(floor((t+T_sample/2)/T_sample))./diff(t);
retval = (signal(t))(2:end).*sample_pulses;
endfunction

#calculate sample-function with real samples at the beginning of the sample period
function retval = real_sample_beg (t)
global T_sample
retval = signal (floor (t(2:end)/T_sample)*T_sample);
endfunction

#calculate sample-function with real samples at the middle of the sample period
function retval = real_sample_mid (t)
global T_sample
retval = signal (floor ((t(2:end)+T_sample/2)/T_sample)*T_sample);
endfunction

#calculate filtered function
function retval = calc_fourier (t_sample, sampled, t)
global f_sample
w_sample=(2*pi())/t_sample;
n_fourier = floor (t_sample*f_sample)/2;
a_0=2/t_sample * trapz(t(2:end), sampled);
a_n = a (n_fourier, w_sample, t(2:end), sampled);
b_n = b (n_fourier, w_sample, t(2:end), sampled);
retval = a_0/2 + a_n + b_n;
endfunction

#define time
increment = 0.000001;
t_sample = sample_time() - increment;
t=0:increment:t_sample;

#sampled = real_sample_beg (t);
#sampled = real_sample_mid (t);
#sampled = ideal_sample_beg (t);
sampled = ideal_sample_mid (t);

plot (t, signal(t), ";Eingangsignal;", t(2:end), sampled, ";abgetastetes Signal;", t, 0.1*floor (t/T_sample), ";Abtastperioden;", t(2:end), 1/2000*calc_fourier(t_sample, sampled, t), "linewidth", 5, ";Ausgangssignal;"); 
axis ([-0.001 (t_sample + 0.001) -0.5 2])
clear all
