FFT, THD and Fourier equations

Hi guys! I have the following code, can I get the Fourier equations like in the "screen2" for maybe 50 harmonics with correct coefficients and get correct fitting in plots? In this code, I get Fourier equations for only one harmonic component and only sin without cos part. I need the form of equations like this: (a(1)*sin(2*pi*f1*t)+a(2)*cos(2*pi*f1*t))+mean(y2); Thanks in advance!

data = readmatrix('tek0005.csv', 'Delimiter', ';', 'HeaderLines', 21);
time = data(:, 1);
voltage = data(:, 2);
current = data(:, 3);
fs = 1 / (time(2) - time(1));
N = length(time);
frequencies = (0:N-1) * fs / N;
voltage_fft = fft(voltage);
current_fft = fft(current);
a1_voltage = abs(voltage_fft(2)) / N;
a1_current = abs(current_fft(2)) / N;
total_harmonic_distortion_voltage = sqrt(sum(abs(voltage_fft(2:end)).^2)) / abs(voltage_fft(2));
total_harmonic_distortion_current = sqrt(sum(abs(current_fft(2:end)).^2)) / abs(current_fft(2));
f1 = frequencies(2);
equation_voltage = @(t) a1_voltage * sin(2 * pi * f1 * t) + mean(voltage);
equation_current = @(t) a1_current * sin(2 * pi * f1 * t) + mean(current);
fprintf('Fourier equation for voltage: %.4f * sin(2 * pi * %.4f * t) + %.4f\n', a1_voltage, f1, mean(voltage));
fprintf('Fourier equation for current: %.4f * sin(2 * pi * %.4f * t) + %.4f\n', a1_current, f1, mean(current));
fprintf('For voltage: a1 = %.4f, THD = %.4f\n', a1_voltage, total_harmonic_distortion_voltage);
fprintf('For current: a1 = %.4f, THD = %.4f\n', a1_current, total_harmonic_distortion_current);
t = linspace(min(time), max(time), 1000); % Create 1000 time points
voltage_fit = equation_voltage(t);
current_fit = equation_current(t);
figure;
subplot(2, 1, 1);
plot(time, voltage, 'b', t, voltage_fit, 'r');
title('Voltage over time');
xlabel('Time [s]');
ylabel('Voltage [V]');
legend('Data', 'Fitting');
subplot(2, 1, 2);
plot(time, current, 'b', t, current_fit, 'r');
title('Current over time');
xlabel('Time [s]');
ylabel('Current [A]');
legend('Data', 'Fitting');
phase_voltage = angle(voltage_fft);
phase_current = angle(current_fft);
phase1_voltage = phase_voltage(2);
phase1_current = phase_current(2);
fprintf('Phase angle for voltage: %.4f radians\n', phase1_voltage);
fprintf('Phase angle for current: %.4f radians\n', phase1_current);

Answer 1

[Sulaymon Eshkabilov]

The solution to your exercise comes from the nonlinear fit model and fft()'s fundamental resonant frequency value, which is 2Hz.

data = readmatrix('tek0005.csv', 'Delimiter', ';', 'HeaderLines', 21);
time = data(:, 1);
voltage = data(:, 2);
current = data(:, 3);
fs = 1 / (time(2) - time(1));
N = length(time);
frequencies = (0:(N/2)) * fs / N;   % It is better to take one half. SInce the other half is just the reflection
voltage_fft = fft(voltage);
current_fft = fft(current);

V = abs(voltage_fft/N);
C = abs(current_fft/N);
% It is better to take one half. Since the other half is just the reflection
V1 = V(1:N/2+1);
V1(2:end-1) = 2*V1(2:end-1);
C1 = C(1:N/2+1);
C1(2:end-1) = 2*C1(2:end-1);
f1  = 2;   % found from fft
Model = @(a, t)(a*sin(2*pi*f1*t));   % two terms can be also considered: (a(1)*sin(2*pi*f1*t)+a(2)*sin(2*pi*f2*t))
t  = time;
y1 = voltage; 
y2 = current; 
a0 = 1;       % Initially estimated guess value for coeff a
Coeff_V = nlinfit(t, y1, Model, a0);
Coeff_C = nlinfit(t, y2, Model, a0);
a1_voltage = Coeff_V; % Found fit model coeff
a1_current = Coeff_C; % Found fit model coeffs
equation_voltage = @(t) a1_voltage * sin(2 * pi * f1 * t) + mean(voltage);
equation_current = @(t) a1_current * sin(2 * pi * f1 * t) + mean(current);
fprintf('Fourier equation for voltage: %.4f * sin(2 * pi * %.4f * t) + (%.4f)\n', a1_voltage, f1, mean(voltage));

Fourier equation for voltage: -2.7014 * sin(2 * pi * 2.0000 * t) + (-0.1160)

fprintf('Fourier equation for current: %.4f * sin(2 * pi * %.4f * t) + (%.4f)\n', a1_current, f1, mean(current));

Fourier equation for current: -0.0335 * sin(2 * pi * 2.0000 * t) + (-0.1750)

