how to generate a Sine wave with changable frequency in mfile?

Dear all,

I would like to generate something similler to a PLL sine signal locked on a measured sine signal. for that i have detected the Zeros of my measured, lets say, voltage signal and calculated the freuqncy from counting the number of points between two zero crossing. (the frequency of the measured signal is not constant and changing with time). by using the following equation that was achieved:

f = 1./(2*diff(Zeros));% where Zeros are my detected zerocorssing times in s.

so by knowing the following:

1- I know the frequency of my measured signal (f vector length is not equal to the length of measured signal)

2- I know my sample rate of my measured signal and my time.

3- sin(2*pi*f*t) will generate my locked sine signal to the measured sine wave signal (Voltage)

my question: why am i not getting a correct sine wave signal with the same frequency of the measured signal?

note: Zeros is a vector containing the interpolated time values when my sine wave signal corss the zero value. (so i only have frequency values at the diff(Zeros) values and not at the same sample rate of the measured signal).

Can anyone help me to condition my calculated frequency signal to have the same lenght as the measured signal (and the correct time steps) in order to make the sin(2*pi*f*t) works correctly?

Best Regards

[dpb]

Can't really tell w/o some actual code and a (small) dataset that illustrates the problem, but --

a) if it's not constant frequency your description seems to only have an average over the whole period so it wouldn't ever match any one portion, and

b) zeros is a builtin (and heavily used) Matlab function; use something else for your data array of crossing points to avoid high confusion/unintended behavior.

[Aubai]

thx for trying to help my code looks like the following:

if true
  % code
function [f,f_time] = Frequency_Calculation_Function(Singal,time)
% This function will locate the zero crossing in your input signal and
% calculate the frequency of that input signal.
% Author: Aubai Al Khatib datum: 19.09.2013 
Zeros = Zeros_finding(Singal,time);
f = 1./(2*diff(Zeros));
f_time = Zeros(1:length(f));
end

function [xc] = Zeros_finding(y,x) n = length(y); ind = 1:(n-1);

% list of intervals with a zero crossing
k = find(((y(ind)<=0) & (y(ind+1)>0)) | ...
         ((y(ind)>=0) & (y(ind+1)<0)));

% list of zero crossings
xc = [];

% intervals where y is zero at both ends of the
% interval are exactly zero are indeterminate.
% the solution may be anywhere on that interval.
% I'll choose to return both endpoints of the
% interval.
L = (y(k)==0) & (y(k+1)==0);
if any(L) 
  xc = x(k(L));
  k(L)=[];
end

% interpolate to find x. I've already removed
% the constant intervals at zero
if ~isempty(k)
  s = (y(k+1)-y(k))./(x(k+1)-x(k));
  xc = [xc,x(k) - y(k)./s];
end

% patch for last element exactly zero
if y(end)==0
 xc(end+1) = x(end);
end
if xc(1) == 0
xc = xc(2:end);
end
end

end

this is the easy version of the code where i show how i calcualted the frequency.

but in the following code i show the complete calculation of the sync. sine wave genertated from that calculated frequency:

if true
  % code
