clear; close all

%% Generate signal

nharodd = 100;
Fs = 1E3;
t = 0:1/Fs:10-1/Fs;
X = nan(length(t), nharodd);
for i = 1:nharodd
    X(:, i) = 1/(2*i-1)*sin(2*pi*(2*i-1)*5*t);
end
x = sum(X, 2);
figure(1); subplot(3, 1, 1); plot(t, x)
xlabel('Time (s)'); ylabel('B (T)')

Y = fftshift(fft(x));
L = size(Y, 1);
P1 = 2*abs(Y/L);
P1((length(P1))/2+1) = P1((length(P1))/2+1)/2;
phs = angle(Y); phs(P1 < 1E-4) = 0;
f = Fs*(-L/2:(L/2)-1)/L;

subplot(3, 1, 2); plot(f, P1, 'r')
xlabel('Frequency (Hz)'); ylabel('B (T)');
subplot(3, 1, 3); plot(f, phs, 'r')
xlabel('Frequency (Hz)'); ylabel('Phase (radians)')

%% Modify signal by adding some even harmonic components

nhareven = 100;
Z = nan(length(t), nhareven);
for i = 0:100
    for j = 1:nhareven
        Z(:, j) = 1/(2*j)*sin(2*pi*(2*j)*5*t);
    end
    z = x + i*0.01*sum(Z, 2);
    figure(2); subplot(3, 1, 1); plot(t, z)
    xlabel('Time (s)'); ylabel('B (T)')
    title(sprintf('Relative strength of even harmonics: %2.2f', 0.01*i))

    Y = fftshift(fft(z));
    L = size(Y, 1);
    P1 = 2*abs(Y/L);
    P1((length(P1))/2+1) = P1((length(P1))/2+1)/2;
    phs = angle(Y); phs(P1 < 1E-4) = 0;
    f = Fs*(-L/2:(L/2)-1)/L;

    subplot(3, 1, 2); plot(f, P1, 'r')
    xlabel('Frequency (Hz)'); ylabel('B (T)');
    subplot(3, 1, 3); plot(f, phs, 'r')
    xlabel('Frequency (Hz)'); ylabel('Phase (radians)')

    pause(0.01)
end

figure; yyaxis left; plot(t, [0; diff(z)*Fs]);
xlabel('Time (s)'); ylabel('H (A/m)')
yyaxis right; plot(t, cumtrapz([0; diff(z)*Fs]))
ylabel('B-hat (T)')

%% Add DC bias to B

q = x + 0.5;

figure; subplot(3, 1, 1); plot(t, q);
xlabel('Time (s)'); ylabel('B (T)')

Y = fftshift(fft(q));
L = size(Y, 1);
P1 = 2*abs(Y/L);
P1((length(P1))/2+1) = P1((length(P1))/2+1)/2;
phs = angle(Y); phs(P1 < 1E-4) = 0;
f = Fs*(-L/2:(L/2)-1)/L;

subplot(3, 1, 2); plot(f, P1, 'r')
xlabel('Frequency (Hz)'); ylabel('B (T)');
subplot(3, 1, 3); plot(f, phs, 'r')
xlabel('Frequency (Hz)'); ylabel('Phase (radians)')

%% Visualize hysteresis loops for a sinusoidal input

% H leads B
h = sin(2*pi*5*t - pi/10);

figure; plot(h, x)
hold on; plot(h, z)
plot(h, q)
xlabel('H (A/m)'); ylabel('B(T)')