%%%%%%%%%%%%%%%%%%%%%%% Rechteck Fourier %%%%%%%%%%%%%%%%%%%%%%%%
% Thomas Ruschival 2004 ruschi@doomgarage.de
% Rechteck Fourier-Reihe  
% 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all
clear 
Ts=.01;			% SampleZeit
perioden =4;	% anzahl der angezeigten Perioden		
harmonische=10;	% anzahl der harmonischne zur Synthese


%%%%%%%%%%%%%%%%%%%%% 
% Rechtecksignal
%%%%%%%%%%%%%%%%%%%%

%Parameter
T0=1;			% Periode (Zeit) 
Ti=.05;			% Anteil angeschaltet  (in bruchteilen)
PW=Ti*T0;		% Pulsweite 
a0   = 0;            % Gleichanteil
N=T0/Ts;		% Samples / periode 
amplitude=3;	% Signalstaerke

%Simples Signal bauen
%datenvektor, sollte ungerade anzahl haben dass man schön in 2 Teile + 0 teilen kann
if mod(perioden*N,2 )==1
  x=zeros(1,perioden*N)+a0;  % Datenvektor bauen
else
  x=zeros(1,perioden*N+1)+a0;  % Datenvektor bauen  
end 

[tmp, samples]=size(x);        % tmp ist unwichtig, aber die anzahl der samples brauche ich

for per=0:perioden-1
  per_anf = per*N+1	;	 				% erste periode faengt bei 1 an
  idf	  = per_anf+floor(PW*N)+1;	 	% index bis negative flanke kommt
  per_next= (per+1)*N;				 	%  Anfang der naechsten periode
  x(per_anf:idf)= amplitude+a0 ;		 	% von periodenanfang bis zur negativen flanke alles 1
  x((idf+1):per_next)=a0; 				% Von der negativen flanke bis zum ende der Periode alles 0
end;

% x(10:40)=-2;  %<-------- test fuer fehlertoleranz wenn das signal nicht mehr
                % periodisch ist.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% grosse maechtige voellig obskure FFT®
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X_r=fft(x);		% <---- Transformation

%%%%
% Betrag & Phase der harmonischen plotten
%%%%

% Umsortieren
% bie der komplexen fourieranalyse lauft det index von -oo bis +oo
% da matlab aber keinen index -oo hat wird das Transformiert signal bei 0 aufgespalten
% und der negative teil gespiegelt angehaengt. der letzte index im array entspricht
% index -1 in der realität
% die hauptfrequenz ist aber der 1.Index des arrays.

% index basteln der der realität entspricht
idx=[-floor(samples/2):1:floor(samples/2)];

A=X_r(ceil(samples/2)+1:end);
A=[A,X_r(1:ceil(samples/2))];
%A=X_r;

abs_X_r = abs(A)./N;
ang_X_r = angle(A)./pi;


figure;
subplot(2,1,1);
title('Betrag und Phase des Transformierten Vektors');
stem(idx,abs_X_r,'ro');
xlabel('sample');
ylabel('Amplitude der Fourierkoeffizienten');
subplot(2,1,2);
stem(idx,ang_X_r,'kd');
xlabel('sample');
ylabel('Phase [pi]');


%%%%
% raussuchen was wirklich wichtig ist, die kleinen anteile weglassen
%%%%

X_r_temp=X_r; 

for h=1:harmonische
  [amp_X_r(h) ind_X_r(h)]=max(X_r_temp);   % index und wert des groessten heraussuchen
  X_r_temp(ind_X_r(h))=0; 			       % wert loeschen
end

X_r_rueck =zeros(1,samples);        % Vektor bauen der Ruecktransformiert wird
								        
for h=1:harmonische
   X_r_rueck(ind_X_r(h))=amp_X_r(h);     % wichtige werte in transformierte einsetzen
end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% inverse Fouriertransformation
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
xn=abs(ifft(X_r_rueck)); 				% <---- Ruecktransformation

figure
title('Fouriersynthese & Signal');
ylabel('Amplitude');
xlabel('sample')
hold on;
stem([0:1:samples-1],xn,'kd')
plot([0:1:samples-1],x,'r');
hold off;

%%
% Betrag & Phase der zur synthese verwendeten harmonischen plotten
%%
B_r=amp_X_r.*exp(i*2*pi/T0.*ind_X_r);
winkel_r = angle(B_r)./(pi);
betrag_r = abs(B_r)./N;
figure
subplot(2,1,1)
stem([1:1:harmonische],betrag_r,'ro');
xlabel('harmonische')
ylabel('Amplitude der Fourierkoeffizienten')
subplot(2,1,2)
stem([1:harmonische],winkel_r,'kd');
xlabel('harmonische')
ylabel('Phase [pi]')
