% Skript fourier_reihe1.m, in dem die Fourier-Reihe
% über die DFT für ein Signal mit mehreren 
% Sinus- und Cosinuskomponenten
clear;
% -------- Parameter der Signale
xm = -10;          % Mittelwert
f1 = 100;          % Frequenz der ersten Komponente
ampl1 = 5;         % Amplitude der Komponente
phi1 = pi/3;       % Nullphase
f2 = 200;          % Frequenz der zweiten Komponente
ampl2 = 12;        % Amplitude der Komponente
phi2 = pi/4;       % Nullphase
f3 = 2300;         % Frequenz der dritten Komponente
%f3 = 2300.4;         % Frequenz der dritten Komponente
ampl3 = 6;         % Amplitude der Komponente
phi3 = -pi/5;      % Nullphase
% Diskretisierung
fs = 1000;         % Abtastfrequenz
Ts = 1/fs;         % Abtastperiode
N = 1000;          % Anzahl der Abtastwerte im Untersuchungsintervall
n = 0:N-1;
k = n;
% ------- Zeitdiskretes Signal
xn = xm + ampl1*cos(2*pi*f1*n/fs+phi1)+ ...
          ampl2*cos(2*pi*f2*n/fs+phi2)+ ...
          ampl3*cos(2*pi*f3*n/fs+phi3);
fm = max([f1, f2, f3]);     Tmin = 1/fm;
dt = Tmin/10;
t = 0:dt:(N-1)*Ts;
xt = xm + ampl1*cos(2*pi*f1*t+phi1)+ ...
          ampl2*cos(2*pi*f2*t+phi2)+ ...
          ampl3*cos(2*pi*f3*t+phi3); % zeitkontinuierliches Signal
figure(1);     clf;
text = ['(f1 = ',num2str(f1),' Hz; f2 = ',num2str(f2),...
    ' Hz; f3 = ',num2str(f3),' Hz )'];
subplot(311), plot(t, xt);
hold on;      stem(n*Ts,xn, 'LineWidth',1.5);
axis tight;
La = axis;    axis([0.1, 0.12, La(3:4)]);
hold off;
title('Kontinuierliches Signal und Abtastwerte (Ausschnitt)');
xlabel('Zeit in s');   grid on;   
% ------- Amplitudenspektrum über die DFT
Xk = fft(xn);
betrag_Xk = abs(Xk);    phase_Xk = angle(Xk);
p = find(abs(real(Xk))<1e-8 & abs(imag(Xk))<1e-8);
phase_Xk(p) = 0;  % Entfernung der Fehler in der Phasenrechnung
subplot(312), stem(k, betrag_Xk/N,  'LineWidth',1.2);
title('Betrag der DFT/N');    
xlabel(['Frequenz in Hz (fs = ',num2str(fs),' Hz)']);
ylabel('|DFT|/N');     grid on;
subplot(313), stem(k, phase_Xk,  'LineWidth',1.2);
title('Phase der DFT');    
xlabel(['Frequenz in Hz (fs = ',num2str(fs),' Hz)']);
ylabel('Rad');     grid on;