y_fit_V = (a1_voltage*sin(2*pi*f1*t))+ mean(voltage);
y_fit_C = (a1_current*sin(2*pi*f1*t))+ mean(current);
figure(1)
plot(t, voltage, 'r')
hold on
plot(t, y_fit_V, 'k--', 'LineWidth',2)
legend('Data: Voltage', 'Fit Model')
ylabel('Voltage')
xlabel('Time')
hold off

...
% Similarly, you can do for the current data
figure(1)
plot(t, current, 'b')
hold on
plot(t, y_fit_C, 'k--', 'LineWidth',2)
legend('Data: Current', 'Fit Model')
ylabel('Current')
xlabel('Time')
hold off

... % The rest is your code for phase and other calcs and displays
    

[Seweryn]

Hi @Sulaymon Eshkabilov, Can you check if the results and plots for these files are correct? Strange results come out I think. The second question is, I wonder if we can get the Fourier equations with more than one coefficient like @(t) a1_voltage * sin(2 * pi * f1 * t)...a10_voltage * sin(2 * pi * f1 * t)?

[Sulaymon Eshkabilov]

The fit model depends on the fundamental and other resonant frequency values. I will check your posted data and get back with my answer.

[Seweryn]

@Sulaymon Eshkabilov You can see the frequencies in the screenshot. I tried to add more files but the daily limit is 10, I deleted a few and still the same problem.

I

Answer 2

[Sulaymon Eshkabilov]

Note that your different data sets has different fundamental frequencies and other resonant freq values. Besides, some of your data has some offsets (e.g. tek0006.csv). Therefore, for each individual set of data the fit model freq values and number of coefficients need to be selected based on fft analysis. Besides, the resonant freq values for voltage and current differs. That might be result of inconsistencies in data collection or something else (a bit difficult to guess without seeing the real procedure of data collection). See the simulation code how to compute the fit models for voltage and current in separate due to the differences of v(t) and i(t) - in in the example of tek0006.csv.

data = readmatrix('tek0006.csv', 'Delimiter', ';', 'HeaderLines', 21);
time = data(:, 1);
voltage = data(:, 2);
current = data(:, 3)-mean(data(:, 3));  % Current data hasa some offset
fs = 1 / (time(2) - time(1));
N = length(time);
frequencies = (0:(N/2)) * fs / N;   % It is better to take one half. SInce the other half is just the reflection
voltage_fft = fft(voltage);
current_fft = fft(current);

V = abs(voltage_fft/N);
C = abs(current_fft/N);
% It is better to take one half. Since the other half is just the reflection
V1 = V(1:N/2+1);
V1(2:end-1) = 2*V1(2:end-1);
C1 = C(1:N/2+1);
C1(2:end-1) = 2*C1(2:end-1);
figure
plot(frequencies, V1)  % See this FFT plot which has five significant resonant freq values for voltage data
xlim([0, 500])
ylabel('|fft(V)|')
xlabel('frequency, [Hz]')
grid on

%
figure
plot(frequencies, C1) % % See this FFT plot which has five significant resonant freq values for voltage data
xlim([0, 2500])
ylabel('|fft(I)|')
xlabel('frequency, [Hz]')
grid on

%% Fit models
f1 = 50;   % found from fft
f2 = 98;   % found from fft
f3 = 102;  % found from fft
f4 = 198;  % found from fft
f5 = 250;  % found from fft
Model = @(a, t)(a(1)*sin(2*pi*f1*t)+a(2)*sin(2*pi*f2*t)+a(3)*sin(2*pi*f3*t)+a(4)*sin(2*pi*f4*t)+...
               +a(5)*sin(2*pi*f5*t));   
t  = time;
y1 = voltage; 
a0 = [1, 1, 1, 1, 1];       % Initially estimated guess value for coeff a
Coeff_V = nlinfit(t, y1, Model, a0);
a1_voltage = Coeff_V; % Found fit model coeff
equation_voltage = @(t) a1_voltage(1) * sin(2 * pi * f1 * t) + a1_voltage(2) * sin(2 * pi * f2 * t) +...
                        a1_voltage(3) * sin(2 * pi * f3 * t) + a1_voltage(4) * sin(2 * pi * f4 * t) + ...
                        a1_voltage(5) *sin(2*pi*f5*t)+ mean(voltage);
y_fit_V = (a1_voltage(1)*sin(2*pi*f1*t))+(a1_voltage(2)*sin(2*pi*f2*t))+...
           a1_voltage(3) * sin(2 * pi * f3 * t) + a1_voltage(4) * sin(2 * pi * f4 * t)+...
           a1_voltage(5) *sin(2*pi*f5*t)+mean(voltage);
figure()
plot(t, voltage, 'r')
hold on
plot(t, y_fit_V, 'k--', 'LineWidth',2)
legend('Data: Voltage', 'Fit Model')
ylabel('Voltage')
xlabel('Time')
hold off
xlim([0, max(time)])

