Getting the Fundamental Frequency from Line to Line Voltage without distorting the signal.

Hello Community,

I have obtained the line to line voltage from the inverter test bench and i want to obtain the fundamental waveform. When the try to Y data from the figure, it only extracts one column data. Can Anyone suggest, how can i fix it?

Any suggestion to the code is welcome.

%%Script to convert data from <CSV
data = readtable('data.csv');
t = data{:,1};
x = data{:,2};
y = data{:,3};
z = data{:,4};

hFig1 = figure;
plot(t,x)
title('Line to Line Voltage')
xlabel('t (seconds)')
ylabel('Voltage(V)')
hold on
plot(t,y)
hold on
plot(t,z)
hold off
figure(hFig1);
ax = hFig1.Children;
ln = ax.Children;

xv = ln.XData;
yv = ln.YData;


%Fast Fourier Tranform to get the Fundamental Frequency 
Ts = 1/mean(diff(xv));                              % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency
Fn = Fs/2;                                          % Nyquist Frequency
L = numel(xv);                                      % Signal Length
Y = fft(yv)/L;                                      % Fouriet Transform (Normalised)
Fv = linspace(0, 1, fix(L/2)+1)*Fn;                 % Frequency Vector
Iv = 1:numel(Fv);                                   % Index Vector
figure
plot(Fv, abs(Y(Iv))*2)
grid
xlabel('Frequency')
ylabel('Amplitude')

Answer 1

[Mathieu NOE]

hello

