%% --- 1.1 Vorgegebene Parameter ---
f_grund = 50;           % Grundfrequenz in Hz
T = 1 / f_grund;        % Periodendauer in Sekunden (20 ms)
U_DC = 100;             % DC-Versorgungsspannung in V
t_schaltpunkte = [0, 1.7e-3, 3.4e-3, 8.3e-3, 10.0e-3, 11.7e-3, 13.4e-3, 18.3e-3, T];
U_out = [0, U_DC/2, U_DC, U_DC/2, 0, -U_DC/2, -U_DC, -U_DC/2];

%% --- 1.2 Periodisches Signal über 40 Perioden erstellen ---
% (Rationaler Ansatz: Vektorisierte Operation)
Anzahl_Perioden = 40;
T_total = Anzahl_Perioden * T;

% Hohe Abtastfrequenz für genaue Auflösung der Schaltflanken
fs = 1e6;   % 1 MHz Abtastfrequenz
ts = 1 / fs;
t = 0:ts:(T_total - ts);

% Vektorisierte Signalerzeugung
t_mod = mod(t, T);
y = interp1(t_schaltpunkte(1:end-1), U_out, t_mod, 'previous', U_out(end));

% --- WICHTIG: DC-Anteil entfernen ---
% Der THD bezieht sich auf die AC-Komponenten.
% Ein DC-Anteil (bei 0 Hz) würde die FFT verfälschen.
y = y - mean(y);

%% --- 2. THD-Berechnung (Standard-MATLAB via FFT) ---
% (Analytischer Ansatz)

fprintf('Berechne THD manuell mittels FFT...\n');

% 2.1 FFT durchführen
L = length(y);          % Länge des Signals (Anzahl Abtastpunkte)
Y = fft(y);             % FFT des Signals berechnen

% 2.2 Einseitiges Amplitudenspektrum berechnen
% Wir brauchen die Amplituden, nicht die komplexen FFT-Werte
P2 = abs(Y/L);          % Zweitseitiges Spektrum (normiert auf Länge)
P1 = P2(1:floor(L/2)+1); % Einseitiges Spektrum
P1(2:end-1) = 2*P1(2:end-1); % Amplituden verdoppeln (außer DC und Nyquist)

% 2.3 Frequenzvektor erstellen
f_vec = fs*(0:floor(L/2))/L;

% 2.4 Frequenzauflösung (Delta-F) der FFT
% Da wir exakt 40 Perioden (0.8s) bei 1MHz Fs analysieren, ist L = 800000
% delta_f = fs / L = 1e6 / 800000 = 1.25 Hz
delta_f = fs / L;

% 2.5 Amplituden der Harmonischen extrahieren
% Wir definieren, wie viele Harmonische wir betrachten (z.B. N=50 -> 2.5 kHz)
N_Harmonische = 50; 
V_amplituden = zeros(1, N_Harmonische);

for k = 1:N_Harmonische
    f_k = k * f_grund;  % Frequenz der k-ten Harmonischen (50, 100, 150...)
    
    % Finde den exakten Index (Bin) im Frequenzvektor
    % +1, da MATLAB 1-indiziert ist und f_vec(1) = 0 Hz (DC)
    idx_k = round(f_k / delta_f) + 1;
    
    % Prüfen, ob der Index im Spektrum liegt (Schutz vor Nyquist-Grenze)
    if idx_k <= length(P1)
        V_amplituden(k) = P1(idx_k);
    else
        % Falls f_k > fs/2, brich ab
        break; 
    end
end

% 2.6 THD berechnen
V_grund = V_amplituden(1); % Amplitude der Grundschwingung (V1)

% Summe der quadrierten Amplituden der Oberschwingungen (V2, V3, ...)
P_harm_sum = sum( V_amplituden(2:end).^2 );

% THD-Formel anwenden
V_harm_rms = sqrt(P_harm_sum);
thd_linear = V_harm_rms / V_grund;

% Umrechnung in Prozent und dB
thd_prozent = thd_linear * 100;
thd_dB = 20 * log10(thd_linear);

fprintf('----------------------------------------------\n');
fprintf('MANUELLE BERECHNUNG (Standard-FFT):\n');
fprintf('Grundschwingung (V1 bei %.1f Hz): %.2f V\n', f_grund, V_grund);
fprintf('RMS der Oberschwingungen (V2-V%d): %.2f V\n', N_Harmonische, V_harm_rms);
fprintf('Total Harmonic Distortion (THD) (in %%):   %.2f %%\n', thd_prozent);
fprintf('Total Harmonic Distortion (THD) (in dB):  %.2f dB\n', thd_dB);
fprintf('----------------------------------------------\n');
THD=thd_prozent;
%% --- 3. Optionale Visualisierung (zur Überprüfung) ---
figure;
stem(f_vec(1:idx_k), P1(1:idx_k)); % Zeige Spektrum bis zur N-ten Harmonischen
title('Einseitiges Amplitudenspektrum (Standard-FFT)');
xlabel('Frequenz (Hz)');
ylabel('Spannung (V)');
xlim([0, f_vec(idx_k) + f_grund]); % Etwas über die letzte Harmonische hinaus zoomen
grid on;