function [f,f_time,f_interp,time_Zeros,Wt,Sine,Cose,Signal_Zeros,Signal_Limits,PLL_time] = PLL_final_cal(Signal,time,tol,option)
% This function will locate the zero crossing in your input signal and
% calculate the frequency of that input signal.
% This function will also generate your input signal 'Wt' function and
% regenerate an pu sin & cos signal just like an PLL_Block in simulink.
% *************************************************************************
% the needed inputs:
% 1- Signal: is your input y-axi data (in case of sine wave those are the
% Amplitude)
% 2- time: is your input x-axi data which is the time space of the input
% signal
% 3- tol: is the tolerance in which zeros points will be selected. and it
% should be given in Hz.
% for example if tol = 1 [Hz] is means any Zeros(i+1) point which distance
% from the point before Zeros(i) is in this tolerance band of
% [1/(fn-tol),1/fn,1/(fn+tol)] in S will be accepted as the next frequency
% of f(i+1).
% 4- option: this variable will allow you to choose between doing the PLL
% evaluation according to the batlab simulink block (option = 1) or to
% generate the Wt signal without intgration (option ~= 1)
% both methods should give the same results yet the second one (option ~=1)
% is mor accurate one as Wt is directly built from the linearly
% interpolated zerocorssing function before
% *************************************************************************
% the available outputs:
% 1- f: is the frequency of the input Signal calculated at every Zeros(i)
% points which means f length is the same as the number of Zeros found in
% the input Signal divided by 2 (in the Sine wave the half wave part is
% neglected). mins 1 (due to the "diff" function effect).
% 2- f_time: is the time signals which can be used to "plot" this f signal
% and returns the times where the Zeros has occured and the frequency were
% calculated.
% f_interp & time_Zeros: are the interpulation of the previous signals in
% order to get the output frequency with the same size of the input Signal.
% 3- Wt, Sine & Cose: are the same output of the PLL function block in
% simulink. which means the are the pu represntation of the input Signal
% where there is no harmonics. yet have the same frequency of the input
% Signal.
% 4- Signal_Zeros: is the cut version of the input Signal (to synchronize
% on the first or second Zero crossing)
% 5- Signal_Limits: is the indx representation of the input signals with
% only two values [Start_Value End_Value], in which the complete output was
% created.
% this function file will need the following function files in order to
% work:
% 1- Fast_Harmonics.m: is one is used to extract the fundamental frequency
% out of the input signal.
% 2-  Zeros_finding_Rowdata_final.m: this function will be used to extract
% the zero crossing of the input signal and then internally takes only the
% zeros points which fulfill the 'tol' argument mentioned before. so for
% example if you have an input signal where too much harmonics are inside
% and some missleading zero crossing in there is function will try to
% remove those unwanted zeros points.
% Author: Aubai Al Khatib 
% 15.01.2014 Created
% 28.01.2013 Modified
Delta_Time_Orginal = time(2)- time(1);
F_Sample = abs(1/Delta_Time_Orginal);
%------------------------------------------------------------------------------------------------------
%---------------------- The Harmonics Spectrum of the complet input signal ----------------------------
%------------------------------------------------------------------------------------------------------
Tp = length(time);
T_S = abs(Tp*Delta_Time_Orginal);
[nyquist_number_I1,Bezugsfrequenz_I1,Signal_amp_I1] = Fast_Harmonics (Signal,F_Sample,T_S);
%plot([0:nyquist_number_I1]*Bezugsfrequenz_I1,Signal_amp_I1);ylabel('amplitude Currents [A]');title('Harmonic spectrom');hold on
F_Spec = [0:nyquist_number_I1]*Bezugsfrequenz_I1; %#ok<*NBRAK>
Fn = round(F_Spec(Signal_amp_I1 == max(Signal_amp_I1)));
[~,Zeros_f,Zeros_f_indx_final] = Zeros_finding_Rowdata_final(Signal,time,Fn,tol);
%[Zeros,k] = Zeros_finding(Signal,time);
Zeros = Zeros_f;
k = Zeros_f_indx_final;
f = 1./(2*diff(Zeros));
f_time = Zeros(1:length(f));
Test_indx = k(1)+2;
if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
  signal_test = Signal;time_test = time;
  signal_test(k) = 0;time_test(k) = Zeros;
  if 0 == 1
      Signal_Zeros = Signal(k(2):k(end)-1);
      time_Zeros = time(k(2):k(end)-1);
      Signal_Final = Signal_Zeros;
      Signal_Final(1) = 0; Signal_Final(end) = 0;
      Signal_Zeros = Signal_Final;
      time_Final = time_Zeros;
      time_Final(1) = Zeros(2);time_Final(end) = Zeros(end);
      f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f(2:end), diff(k(2:end)),'uniformoutput',0));
      time_f_interp = time_Final;
  else
      Signal_Zeros = signal_test(k(2):k(end)-1);
      time_Zeros = time_test(k(2):k(end)-1);
      time_Final = time_Zeros;
      f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f(2:end), diff(k(2:end)),'uniformoutput',0));
      time_f_interp = time_Final;
  end