I have not really understood what the purpose of Y data from the figure was, but I guessed that the plan was to get the "analog" voltage waveform (and not the PWM switched one, so I lowpass filtered the data and used some extra function to perform cycle to cycle time measurement (the frequency is basically it's inverse)

does it help ?

all the best

Modified code :

%%Script to convert data from <CSV
data = readtable('data.csv');
t = data{:,1};
x = data{:,2};
y = data{:,3};
z = data{:,4};
hFig1 = figure;
% plot(t,x)
% title('Line to Line Voltage')
% xlabel('t (seconds)')
% ylabel('Voltage(V)')

%%%%%%%%%%%%%%%%
Ts = mean(diff(t));                              % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency

NN = 4;
Wn = 0.005;  % normalized to Nyquist frequency
[B,A] = butter(NN,Wn);

figure(1)
xs = filtfilt(B,A,x);
plot(t,x,t,xs);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Filtered with butterworth LP' ]);
grid on

% zero crossing detection => mesure cycle to cycle time intervals (and derive frequency from time intervals)
threshold = 0;
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(xs,t,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(2)
plot(t,xs,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid on
legend('signal','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');

time_interval = diff(t0_pos);  % cycle to cycle time intervals
freq = [0 1./time_interval];    % frequency from time intervals

figure(3)
plot(t0_pos,freq,'linewidth',2,'markersize',12);grid on
xlabel('Time (s)');
ylabel('Frequency (Hz)');



function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% M Noe
% added positive or negative slope condition

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            slope_sign(ii) = sign(slope);
            t0(ii) = t0(ii) - S(ind(ii)) * slope;
            s0(ii) = level;
        end
    end
end

% extract the positive slope crossing points 
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);

% extract the negative slope crossing points 
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);

end

[ragnor]

@Mathieu NOE Thank you for your time and effort to solve the issue. The idea is not to use a filter but perform a FFT to get the fundamental waveform where i could see the frequency peaks in the spectrum. In my previous code, i could only extract one phase data for the Fourier analysis whereas i needed all the three column data to form a matrix and then perform the FFT.

However, my problem is solved now.

Thanks a lot

[ragnor]

@Mathieu NOE Can you please let me know if i can reconstruct a pure sine wave with the given zero crossing points you calculated in the script above with a given frequency and amplitude. That would be really helpful.

[Mathieu NOE]

here you are , my friend :

%%Script to convert data from <CSV
data = readtable('data.csv');
t = data{:,1};
x = data{:,2};
y = data{:,3};
z = data{:,4};
hFig1 = figure;
% plot(t,x)
% title('Line to Line Voltage')
% xlabel('t (seconds)')
% ylabel('Voltage(V)')

%%%%%%%%%%%%%%%%
Ts = mean(diff(t));                              % Sampling Interval
Fs = 1/Ts;                                          % Sampling Frequency

NN = 4;
Wn = 0.005;  % normalized to Nyquist frequency
[B,A] = butter(NN,Wn);

figure(1)
xs = filtfilt(B,A,x);
plot(t,x,t,xs);legend('Raw','Smoothed');
title(['Data samples at Fs = ' num2str(round(Fs)) ' Hz / Filtered with butterworth LP' ]);
grid on

% zero crossing detection => mesure cycle to cycle time intervals (and derive frequency from time intervals)
threshold = 0;
[t0_pos,s0_pos,t0_neg,s0_neg]= crossing_V7(xs,t,threshold,'linear'); % positive (pos) and negative (neg) slope crossing points 
% t0 => corresponding time (x) values 
% s0 => corresponding function (y) values , obviously they must be equal to "threshold"

figure(2)
plot(t,xs,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid on
legend('signal','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');

time_interval = diff(t0_pos);  % cycle to cycle time intervals
freq = [0 1./time_interval];    % frequency from time intervals

figure(3)
plot(t0_pos,freq,'linewidth',2,'markersize',12);grid on
xlabel('Time (s)');
ylabel('Frequency (Hz)')

%%  get a pseudo sinus waveform

% select crossing points for t > 0.05 s (don't bother with initial
% transient)
t0_pos2 = t0_pos(t0_pos>0.05);
t0_neg2 = t0_neg(t0_neg>0.05);

t0_all = unique(sort([t0_pos2 t0_neg2])); % merge all crossing points , sort them ascending order and keep only unique values
amplitude = max(abs(xs(t>0.3)));  % find raw signal max amplitude (once stabilized)
sign = 1;
xss  = [];
ttt  = [];

for ci = 1:length(t0_all)-1
    t_start = t0_all(ci);
    t_end =   t0_all(ci+1);
    ind = find(t>=t_start & t<t_end);
    tt = t(ind);
    angl = linspace(0,pi,length(ind));
    
    xss_tmp = sign*amplitude*sin(angl);
    xss  = [xss xss_tmp];
    ttt  = [ttt; tt];
    
    sign = -1*sign; % alternate half wave sinus absolute phase
    
end

figure(4)
plot(t,xs,ttt,xss,t0_pos,s0_pos,'+r',t0_neg,s0_neg,'+g','linewidth',2,'markersize',12);grid on
legend('signal','sinus approx','positive slope crossing points','negative slope crossing points');
xlabel('Time (s)');







function [t0_pos,s0_pos,t0_neg,s0_neg] = crossing_V7(S,t,level,imeth)

% [ind,t0,s0,t0close,s0close] = crossing_V6(S,t,level,imeth,slope_sign) % older format

% CROSSING find the crossings of a given level of a signal
%   ind = CROSSING(S) returns an index vector ind, the signal
%   S crosses zero at ind or at between ind and ind+1
%   [ind,t0] = CROSSING(S,t) additionally returns a time
%   vector t0 of the zero crossings of the signal S. The crossing
%   times are linearly interpolated between the given times t
%   [ind,t0] = CROSSING(S,t,level) returns the crossings of the
%   given level instead of the zero crossings
%   ind = CROSSING(S,[],level) as above but without time interpolation
%   [ind,t0] = CROSSING(S,t,level,par) allows additional parameters
%   par = {'none'|'linear'}.
%	With interpolation turned off (par = 'none') this function always
%	returns the value left of the zero (the data point thats nearest
%   to the zero AND smaller than the zero crossing).
%
%	[ind,t0,s0] = ... also returns the data vector corresponding to 
%	the t0 values.
%
%	[ind,t0,s0,t0close,s0close] additionally returns the data points
%	closest to a zero crossing in the arrays t0close and s0close.
%
%	This version has been revised incorporating the good and valuable
%	bugfixes given by users on Matlabcentral. Special thanks to
%	Howard Fishman, Christian Rothleitner, Jonathan Kellogg, and
%	Zach Lewis for their input. 

% Steffen Brueckner, 2002-09-25
% Steffen Brueckner, 2007-08-27		revised version

% Copyright (c) Steffen Brueckner, 2002-2007
% [e-mail]

% M Noe
% added positive or negative slope condition

% check the number of input arguments
error(nargchk(1,4,nargin));

% check the time vector input for consistency
if nargin < 2 | isempty(t)
	% if no time vector is given, use the index vector as time
    t = 1:length(S);
elseif length(t) ~= length(S)
	% if S and t are not of the same length, throw an error
    error('t and S must be of identical length!');    
end

% check the level input
if nargin < 3
	% set standard value 0, if level is not given
    level = 0;
end

% check interpolation method input
if nargin < 4
    imeth = 'linear';
end


% make row vectors
t = t(:)';
S = S(:)';

% always search for zeros. So if we want the crossing of 
% any other threshold value "level", we subtract it from
% the values and search for zeros.
S   = S - level;

% first look for exact zeros
ind0 = find( S == 0 ); 

% then look for zero crossings between data points
S1 = S(1:end-1) .* S(2:end);
ind1 = find( S1 < 0 );

% bring exact zeros and "in-between" zeros together 
ind = sort([ind0 ind1]);

% and pick the associated time values
t0 = t(ind); 
s0 = S(ind);

if strcmp(imeth,'linear')
    % linear interpolation of crossing
    for ii=1:length(t0)
        %if abs(S(ind(ii))) > eps(S(ind(ii)))    % MATLAB V7 et +
        if abs(S(ind(ii))) > eps*abs(S(ind(ii)))    % MATLAB V6 et -    EPS * ABS(X)
            
            % interpolate only when data point is not already zero
            NUM = (t(ind(ii)+1) - t(ind(ii)));
            DEN = (S(ind(ii)+1) - S(ind(ii)));
            slope =  NUM / DEN;
            slope_sign(ii) = sign(slope);
            t0(ii) = t0(ii) - S(ind(ii)) * slope;
            s0(ii) = level;
        end
    end
end

% extract the positive slope crossing points 
ind_pos = find(sign(slope_sign)>0);
t0_pos = t0(ind_pos);
s0_pos = s0(ind_pos);

% extract the negative slope crossing points 
ind_neg = find(sign(slope_sign)<0);
t0_neg = t0(ind_neg);
s0_neg = s0(ind_neg);

end

[ragnor]

@Mathieu NOE Thank you very much. It was indeed very kind of you. I tried to adapt your code and do what i actually wanted but somehow i am not able to do it. My main intention is use a data set unlike the one I gave before without the initial noise. I want to use the following equation to generate the sine wave.

signal = Amplitude*sin(2*pi*F*t);
and t can be calculated by sampling frequency
fs = xxx; % Sampling frequency (samples per second) 
 dt = 1/fs; % seconds per sample 
 StopTime = zzz; % seconds 
 t = (0:dt:StopTime)'; % seconds 

The amplitude and the value of F is known which is (221.1 and 12.56Hz respectively) which i have precisely calculated from Fourier Transformation of the same data set. Now "t" should be such that it crosses over the zero crossings you have generated. The rest of the paramaters are known. This is the script i used for the same data set to generate the frequency and amplitude.

Can you please have a look ? Now i have two scripts which uses the same data set one calculates the amplitude and frequency and second through ur script zero crossing is known. Just need to use both the all together to reconstruct the sine wave.

[Mathieu NOE]

hello again

as you can see from my code the frequency is not 100% stationnary, and if you try to plot a sinewave with only one value for f = 12.56 Hz, there are few chances that your signal will go accross all zero crossing points

FYI, I already tried to plot a sinewave using my own frequency computation, but it was not 100% satisfactory as I had to add another loop to calculate phase vs time so that the sinewave would exactly cross at the desired points

I found this method too complicated so I went to generate more simply the sine wave directly ftom the time values of the zero crossing points

I cannot be more simple than that , and I like the idea to have simple and efficient code !