...
% Similarly, you can do for the current data
f1 = 569;   % found from fft
f2 = 1137;  % found from fft
f3 = 1706;  % found from fft
f4 = 2275;  % found from fft
t  = time;
y2 = current; 
a0 = [1, 1, 1, 1];       % Initially estimated guess value for coeff a
Model = @(a, t)(a(1)*sin(2*pi*f1*t)+a(2)*sin(2*pi*f2*t)+a(3)*sin(2*pi*f3*t)+a(4)*sin(2*pi*f4*t));   
Coeff_C = nlinfit(t, y2, Model, a0);
a1_current = Coeff_C; % Found fit model coeffs

equation_current = @(t) a1_current(1) * sin(2 * pi * f1 * t) + a1_current(2) * sin(2 * pi * f2 * t)+...
                        a1_current(3) * sin(2 * pi * f2 * t) + a1_current(4) * sin(2 * pi * f2 * t);

y_fit_C = (a1_current(1)*sin(2*pi*f1*t))+(a1_current(2)*sin(2*pi*f2*t))+...
           a1_current(3) * sin(2 * pi * f3 * t) + a1_current(4) * sin(2 * pi * f4 * t) ;

figure()
plot(t, current, 'b')
hold on
plot(t, y_fit_C, 'k--', 'LineWidth',2)
legend('Data: Current', 'Fit Model')
ylabel('Current')
xlabel('Time')
hold off
xlim([0, max(time)/10])

%... % Display of of results and phase calcs you can figure out

Answer 3

[Sulaymon Eshkabilov]

You need to have a look at tek0025.csv a bit more careful. Both voltage and current seem to have the same fundamental frequency of ~5Hz. However, consiistenciy of the meaured data doesn't appear to be sufficiently strong.

data = readmatrix('tek0025.csv', 'Delimiter', ';', 'HeaderLines', 21);
time = data(:, 1);
voltage = data(:, 2);
current = data(:, 3);
fs = 1 / (time(2) - time(1));
N = length(time);
frequencies = (0:(N/2)) * fs / N;   % It is better to take one half. SInce the other half is just the reflection
voltage_fft = fft(voltage);
current_fft = fft(current);

V = abs(voltage_fft/N);
C = abs(current_fft/N);
% It is better to take one half. Since the other half is just the reflection
V1 = V(1:N/2+1);
V1(2:end-1) = 2*V1(2:end-1);
C1 = C(1:N/2+1);
C1(2:end-1) = 2*C1(2:end-1);
figure
plot(frequencies, V1)  % See this FFT plot which has five significant resonant freq values for voltage data
xlim([0, 500])
ylabel('|fft(V)|')
xlabel('frequency, [Hz]')
grid on

%
figure
plot(frequencies, C1) % % See this FFT plot which has five significant resonant freq values for voltage data
xlim([0, 500])
ylabel('|fft(I)|')
xlabel('frequency, [Hz]')
grid on

f1  = 5;   % found from fft
Model = @(a, t)(a*sin(2*pi*f1*t));   
t  = time;
y1 = voltage; 
y2 = current; 
a0 = 1;       % Initially estimated guess value for coeff a
Coeff_V = nlinfit(t, y1, Model, a0);
Coeff_C = nlinfit(t, y2, Model, a0);
a1_voltage = Coeff_V; % Found fit model coeff
a1_current = Coeff_C; % Found fit model coeff
equation_voltage = @(t) a1_voltage * sin(2 * pi * f1 * t) ;
equation_current = @(t) a1_current * sin(2 * pi * f1 * t) ;
fprintf('Fourier equation for voltage: %.4f * sin(2 * pi * %.4f * t)+ (%.4f) \n', a1_voltage, f1, mean(voltage));

Fourier equation for voltage: -0.6070 * sin(2 * pi * 5.0000 * t)+ (-0.9302) 

fprintf('Fourier equation for current: %.4f * sin(2 * pi * %.4f * t)+ (%.4f) \n', a1_current, f1, mean(current));

Fourier equation for current: -0.0075 * sin(2 * pi * 5.0000 * t)+ (-0.1821) 

y_fit_V = (a1_voltage*sin(2*pi*f1*t))+mean(voltage);
y_fit_C = (a1_current*sin(2*pi*f1*t))+mean(current);
figure()
plot(t, voltage, 'r')
hold on
plot(t, y_fit_V, 'k--', 'LineWidth',2)
legend('Data: Voltage', 'Fit Model')
ylabel('Voltage')
xlabel('Time')
hold off

...
% Similarly, you can do for the current data
figure()
plot(t, current, 'b')
hold on
plot(t, y_fit_C, 'k--', 'LineWidth',2)
legend('Data: Current', 'Fit Model')
ylabel('Current')
xlabel('Time')
hold off

%... % The rest is your code for phase and other calcs and displays
    

[Seweryn]

@Sulaymon Eshkabilov Thank you very much, I really appreciate your work. I will try to analyse the rest of the files. Best regards!

[Sulaymon Eshkabilov]

Most welcome. It is my pleasure. All the best with your studies!