else
  signal_test = Signal;time_test = time;
  signal_test(k) = 0;time_test(k) = Zeros;
  if 0 == 1
      Signal_Zeros = Signal(k(1):k(end)-1);
      time_Zeros = time(k(1):k(end)-1);
      Signal_Final = Signal_Zeros;
      Signal_Final(1) = 0; Signal_Final(end) = 0;
      Signal_Zeros = Signal_Final;
      time_Final = time_Zeros;
      time_Final(1) = Zeros(1);time_Final(end) = Zeros(end);
      f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
      time_f_interp = time_Final;
  else
      Signal_Zeros = signal_test(k(1):k(end)-1);
      time_Zeros = time_test(k(1):k(end)-1);
      time_Final = time_Zeros;
      f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
      time_f_interp = time_Final;
  end
end
%f_interp = cell2mat(arrayfun(@(x,nx) repmat(x,1,nx), f, diff(k),'uniformoutput',0));
%time_f_interp = time_Final;
%time_f_interp = time(k(1):k(end)-1);
time_Zeros_1 = time_Zeros;
time_Zeros = time(k(1):k(end)-1);%time_f_interp;
if option
  W = 2*pi*f_interp;
  Wt = cumtrapzt(time_Zeros,W);
  Wt = mod(Wt,2*pi);
  if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
      Start_Indx = k(2);
      End_Indx = k(end)-1;
      Sine = -sin(Wt);
      Cose = -cos(Wt);
  else
      Start_Indx = k(1);
      End_Indx = k(end)-1;
      Sine = sin(Wt);
      Cose = cos(Wt);
  end
  Signal_Limits = [Start_Indx, End_Indx];
  PLL_time = time_f_interp;
else
  Zeros_one_period = Zeros(1:2:end);
  D_Zeros = kron(Zeros_one_period,[1 1])';
  Test = 2:2:length(D_Zeros);
  D_Zeros(Test) = D_Zeros(Test)+0.00000001;
  A_max = 0;
  A_min = 2*pi;
  A = [A_min; A_max];%*(pi/180);
  R_Z_G = repmat(A,length(D_Zeros),1);
  R_Z_G_1 = R_Z_G(1:length(D_Zeros));
  %R_Z_G_2 = R_Z_G(1:length(Zeros));
  R_Z_G_interp = interp1(D_Zeros,R_Z_G_1,time,'linear');
  %R_Z_G_interp = interp1(R_Z_G_1,time,'linear');
  R_Z_G_interp_test = R_Z_G_interp(~isnan(R_Z_G_interp));
  Wt = R_Z_G_interp_test;
  Test_indx = k(1)+2;
  if Zeros(1) < 1/(Fn*2) && Signal(Test_indx) < 0
      Sine = -sin(R_Z_G_interp_test);
      Cose = -cos(R_Z_G_interp_test);
  else
      Sine = sin(R_Z_G_interp_test);
      Cose = cos(R_Z_G_interp_test);
  end
  Time_test = time(isnan(R_Z_G_interp) == 0);
  PLL_time = Time_test;
  Start_Indx = find(time == Time_test(1));
  End_Indx = find(time == Time_test(end));
  Signal_Limits = [Start_Indx, End_Indx];
  Signal_Zeros = Signal(isnan(R_Z_G_interp) == 0);
