clear all; %resetet Workspace, alle Variablen werden auf null gestetzt
clc %löscht alte Fehlermeldungen im commando window

%--- Variablen und Konstanten
f2=120; % Oberschwingungsfrequenz der Testfunktion
f1=50; % Grundschwingungsfrequenz der Testfunktion.

T1=3/f1         % Anzahl der Periodendauer, des Testsignals
fa=2*f1*100;    % freie Wahl einer Abtastfrequenz unter berücksichtigung des Abtasttherorems;
ta=1/fa;        % Abtastdauer
N=T1/ta;        % Anzahl der Stützstellen pro periode, je größer N, desto besser die Auflösung (weniger leckeffekt)
               

%--- Messreihe der Testfunktion
t=0:ta:T1; %Zeitvektor
y=40.*sin(2*pi*f1.*t)+10.*cos(2*pi*f2.*t);


%--- Plotten des Testsignals und deren Abtastung
figure(1)
subplot(3,1,1)
plot(t,y)
hold on
stem(t(1:5:length(t)+1),y(1:5:length(y)+1))
grid on;
hold off


%--- Aufstellen der FFT
X_n=fft(y,N);
% Entspricht: sum(x(k)*exp(-j*2*pi*n*k/N))
% über die Indizes k=0..N-1, für die Ordnung (n) -inf<n<inf; 1/N=ta/T_m

c_n=abs(X_n)./N;
% Matlab löst die FFT nummerisch. Unter der Voraussetzung, das Mess- und
% Grundschwingungsperiodendauer übereinstimmen (T_m=T_1, sog.
% synchronisierte Messdauer oder kohärente Abtastung), lässt sich die
% Koeffizientenberechnung annähern zu c_n=1/N*sum(x(k)*exp(-j*2*pi*n*k/N))
% über die Indizes k=0..N-1, für die Ordnung (n) -inf<n<inf; 1/N=ta/T_m

Y=c_n.*2; 
% Die Amplitude(Y) entspricht für Ordnungen(n)=!0 zwei mal den Koeffizienten 2*c_n; 
% (für A_n=0=c_n)          
 


%--- Plotten der FFT, des Testsignals
%Frequenzvektor, es reicht nur die Hälfte des Ergebnisvektors der FFT darzustellen, 
% aufgrund von Symetrie
f=(1/T1).*(0:(N/2)); 
subplot(3,1,2)
stem(f,Y(1:(N/2+1)));
hold on
plot(f,Y(1:(N/2+1)));
axis([0 200 0 max(Y)+1]);
grid on;
hold off


%--- Rücktrafo
yr=ifft(X_n,N); %Rückwärtstrafo
subplot(3,1,3)
plot((1:N).*ta,yr);
%axis([0 200 -10 10])
grid on;



%--- Spectralanalyse
figure(2)
specgram(y,[],240);
colorbar
colormap;