end
plots = 0;
if plots
  if 0 == 1
      figure(); %#ok<*UNRCH>
      if 1 == 1
          x_min = 1;
          x_max = x_min + .01;%time(end);%0.1;
      elseif 1 == 0
          x_min = time_Zeros(end)-6.5;
          x_max =  x_min + 0.1;%time(end);%0.1;
      elseif 0 == 1
          x_min = time_Zeros(1);
          x_max = x_min + 0.1;%time(end);%0.1;
      else
          x_min = time_Zeros(1);
          x_max = time_Zeros(end);%time(end);%0.1;
      end
      subplot(4,1,1);
      [AX,ha2,ha3] = plotyy(time,Signal,time_Zeros,Sine);grid on;
      set(AX(1),'xlim',[x_min x_max])
      set(AX(2),'xlim',[x_min x_max])
      set(AX(1),'ylim',[-max(Signal) max(Signal)])
      subplot(4,1,2);
      plot(time_Zeros,f_interp);grid on;xlim([x_min x_max]);
      subplot(4,1,3);
      plot(time_Zeros,Wt);grid on;xlim([x_min x_max]);
      subplot(4,1,4);
      plot(time_Zeros,Sine);grid on;xlim([x_min x_max]);hold on
      plot(time_Zeros,Cose,'r');hold off
  else
      figure(); %#ok<*UNRCH>
      if 1 == 1
          x_min = 80;
          x_max = x_min + 1;%time(end);%0.1;
      elseif 1 == 0
          x_min = time_Zeros(end)-6.5;
          x_max =  x_min + 0.5;%time(end);%0.1;
      elseif 1 == 0
          x_min = time_Zeros(1);
          x_max = x_min + 0.1;%time(end);%0.1;
      elseif 1 == 0
          x_min = time_Zeros(end)-5;
          x_max =  x_min + 0.1;%time(end);%0.1;
      else
          x_min = time_Zeros(1);
          x_max = time_Zeros(end);%time(end);%0.1;
      end
      subplot(4,1,1);
      [AX,ha2,ha3] = plotyy(time,Signal,time_Final,Sine);grid on;
      set(AX(1),'xlim',[x_min x_max])
      set(AX(2),'xlim',[x_min x_max])
      set(AX(1),'ylim',[-max(Signal) max(Signal)])
      subplot(4,1,2);
      plot(time_Final,f_interp);grid on;xlim([x_min x_max]);
      subplot(4,1,3);
      plot(time_Final,Wt);grid on;xlim([x_min x_max]);
      subplot(4,1,4);
      plot(time_Final,Sine);grid on;xlim([x_min x_max]);hold on
      plot(time_Final,Cose,'r');hold off
  end
end
end
end

Answer 1

[Daniel kiracofe]

the short answer is that you cannot just multiply f by t to get what you want. You need to integrate f with respect to time to get theta. i.e. instead of sin(2 * pi * f * t) you want sin(2 * pi * cumtrapz(t, f)).

this short tutorial I wrote goes into more detail: http://mechanicalvibration.com/Generating_signal_from_freq.html

[Aubai]

thx very much for the information it is helpfull and easy to understand :)

as i can see you have a good idea about the subject so if you try to follow my code u will find that i have integrated my frequency and my time in order to obtain my wished generated sine wave those steps are done in the line code:

if true
  % code
W = 2*pi*f_interp;
Wt = cumtrapzt(time_Zeros,W);
Wt = mod(Wt,2*pi);
Sine = sin(Wt);
end

note: cumtrapzt is the same build in cumtrapz function used by matlab i only did some changes in order to get a non-delayed effect integrational results (not the subject of this question).

so my question to you is: how can i generate a correct sine wave in which my time changes are not constant which means the sampling rate is changable and not constant the complete time period do you have an idea?

Thx for you help and understanding

[Daniel kiracofe]

I think the easiest way is to just use interp1(). i.e.

Wtnew = interp1( time_Zeros, W, t); Sine = sin(Wtnew);

in this way, Sine will have a sample at every value of the vector t (which could have any sampling rate you want, even a non-constant one). you could also use cubic spline if linear interpolation is not good enough.

note the mod(Wt, 2*pi) is not necessary. it doesn't hurt, but it also is not doing anything for you.